LCOV - code coverage report
Current view: | top level - chemistry/mozart - mo_srf_emissions.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 133 | 171 | 77.8 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 3 | 5 | 60.0 %

              Line data    Source code
           1             : module mo_srf_emissions
           2             :   !---------------------------------------------------------------
           3             :   !     ... surface emissions module
           4             :   !---------------------------------------------------------------
           5             :
           6             :   use shr_kind_mod,  only : r8 => shr_kind_r8
           7             :   use chem_mods,     only : gas_pcnst
           8             :   use spmd_utils,    only : masterproc
           9             :   use cam_abortutils,only : endrun
          10             :   use ioFileMod,     only : getfil
          11             :   use cam_logfile,   only : iulog
          12             :   use tracer_data,   only : trfld,trfile
          13             :
          14             :   implicit none
          15             :
          16             :   type :: emission
          17             :      integer           :: spc_ndx
          18             :      real(r8)          :: mw
          19             :      real(r8)          :: scalefactor
          20             :      character(len=256):: filename
          21             :      character(len=16) :: species
          22             :      character(len=8)  :: units
          23             :      integer                   :: nsectors
          24             :      character(len=32),pointer :: sectors(:)
          25             :      type(trfld), pointer      :: fields(:)
          26             :      type(trfile)              :: file
          27             :   end type emission
          28             :
          29             :   private
          30             :
          31             :   public  :: srf_emissions_inti, set_srf_emissions, set_srf_emissions_time
          32             :
          33             :   logical, public, protected :: has_emis(gas_pcnst) = .false.
          34             :
          35             :   real(r8), parameter :: amufac = 1.65979e-23_r8         ! 1.e4* kg / amu
          36             :   type(emission), allocatable :: emissions(:)
          37             :   integer                     :: n_emis_files
          38             :   integer :: c10h16_ndx, isop_ndx
          39             :
          40             : contains
          41             :
          42        1536 :   subroutine srf_emissions_inti( srf_emis_specifier, emis_type_in, emis_cycle_yr, emis_fixed_ymd, emis_fixed_tod )
          43             :
          44             :     !-----------------------------------------------------------------------
          45             :     !   ... initialize the surface emissions
          46             :     !-----------------------------------------------------------------------
          47             :
          48             :     use chem_mods,        only : adv_mass
          49             :     use mo_chem_utls,     only : get_spc_ndx
          50             :     use tracer_data,      only : trcdata_init
          51             :     use cam_pio_utils,    only : cam_pio_openfile
          52             :     use pio,              only : pio_inquire, pio_nowrite, pio_closefile, pio_inq_varndims
          53             :     use pio,              only : pio_inq_varname, pio_inq_vardimid, pio_inq_dimid
          54             :     use pio,              only : file_desc_t, pio_get_att, PIO_NOERR, PIO_GLOBAL
          55             :     use pio,              only : pio_seterrorhandling, PIO_BCAST_ERROR,PIO_INTERNAL_ERROR
          56             :     use chem_surfvals,    only : flbc_list
          57             :     use string_utils,     only : GLC
          58             :     use m_MergeSorts,     only : IndexSort
          59             :
          60             :     implicit none
          61             :
          62             :     !-----------------------------------------------------------------------
          63             :     !   ... dummy arguments
          64             :     !-----------------------------------------------------------------------
          65             :     character(len=*), intent(in) :: srf_emis_specifier(:)
          66             :     character(len=*), intent(in) :: emis_type_in
          67             :     integer,          intent(in) :: emis_cycle_yr
          68             :     integer,          intent(in) :: emis_fixed_ymd
          69             :     integer,          intent(in) :: emis_fixed_tod
          70             :
          71             :     !-----------------------------------------------------------------------
          72             :     !   ... local variables
          73             :     !-----------------------------------------------------------------------
          74             :     integer  :: astat
          75             :     integer  :: j, l, m, n, i, nn                     ! Indices
          76             :     character(len=16)  :: spc_name
          77             :     character(len=256) :: filename
          78             :
          79        3072 :     character(len=16)  :: emis_species(size(srf_emis_specifier))
          80        3072 :     character(len=256) :: emis_filenam(size(srf_emis_specifier))
          81        3072 :     integer  :: emis_indexes(size(srf_emis_specifier))
          82        3072 :     integer  :: indx(size(srf_emis_specifier))
          83        3072 :     real(r8) :: emis_scalefactor(size(srf_emis_specifier))
          84             :
          85             :     integer :: vid, nvars, isec, num_dims_emis
          86             :     integer :: vndims
          87        1536 :     logical, allocatable :: is_sector(:)
          88             :     type(file_desc_t) :: ncid
          89             :     character(len=32)  :: varname
          90             :     character(len=256) :: locfn
          91             :     integer :: ierr
          92             :     character(len=1), parameter :: filelist = ''
          93             :     character(len=1), parameter :: datapath = ''
          94             :     logical         , parameter :: rmv_file = .false.
          95             :     logical :: unstructured
          96             :     character(len=32) :: emis_type = ' '
          97             :     character(len=80) :: file_interp_type = ' '
          98             :     character(len=256) :: tmp_string = ' '
          99             :     character(len=32) :: xchr = ' '
         100             :     real(r8) :: xdbl
         101             :     integer :: time_dimid, ncol_dimid
         102        1536 :     integer, allocatable :: dimids(:)
         103             :
         104        1536 :     has_emis(:) = .false.
         105        1536 :     nn = 0
         106      155136 :     indx(:) = 0
         107             :
         108       47616 :     count_emis: do n=1,size(srf_emis_specifier)
         109       47616 :        if ( len_trim(srf_emis_specifier(n) ) == 0 ) then
         110             :           exit count_emis
         111             :        endif
         112             :
         113       46080 :        i = scan(srf_emis_specifier(n),'->')
         114       46080 :        spc_name = trim(adjustl(srf_emis_specifier(n)(:i-1)))
         115             :
         116             :        ! need to parse out scalefactor ...
         117       46080 :        tmp_string = adjustl(srf_emis_specifier(n)(i+2:))
         118       46080 :        j = scan( tmp_string, '*' )
         119       46080 :        if (j>0) then
         120       18432 :           xchr = tmp_string(1:j-1) ! get the multipler (left of the '*')
         121       18432 :           read( xchr, * ) xdbl   ! convert the string to a real
         122       18432 :           tmp_string = adjustl(tmp_string(j+1:)) ! get the filepath name (right of the '*')
         123             :        else
         124       27648 :           xdbl = 1._r8
         125             :        endif
         126       46080 :        filename = trim(tmp_string)
         127             :
         128       46080 :        m = get_spc_ndx(spc_name)
         129             :
         130       46080 :        if (m > 0) then
         131       46080 :           has_emis(m) = .true.
         132             :        else
         133           0 :           write(iulog,*) 'srf_emis_inti: spc_name ',spc_name,' is not included in the simulation'
         134           0 :           call endrun('srf_emis_inti: invalid surface emission specification')
         135             :        endif
         136             :
         137     1935360 :        if (any( flbc_list == spc_name )) then
         138             :           call endrun('srf_emis_inti: ERROR -- cannot specify both fixed LBC ' &
         139           0 :                     //'and emissions for the same species: '//trim(spc_name))
         140             :        endif
         141             :
         142       46080 :        nn = nn+1
         143       46080 :        emis_species(nn) = spc_name
         144       46080 :        emis_filenam(nn) = filename
         145       46080 :        emis_indexes(nn) = m
         146       46080 :        emis_scalefactor(nn) = xdbl
         147             :
         148       47616 :        indx(n)=n
         149             :
         150             :     enddo count_emis
         151             :
         152        1536 :     n_emis_files = nn
         153             :
         154        1536 :     if (masterproc) write(iulog,*) 'srf_emis_inti: n_emis_files = ',n_emis_files
         155             :
         156       56832 :     allocate( emissions(n_emis_files), stat=astat )
         157        1536 :     if( astat/= 0 ) then
         158           0 :        write(iulog,*) 'srf_emis_inti: failed to allocate emissions array; error = ',astat
         159           0 :        call endrun('srf_emis_inti: failed to allocate emissions array')
         160             :     end if
         161             :
         162             :     !-----------------------------------------------------------------------
         163             :     ! Sort the input files so that the emissions sources are summed in the
         164             :     ! same order regardless of the order of the input files in the namelist
         165             :     !-----------------------------------------------------------------------
         166        1536 :     if (n_emis_files > 0) then
         167        1536 :       call IndexSort(n_emis_files, indx, emis_filenam)
         168             :     end if
         169             :
         170             :     !-----------------------------------------------------------------------
         171             :     !   ... setup the emission type array
         172             :     !-----------------------------------------------------------------------
         173       47616 :     do m=1,n_emis_files
         174       46080 :        emissions(m)%spc_ndx          = emis_indexes(indx(m))
         175       46080 :        emissions(m)%units            = 'Tg/y'
         176       46080 :        emissions(m)%species          = emis_species(indx(m))
         177       46080 :        emissions(m)%mw               = adv_mass(emis_indexes(indx(m)))                     ! g / mole
         178       46080 :        emissions(m)%filename         = emis_filenam(indx(m))
         179       47616 :        emissions(m)%scalefactor      = emis_scalefactor(indx(m))
         180             :     enddo
         181             :
         182             :     !-----------------------------------------------------------------------
         183             :     ! read emis files to determine number of sectors
         184             :     !-----------------------------------------------------------------------
         185       47616 :     spc_loop: do m = 1, n_emis_files
         186             :
         187       46080 :        emissions(m)%nsectors = 0
         188             :
         189       46080 :        if (masterproc) then
         190          60 :           write(iulog,'(a,i3,a)') 'srf_emissions_inti m: ',m,' init file : '//trim(emissions(m)%filename)
         191             :        endif
         192             :
         193       46080 :        call getfil (emissions(m)%filename, locfn, 0)
         194       46080 :        call cam_pio_openfile ( ncid, trim(locfn), PIO_NOWRITE)
         195       46080 :        ierr = pio_inquire (ncid, nVariables=nvars)
         196             :
         197       46080 :        call pio_seterrorhandling(ncid, PIO_BCAST_ERROR)
         198       46080 :        ierr = pio_inq_dimid( ncid, 'ncol', ncol_dimid )
         199       46080 :        unstructured = ierr==PIO_NOERR
         200       46080 :        call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)
         201             :
         202      138240 :        allocate(is_sector(nvars))
         203      328704 :        is_sector(:) = .false.
         204             :
         205       46080 :        if (unstructured) then
         206       46080 :           ierr = pio_inq_dimid( ncid, 'time', time_dimid )
         207             :        end if
         208             :
         209      328704 :        do vid = 1,nvars
         210             :
         211      282624 :           ierr = pio_inq_varndims (ncid, vid, vndims)
         212             :
         213      282624 :           if (unstructured) then
         214             :              num_dims_emis = 2
         215             :           else
         216           0 :              num_dims_emis = 3
         217             :           endif
         218             :
         219      282624 :           if( vndims < num_dims_emis ) then
         220             :              cycle
         221       52224 :           elseif( vndims > num_dims_emis ) then
         222           0 :              ierr = pio_inq_varname (ncid, vid, varname)
         223           0 :              write(iulog,*) 'srf_emis_inti: Skipping variable ', trim(varname),', ndims = ',vndims, &
         224           0 :                   ' , species=',trim(emissions(m)%species)
         225           0 :              cycle
         226             :           end if
         227             :
         228       98304 :           if (unstructured) then
         229      156672 :              allocate( dimids(vndims) )
         230       52224 :              ierr = pio_inq_vardimid( ncid, vid, dimids )
         231      104448 :              if ( any(dimids(:)==ncol_dimid) .and. any(dimids(:)==time_dimid) ) then
         232       52224 :                 emissions(m)%nsectors = emissions(m)%nsectors+1
         233       52224 :                 is_sector(vid)=.true.
         234             :              endif
         235       52224 :              deallocate(dimids)
         236             :           else
         237           0 :              emissions(m)%nsectors = emissions(m)%nsectors+1
         238           0 :              is_sector(vid)=.true.
         239             :           end if
         240             :
         241             :        enddo
         242             :
         243      138240 :        allocate( emissions(m)%sectors(emissions(m)%nsectors), stat=astat )
         244       46080 :        if( astat/= 0 ) then
         245           0 :          write(iulog,*) 'srf_emis_inti: failed to allocate emissions(m)%sectors array; error = ',astat
         246           0 :          call endrun
         247             :        end if
         248             :
         249       46080 :        isec = 1
         250             :
         251      328704 :        do vid = 1,nvars
         252      328704 :           if( is_sector(vid) ) then
         253       52224 :              ierr = pio_inq_varname(ncid, vid, emissions(m)%sectors(isec))
         254       52224 :              isec = isec+1
         255             :           endif
         256             :        enddo
         257       46080 :        deallocate(is_sector)
         258             :
         259             :        ! Global attribute 'input_method' overrides the srf_emis_type namelist setting on
         260             :        ! a file-by-file basis.  If the emis file does not contain the 'input_method'
         261             :        ! attribute then the srf_emis_type namelist setting is used.
         262       46080 :        call pio_seterrorhandling(ncid, PIO_BCAST_ERROR)
         263       46080 :        ierr = pio_get_att(ncid, PIO_GLOBAL, 'input_method', file_interp_type)
         264       46080 :        call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)
         265       46080 :        if ( ierr == PIO_NOERR) then
         266           0 :           l = GLC(file_interp_type)
         267           0 :           emis_type(1:l) = file_interp_type(1:l)
         268           0 :           emis_type(l+1:) = ' '
         269             :        else
         270       46080 :           emis_type = trim(emis_type_in)
         271             :        endif
         272             :
         273       46080 :        call pio_closefile (ncid)
         274             :
         275      138240 :        allocate(emissions(m)%file%in_pbuf(size(emissions(m)%sectors)))
         276       98304 :        emissions(m)%file%in_pbuf(:) = .false.
         277             :
         278           0 :        call trcdata_init( emissions(m)%sectors, &
         279             :                           emissions(m)%filename, filelist, datapath, &
         280             :                           emissions(m)%fields,  &
         281             :                           emissions(m)%file, &
         282       47616 :                           rmv_file, emis_cycle_yr, emis_fixed_ymd, emis_fixed_tod, trim(emis_type) )
         283             :
         284             :     enddo spc_loop
         285             :
         286        1536 :     c10h16_ndx = get_spc_ndx('C10H16')
         287        1536 :     isop_ndx = get_spc_ndx('ISOP')
         288             :
         289        3072 :   end subroutine srf_emissions_inti
         290             :
         291      741888 :   subroutine set_srf_emissions_time( pbuf2d, state )
         292             :     !-----------------------------------------------------------------------
         293             :     !       ... check serial case for time span
         294             :     !-----------------------------------------------------------------------
         295             :
         296        1536 :     use physics_types,only : physics_state
         297             :     use ppgrid,       only : begchunk, endchunk
         298             :     use tracer_data,  only : advance_trcdata
         299             :     use physics_buffer, only : physics_buffer_desc
         300             :
         301             :     implicit none
         302             :
         303             :     type(physics_state), intent(in):: state(begchunk:endchunk)
         304             :     type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         305             :
         306             :     !-----------------------------------------------------------------------
         307             :     !       ... local variables
         308             :     !-----------------------------------------------------------------------
         309             :     integer :: m
         310             :
         311    11499264 :     do m = 1,n_emis_files
         312    11499264 :        call advance_trcdata( emissions(m)%fields, emissions(m)%file, state, pbuf2d  )
         313             :     end do
         314             :
         315      370944 :   end subroutine set_srf_emissions_time
         316             :
         317             :   ! adds surf flux specified in file to sflx
         318     1489176 :   subroutine set_srf_emissions( lchnk, ncol, sflx )
         319             :     !--------------------------------------------------------
         320             :     !   ... form the surface fluxes for this latitude slice
         321             :     !--------------------------------------------------------
         322             :
         323      370944 :     use mo_constants, only : pi
         324             :     use time_manager, only : get_curr_calday
         325             :     use string_utils, only : to_lower, GLC
         326             :     use phys_grid,    only : get_rlat_all_p, get_rlon_all_p
         327             :
         328             :     implicit none
         329             :
         330             :     !--------------------------------------------------------
         331             :     !   ... Dummy arguments
         332             :     !--------------------------------------------------------
         333             :     integer,  intent(in)  :: ncol                  ! columns in chunk
         334             :     integer,  intent(in)  :: lchnk                 ! chunk index
         335             :     real(r8), intent(out) :: sflx(:,:) ! surface emissions ( kg/m^2/s )
         336             :
         337             :     !--------------------------------------------------------
         338             :     !   ... local variables
         339             :     !--------------------------------------------------------
         340             :     integer  ::  i, m, n
         341             :     real(r8) ::  factor
         342             :     real(r8) ::  dayfrac            ! fration of day in light
         343             :     real(r8) ::  iso_off            ! time iso flux turns off
         344             :     real(r8) ::  iso_on             ! time iso flux turns on
         345             :
         346             :     logical  :: polar_day,polar_night
         347             :     real(r8) :: doy_loc
         348             :     real(r8) :: sunon,sunoff
         349             :     real(r8) :: loc_angle
         350             :     real(r8) :: latitude
         351             :     real(r8) :: declination
         352             :     real(r8) :: tod
         353             :     real(r8) :: calday
         354             :
         355             :     real(r8), parameter :: dayspy = 365._r8
         356             :     real(r8), parameter :: twopi = 2.0_r8 * pi
         357             :     real(r8), parameter :: pid2  = 0.5_r8 * pi
         358             :     real(r8), parameter :: dec_max = 23.45_r8 * pi/180._r8
         359             :
         360     2978352 :     real(r8) :: flux(ncol)
         361             :     real(r8) :: mfactor
         362             :     integer  :: isec
         363             :
         364             :     character(len=12),parameter :: mks_units(4) = (/ "kg/m2/s     ", &
         365             :                                                      "kg/m2/sec   ", &
         366             :                                                      "kg/m^2/s    ", &
         367             :                                                      "kg/m^2/sec  " /)
         368             :     character(len=12) :: units
         369             :
         370     2978352 :     real(r8), dimension(ncol) :: rlats, rlons
         371             :
         372   786284928 :     sflx(:,:) = 0._r8
         373             :
         374             :     !--------------------------------------------------------
         375             :     !   ... set non-zero emissions
         376             :     !--------------------------------------------------------
         377    46164456 :     emis_loop : do m = 1,n_emis_files
         378             :
         379    44675280 :        n = emissions(m)%spc_ndx
         380             :
         381   745973280 :        flux(:) = 0._r8
         382    95307264 :        do isec = 1,emissions(m)%nsectors
         383   890111664 :           flux(:ncol) = flux(:ncol) + emissions(m)%scalefactor*emissions(m)%fields(isec)%data(:ncol,1,lchnk)
         384             :        enddo
         385             :
         386    44675280 :        units = to_lower(trim(emissions(m)%fields(1)%units(:GLC(emissions(m)%fields(1)%units))))
         387             :
         388   224865576 :        if ( any( mks_units(:) == units ) ) then
         389           0 :           sflx(:ncol,n) = sflx(:ncol,n) + flux(:ncol)
         390             :        else
         391    44675280 :           mfactor = amufac * emissions(m)%mw
         392   745973280 :           sflx(:ncol,n) = sflx(:ncol,n) + flux(:ncol) * mfactor
         393             :        endif
         394             :
         395             :     end do emis_loop
         396             :
         397     1489176 :     call get_rlat_all_p( lchnk, ncol, rlats )
         398     1489176 :     call get_rlon_all_p( lchnk, ncol, rlons )
         399             :
         400     1489176 :     calday = get_curr_calday()
         401     1489176 :     doy_loc     = aint( calday )
         402     1489176 :     declination = dec_max * cos((doy_loc - 172._r8)*twopi/dayspy)
         403     1489176 :     tod = (calday - doy_loc) + .5_r8
         404             :
         405    24865776 :     do i = 1,ncol
         406             :        !
         407    23376600 :        polar_day   = .false.
         408    23376600 :        polar_night = .false.
         409             :        !
         410    23376600 :        loc_angle = tod * twopi + rlons(i)
         411    23376600 :        loc_angle = mod( loc_angle,twopi )
         412    23376600 :        latitude =  rlats(i)
         413             :        !
         414             :        !------------------------------------------------------------------
         415             :        !        determine if in polar day or night
         416             :        !        if not in polar day or night then
         417             :        !        calculate terminator longitudes
         418             :        !------------------------------------------------------------------
         419    23376600 :        if( abs(latitude) >= (pid2 - abs(declination)) ) then
         420     1575552 :           if( sign(1._r8,declination) == sign(1._r8,latitude) ) then
         421             :              polar_day = .true.
         422             :              sunoff = 2._r8*twopi
         423             :              sunon  = -twopi
         424             :           else
         425      787776 :              polar_night = .true.
         426             :           end if
         427             :        else
         428    21801048 :           sunoff = acos( -tan(declination)*tan(latitude) )
         429    21801048 :           sunon  = twopi - sunoff
         430             :        end if
         431             :
         432             :        !--------------------------------------------------------
         433             :        !        ... adjust alpha-pinene for diurnal variation
         434             :        !--------------------------------------------------------
         435    23376600 :        if( c10h16_ndx > 0 ) then
         436           0 :           if( has_emis(c10h16_ndx) ) then
         437           0 :              if( .not. polar_night .and. .not. polar_day ) then
         438           0 :                 dayfrac = sunoff / pi
         439           0 :                 sflx(i,c10h16_ndx) = sflx(i,c10h16_ndx) / (.7_r8 + .3_r8*dayfrac)
         440           0 :                 if( loc_angle >= sunoff .and. loc_angle <= sunon ) then
         441           0 :                    sflx(i,c10h16_ndx) = sflx(i,c10h16_ndx) * .7_r8
         442             :                 endif
         443             :              end if
         444             :           end if
         445             :        end if
         446             :
         447             :        !--------------------------------------------------------
         448             :        !        ... adjust isoprene for diurnal variation
         449             :        !--------------------------------------------------------
         450    24865776 :        if( isop_ndx > 0 ) then
         451           0 :           if( has_emis(isop_ndx) ) then
         452           0 :              if( .not. polar_night ) then
         453           0 :                 if( polar_day ) then
         454             :                    iso_off = .8_r8 * pi
         455             :                    iso_on  = 1.2_r8 * pi
         456             :                 else
         457           0 :                    iso_off = .8_r8 * sunoff
         458           0 :                    iso_on  = 2._r8 * pi - iso_off
         459             :                 end if
         460           0 :                 if( loc_angle >= iso_off .and. loc_angle <= iso_on ) then
         461           0 :                    sflx(i,isop_ndx) = 0._r8
         462             :                 else
         463           0 :                    factor = loc_angle - iso_on
         464           0 :                    if( factor <= 0._r8 ) then
         465           0 :                       factor = factor + 2._r8*pi
         466             :                    end if
         467           0 :                    factor = factor / (2._r8*iso_off + 1.e-6_r8)
         468           0 :                    sflx(i,isop_ndx) = sflx(i,isop_ndx) * 2._r8 / iso_off * pi * (sin(pi*factor))**2
         469             :                 end if
         470             :              else
         471           0 :                 sflx(i,isop_ndx) = 0._r8
         472             :              end if
         473             :           end if
         474             :        end if
         475             :
         476             :     end do
         477             :
         478     1489176 :   end subroutine set_srf_emissions
         479             :
         480     1489176 : end module mo_srf_emissions

Generated by: LCOV version 1.14