LCOV - code coverage report
Current view: | top level - chemistry/mozart - mo_airplane.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 19 | 155 | 12.3 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 2 | 2 | 100.0 %

              Line data    Source code
           1             : module mo_airplane
           2             :   !--------------------------------------------------------------------
           3             :   !     ... Airplane insitu emission sources
           4             :   !--------------------------------------------------------------------
           5             :
           6             :   use shr_kind_mod,     only : r8 => shr_kind_r8
           7             :   use cam_abortutils,   only : endrun
           8             :   use pio,              only : pio_inq_dimid, pio_inq_dimlen, pio_get_var, &
           9             :                                file_desc_t, var_desc_t, pio_inq_vardimid, pio_inq_varndims, pio_nowrite, &
          10             :                                pio_inq_varid, pio_closefile
          11             :
          12             :   use cam_pio_utils,    only : cam_pio_openfile
          13             :   use cam_logfile,      only : iulog
          14             :   implicit none
          15             :
          16             :   private
          17             :
          18             :   save
          19             :
          20             :   real(r8), allocatable :: &
          21             :        pno(:,:,:), &
          22             :        pco(:,:,:), &
          23             :        air_altitude(:)
          24             :   public :: airpl_set, airpl_src
          25             :
          26             :
          27             :
          28             :   logical :: has_airpl_src = .false.
          29             :
          30             : contains
          31             :
          32     1489176 :   subroutine airpl_set( lchnk, ncol, no_ndx, co_ndx, xno_ndx, cldtop, zint_abs, extfrc)
          33             :     use ppgrid,       only : pver
          34             :     use cam_history,  only : outfld
          35             :
          36             :     implicit none
          37             :     integer, intent(in) :: lchnk, ncol, no_ndx, co_ndx, xno_ndx
          38             :     real(r8), intent(in) :: cldtop(:), zint_abs(:,:)
          39             :     real(r8), intent(inout) :: extfrc(:,:,:)
          40             :
          41             :
          42             : ! Local Variables
          43     1489176 :     real(r8), dimension(ncol,pver) :: no_air, co_air
          44             :     real(r8) :: ztab_top, ztab_bot, zdel, zdeli, frac, zlev_top, zlev_bot
          45             :     integer :: nlev
          46             :     integer :: cldind, kk, i, k
          47             :
          48  2314006344 :     no_air(:,:) = 0._r8
          49  2314006344 :     co_air(:,:) = 0._r8
          50             :
          51     1489176 :     if(has_airpl_src) then
          52             :        !---------------------------------------------------------------------
          53             :        !     ... Add the airplane emissions; must do vertical interpolation
          54             :        !---------------------------------------------------------------------
          55           0 :        ztab_top = maxval( air_altitude )
          56           0 :        ztab_bot = minval( air_altitude )
          57           0 :        nlev     = size(air_altitude) - 1
          58             :
          59             :        !---------------------------------------------------------------------
          60             :        !     ... add the airplane emissions; must do vertical interpolation
          61             :        !         Note: the interpolation code is conserving and assumes the
          62             :        !               aircraft emission vertical grid is uniform with a
          63             :        !               one kilometer spacing
          64             :        !---------------------------------------------------------------------
          65           0 :        level_loop : do k = 1,pver
          66           0 :           long_loop : do i = 1,ncol
          67           0 :              zlev_top = zint_abs(i,k)                 ! altitude at top of model level (km)
          68           0 :              zlev_bot = zint_abs(i,k+1)               ! altitude at bottom of model level (km)
          69           0 :              zdel = (zlev_top - zlev_bot) * 1.e5_r8  ! model level thickness (cm)
          70           0 :              zdeli = 1._r8/zdel
          71           0 :              if( zlev_bot <= ztab_top .and. zlev_top >= ztab_bot ) then
          72           0 :                 do kk = 1,nlev
          73           0 :                    if( zlev_bot <= air_altitude(kk+1) .and. zlev_top >= air_altitude(kk) ) then
          74             :                       frac = (min( zlev_top, air_altitude(kk+1) ) - max( zlev_bot, air_altitude(kk) )) &
          75           0 :                            /(air_altitude(kk+1) - air_altitude(kk)) ! *del_alti(kk)
          76           0 :                       if( no_ndx > 0 ) then
          77           0 :                          extfrc(i,k,no_ndx) = extfrc(i,k,no_ndx) + frac * pno(i,kk,lchnk) * zdeli
          78           0 :                          no_air(i,k) = frac * pno(i,kk,lchnk) * zdeli
          79             :                       end if
          80           0 :                       if( xno_ndx > 0 ) then
          81           0 :                          extfrc(i,k,xno_ndx) = extfrc(i,k,xno_ndx) + frac * pno(i,kk,lchnk) * zdeli
          82             :                       end if
          83           0 :                       if( co_ndx > 0 ) then
          84           0 :                          extfrc(i,k,co_ndx) = extfrc(i,k,co_ndx) + frac * pco(i,kk,lchnk) * zdeli
          85           0 :                          co_air(i,k) = frac * pco(i,kk,lchnk) * zdeli
          86             :                       end if
          87             :                    end if
          88             :                 end do
          89             :              end if
          90           0 :              if( k == pver ) then
          91           0 :                 do kk = 1,nlev
          92           0 :                    if( zlev_bot > air_altitude(kk) ) then
          93           0 :                       frac = (min( zlev_bot, air_altitude(kk+1) ) - air_altitude(kk)) &
          94           0 :                            /(air_altitude(kk+1) - air_altitude(kk)) ! *del_alti(kk)
          95           0 :                       if( no_ndx > 0 ) then
          96           0 :                          extfrc(i,k,no_ndx) = extfrc(i,k,no_ndx) + frac * pno(i,kk,lchnk) * zdeli
          97           0 :                          no_air(i,k) = frac * pno(i,kk,lchnk) * zdeli
          98             :                       end if
          99           0 :                       if( xno_ndx > 0 ) then
         100           0 :                          extfrc(i,k,xno_ndx) = extfrc(i,k,xno_ndx) + frac * pno(i,kk,lchnk) * zdeli
         101             :                       end if
         102           0 :                       if( co_ndx > 0 ) then
         103           0 :                          extfrc(i,k,co_ndx) = extfrc(i,k,co_ndx) + frac * pco(i,kk,lchnk) * zdeli
         104           0 :                          co_air(i,k) = frac * pco(i,kk,lchnk) * zdeli
         105             :                       end if
         106             :                    else
         107             :                       exit
         108             :                    end if
         109             :                 end do
         110             :              end if
         111             :           end do long_loop
         112             :        end do level_loop
         113             :
         114             :     end if
         115     1489176 :     call outfld( 'NO_Aircraft',  no_air(:ncol,:), ncol, lchnk )
         116     1489176 :     call outfld( 'CO_Aircraft',  co_air(:ncol,:), ncol, lchnk )
         117             :
         118     1489176 :   end subroutine airpl_set
         119             :
         120             :
         121        1536 :   subroutine airpl_src( airpl_emis_file )
         122             :     !-----------------------------------------------------------------------
         123             :     !   ... Initialize airplane emissions
         124             :     !       Note: The emissions are read in in units of molecules/cm**2/s
         125             :     !             on a vertically resolved grid.
         126             :     !             Conversion to units of molec/cm**3/s is done in SETEXT
         127             :     !-----------------------------------------------------------------------
         128     1489176 :     use spmd_utils,    only : masterproc
         129             :     use interpolate_data, only : lininterp_init, lininterp, lininterp_finish, &
         130             :          interp_type
         131             :     use chem_mods,     only : adv_mass
         132             :     use ioFileMod,     only : getfil
         133             :     use mo_chem_utls,  only : get_spc_ndx, get_extfrc_ndx
         134             :     use phys_grid,     only : get_ncols_p, get_rlat_all_p, get_rlon_all_p
         135             :     use ppgrid,        only : begchunk, endchunk, pcols
         136             :     use mo_constants,  only : pi, d2r, rearth
         137             :     use gmean_mod,     only : gmean
         138             :     implicit none
         139             :
         140             :     !-----------------------------------------------------------------------
         141             :     !   ... Dummy args
         142             :     !-----------------------------------------------------------------------
         143             :     character(len=*), intent(in) :: airpl_emis_file
         144             :
         145             :     !-----------------------------------------------------------------------
         146             :     !   ... Local variables
         147             :     !-----------------------------------------------------------------------
         148             :     real(r8), parameter :: msq2cmsq = 1.e4_r8, zero=0._r8, twopi=2._r8*pi
         149             :     integer  :: k
         150             :     integer  :: nlat, nlon, nlev, ndims
         151             :     integer  :: ierr
         152             :     type(file_desc_t) :: piofile
         153             :     type(var_desc_t) :: vid
         154             :     integer  :: dimid_lat, dimid_lon, dimid_lev
         155             :     integer  :: dimid(3)
         156        1536 :     real(r8), allocatable :: lat(:), lon(:)
         157        1536 :     real(r8), allocatable :: pno_in(:,:,:), pco_in(:,:,:)
         158             :     real(r8) :: total(2), tmp
         159             :     real(r8) :: factor
         160             :     character(len=256) :: locfn
         161             :     integer :: co_ndx, no_ndx
         162             :     type(interp_type) :: lon_wgts, lat_wgts
         163             :     real(r8) :: to_lats(pcols), to_lons(pcols)
         164             :     integer :: ncols, c
         165             :
         166        1536 :     co_ndx = get_extfrc_ndx('CO')
         167        1536 :     no_ndx = get_extfrc_ndx('NO')
         168             :
         169        1536 :     if ( co_ndx < 0 .and. no_ndx < 0 ) then
         170        1536 :        if( masterproc ) then
         171           2 :           write(iulog,*) 'airpl_src: NO and CO do not have external source --> no aircraft sources will be applied'
         172             :        endif
         173        1536 :        return
         174             :     endif
         175             :
         176           0 :     if ( len_trim(airpl_emis_file) == 0 ) then
         177             :        return
         178             :     endif
         179             :
         180           0 :     has_airpl_src = .true.
         181             :
         182           0 :     co_ndx = get_spc_ndx('CO')
         183           0 :     no_ndx = get_spc_ndx('NO')
         184             :
         185             :
         186             :     !-----------------------------------------------------------------------
         187             :     !   ... Open NetCDF file
         188             :     !-----------------------------------------------------------------------
         189           0 :     call getfil (airpl_emis_file, locfn, 0)
         190           0 :     call cam_pio_openfile (piofile, trim(locfn), PIO_NOWRITE)
         191             :
         192             :     !-----------------------------------------------------------------------
         193             :     !       ... Get grid dimensions from file
         194             :     !-----------------------------------------------------------------------
         195           0 :     ierr = pio_inq_dimid( piofile, 'lat', dimid_lat )
         196           0 :     ierr = pio_inq_dimlen( piofile, dimid_lat, nlat )
         197           0 :     allocate( lat(nlat), stat=ierr )
         198           0 :     if( ierr /= 0 ) then
         199           0 :        write(iulog,*) 'airpl_src: lat allocation error = ',ierr
         200           0 :        call endrun
         201             :     end if
         202           0 :     ierr = pio_inq_varid( piofile, 'lat', vid )
         203           0 :     ierr = pio_get_var( piofile, vid, lat )
         204           0 :     lat(:nlat) = lat(:nlat) * d2r
         205             :
         206           0 :     ierr = pio_inq_dimid( piofile, 'lon', dimid_lon )
         207           0 :     ierr = pio_inq_dimlen( piofile, dimid_lon, nlon )
         208           0 :     allocate( lon(nlon), stat=ierr )
         209           0 :     if( ierr /= 0 ) then
         210           0 :        write(iulog,*) 'airpl_src: lon allocation error = ',ierr
         211           0 :        call endrun
         212             :     end if
         213           0 :     ierr = pio_inq_varid( piofile, 'lon', vid )
         214           0 :     ierr = pio_get_var( piofile, vid, lon )
         215           0 :     lon(:nlon) = lon(:nlon) * d2r
         216             :
         217           0 :     ierr = pio_inq_dimid( piofile, 'altitude', dimid_lev )
         218           0 :     ierr = pio_inq_dimlen( piofile, dimid_lev, nlev )
         219           0 :     allocate( air_altitude(nlev+1), stat=ierr )
         220           0 :     if( ierr /= 0 ) then
         221           0 :        write(iulog,*) 'airpl_src: air_altitude allocation error = ',ierr
         222           0 :        call endrun
         223             :     end if
         224           0 :     ierr = pio_inq_varid( piofile, 'altitude', vid )
         225           0 :     ierr = pio_get_var( piofile, vid, air_altitude(1:nlev) )
         226           0 :     air_altitude(nlev+1) = air_altitude(nlev) + (air_altitude(nlev) - air_altitude(nlev-1))
         227             :
         228             :     !-----------------------------------------------------------------------
         229             :     !       ... Set up regridding
         230             :     !-----------------------------------------------------------------------
         231             :
         232           0 :     allocate( pno_in(nlon,nlat,nlev), stat=ierr )
         233           0 :     if( ierr /= 0 ) then
         234           0 :        write(iulog,*) 'airpl_src: pno_in allocation error = ',ierr
         235           0 :        call endrun
         236             :     end if
         237           0 :     allocate( pco_in(nlon,nlat,nlev), stat=ierr )
         238           0 :     if( ierr /= 0 ) then
         239           0 :        write(iulog,*) 'airpl_src: pco_in allocation error = ',ierr
         240           0 :        call endrun
         241             :     end if
         242           0 :     allocate(pno(pcols,nlev,begchunk:endchunk), stat=ierr )
         243           0 :     if( ierr /= 0 ) then
         244           0 :        write(iulog,*) 'airpl_src: pno allocation error = ',ierr
         245           0 :        call endrun
         246             :     end if
         247           0 :     allocate( pco(pcols,nlev,begchunk:endchunk), stat=ierr )
         248           0 :     if( ierr /= 0 ) then
         249           0 :        write(iulog,*) 'airpl_src: pco allocation error = ',ierr
         250           0 :        call endrun
         251             :     end if
         252             :
         253             :     !-----------------------------------------------------------------------
         254             :     !   ... Read emissions
         255             :     !-----------------------------------------------------------------------
         256           0 :     ierr = pio_inq_varid( piofile, 'nox', vid )
         257           0 :     ierr = pio_inq_varndims( piofile, vid, ndims )
         258           0 :     if( ndims /= 3 ) then
         259           0 :        write(iulog,*) 'airpl_src: variable nox has ndims = ',ndims,', expecting 3'
         260           0 :        call endrun
         261             :     end if
         262           0 :     ierr = pio_inq_vardimid( piofile, vid, dimid )
         263           0 :     if( dimid(1) /= dimid_lon  .or. dimid(2) /= dimid_lat .or.  dimid(3) /= dimid_lev ) then
         264           0 :        write(iulog,*) 'airpl_src: Dimensions in wrong order for variable nox'
         265           0 :        write(iulog,*) '...      Expecting (lon, lat, lev)'
         266           0 :        call endrun
         267             :     end if
         268             :     ierr = pio_get_var( piofile, vid, &
         269             :          (/ 1, 1, 1/), &                    ! start
         270             :          (/ nlon, nlat, nlev /), &   ! count
         271           0 :          pno_in )
         272             :
         273           0 :     ierr = pio_inq_varid( piofile, 'co', vid )
         274           0 :     ierr = pio_inq_varndims( piofile, vid, ndims )
         275             :
         276           0 :     if( ndims /= 3 ) then
         277           0 :        write(iulog,*) 'READ_SFLX: variable co has ndims = ',ndims,', expecting 3'
         278           0 :        call endrun
         279             :     end if
         280           0 :     ierr = pio_inq_vardimid( piofile, vid, dimid )
         281           0 :     if( dimid(1) /= dimid_lon .or. dimid(2) /= dimid_lat .or. dimid(3) /= dimid_lev ) then
         282           0 :        write(iulog,*) 'airpl_src: Dimensions in wrong order for variable co'
         283           0 :        write(iulog,*) '...      Expecting (lon, lat, lev)'
         284           0 :        call endrun
         285             :     end if
         286             :     ierr = pio_get_var( piofile, vid, &
         287             :          (/ 1, 1, 1/), &                    ! start
         288             :          (/ nlon, nlat, nlev /), &   ! count
         289           0 :          pco_in )
         290           0 :     call pio_closefile( piofile )
         291             :
         292             :     !-----------------------------------------------------------------------
         293             :     !   ... Regrid emissions
         294             :     !-----------------------------------------------------------------------
         295           0 :     do c=begchunk,endchunk
         296           0 :        ncols = get_ncols_p(c)
         297           0 :        call get_rlat_all_p(c, pcols, to_lats)
         298           0 :        call get_rlon_all_p(c, pcols, to_lons)
         299           0 :        call lininterp_init(lon, nlon, to_lons, ncols, 2, lon_wgts, zero, twopi)
         300           0 :        call lininterp_init(lat, nlat, to_lats, ncols, 1, lat_wgts)
         301             :
         302           0 :        do k = 1,nlev
         303           0 :           call lininterp(pno_in(:,:,k), nlon, nlat, pno(:,k,c), ncols, lon_wgts, lat_wgts)
         304           0 :           call lininterp(pco_in(:,:,k), nlon, nlat, pco(:,k,c), ncols, lon_wgts, lat_wgts)
         305             :        enddo
         306           0 :        call lininterp_finish(lon_wgts)
         307           0 :        call lininterp_finish(lat_wgts)
         308             :     enddo
         309             :
         310           0 :     deallocate( pno_in, pco_in, lon, lat, stat=ierr )
         311           0 :     if( ierr /= 0 ) then
         312           0 :        write(iulog,*) 'airpl_src: Failed to deallocate pno_in,pco_in; ierr = ',ierr
         313           0 :        call endrun
         314             :     end if
         315             :     !-----------------------------------------------------------------------
         316             :     !       ... Get global emission from this source
         317             :     !-----------------------------------------------------------------------
         318           0 :     total = zero
         319           0 :     do k=1,nlev
         320           0 :        call gmean(pno(:,k,:), tmp)
         321           0 :        total(1)= total(1)+tmp
         322           0 :        call gmean(pco(:,k,:), tmp)
         323           0 :        total(2)= total(2)+tmp
         324             :     end do
         325             :
         326           0 :     if(masterproc) then
         327             :
         328             :        factor = 86400._r8 * 365._r8 &   ! sec / year
         329             :             / 6.022e23_r8 &           ! molec / mole
         330             :             * 1.e-12_r8   &            ! Tg / g
         331             :             * msq2cmsq    &             ! meters**2 to cm**2
         332           0 :             * 4._r8*pi*rearth*rearth      ! global mean to global total
         333             :
         334           0 :        write(iulog,*) 'airpl_src: nlev = ',nlev
         335             :        !-----------------------------------------------------------------------
         336             :        !       ... Convert totals from molec cm^-2 s^-1 to Tg y^-1
         337             :        !-----------------------------------------------------------------------
         338           0 :        if (no_ndx .gt. 0) then
         339           0 :           total(1) = total(1) * adv_mass(no_ndx) * factor
         340           0 :           write(iulog,'('' airpl_src Aircraft emissions: '',a6,'' = '',f10.3,1X,a6)') 'NO',total(1),'TgN/y'
         341             :        endif
         342           0 :        if (co_ndx .gt. 0) then
         343           0 :           total(2) = total(2) * adv_mass(co_ndx) * factor
         344           0 :           write(iulog,'('' airpl_src Aircraft emissions: '',a6,'' = '',f10.3,1X,a6)') 'CO',total(2),'Tg/y'
         345             :        endif
         346             :     end if
         347             :
         348        3072 :   end subroutine airpl_src
         349             :
         350             : end module mo_airplane

Generated by: LCOV version 1.14