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

              Line data    Source code
           1             : module mo_extfrc
           2             :   !---------------------------------------------------------------
           3             :   !     ... insitu forcing module
           4             :   !---------------------------------------------------------------
           5             :
           6             :   use shr_kind_mod,  only : r8 => shr_kind_r8
           7             :   use ppgrid,        only : pver, pverp
           8             :   use chem_mods,     only : extcnt, extfrc_lst, frc_from_dataset, adv_mass
           9             :   use spmd_utils,    only : masterproc
          10             :   use cam_abortutils,only : endrun
          11             :   use cam_history,   only : addfld, outfld, add_default, horiz_only
          12             :   use cam_history_support,only : max_fieldname_len
          13             :   use cam_logfile,   only : iulog
          14             :   use tracer_data,   only : trfld,trfile
          15             :   use mo_constants,  only : avogadro
          16             :   use ioFileMod,     only : getfil
          17             :
          18             :   implicit none
          19             :
          20             :   type :: forcing
          21             :      integer           :: frc_ndx
          22             :      real(r8)          :: scalefactor
          23             :      character(len=265):: filename
          24             :      character(len=16) :: species
          25             :      integer                   :: nsectors
          26             :      character(len=32),pointer :: sectors(:)
          27             :      type(trfld), pointer      :: fields(:)
          28             :      type(trfile)              :: file
          29             :   end type forcing
          30             :
          31             :   private
          32             :   public  :: extfrc_inti
          33             :   public  :: extfrc_set
          34             :   public  :: extfrc_timestep_init
          35             :
          36             :   save
          37             :
          38             :   integer, parameter :: time_span = 1
          39             :
          40             :   character(len=256) ::   filename
          41             :
          42             :   type(forcing), allocatable  :: forcings(:)
          43             :   integer :: n_frc_files = 0
          44             :
          45             : contains
          46             :
          47        1536 :   subroutine extfrc_inti( extfrc_specifier, extfrc_type_in, extfrc_cycle_yr, extfrc_fixed_ymd, extfrc_fixed_tod)
          48             :
          49             :     !-----------------------------------------------------------------------
          50             :     !   ... initialize the surface forcings
          51             :     !-----------------------------------------------------------------------
          52             :     use cam_pio_utils, only : cam_pio_openfile, cam_pio_closefile
          53             :     use pio,           only : pio_inquire, pio_inq_varndims, pio_inq_dimid
          54             :     use pio,           only : pio_inq_varname, pio_nowrite, file_desc_t
          55             :     use pio,           only : pio_get_att, PIO_NOERR, PIO_GLOBAL
          56             :     use pio,           only : pio_seterrorhandling, PIO_BCAST_ERROR,PIO_INTERNAL_ERROR
          57             :     use mo_chem_utls,  only : get_extfrc_ndx
          58             :     use chem_mods,     only : frc_from_dataset
          59             :     use tracer_data,   only : trcdata_init
          60             :     use phys_control,  only : phys_getopts
          61             :     use string_utils,  only : GLC
          62             :     use m_MergeSorts,  only : IndexSort
          63             :
          64             :     implicit none
          65             :
          66             :     !-----------------------------------------------------------------------
          67             :     !   ... dummy arguments
          68             :     !-----------------------------------------------------------------------
          69             :     character(len=*), dimension(:), intent(in) :: extfrc_specifier
          70             :     character(len=*), intent(in) :: extfrc_type_in
          71             :     integer  , intent(in)        :: extfrc_cycle_yr
          72             :     integer  , intent(in)        :: extfrc_fixed_ymd
          73             :     integer  , intent(in)        :: extfrc_fixed_tod
          74             :
          75             :     !-----------------------------------------------------------------------
          76             :     !   ... local variables
          77             :     !-----------------------------------------------------------------------
          78             :     integer :: astat
          79             :     integer :: j, l, m, n, i,mm                          ! Indices
          80             :     character(len=16)  :: spc_name
          81        3072 :     character(len=256) :: frc_fnames(size(extfrc_specifier))
          82        3072 :     real(r8)           :: frc_scalefactor(size(extfrc_specifier))
          83        3072 :     character(len=16)  :: frc_species(size(extfrc_specifier))
          84        3072 :     integer            :: frc_indexes(size(extfrc_specifier))
          85        3072 :     integer            :: indx(size(extfrc_specifier))
          86             :
          87             :     integer ::  vid, ndims, nvars, isec, ierr, num_dims_xfrc, dimid
          88        1536 :     logical, allocatable :: is_sector(:)
          89             :     type(file_desc_t) :: ncid
          90             :     character(len=32)  :: varname
          91             :     logical :: unstructured
          92             :
          93             :     character(len=1), parameter :: filelist = ''
          94             :     character(len=1), parameter :: datapath = ''
          95             :     logical         , parameter :: rmv_file = .false.
          96             :     logical  :: history_aerosol
          97             :     logical  :: history_chemistry
          98             :     logical  :: history_cesm_forcing
          99             :
         100             :     character(len=32) :: extfrc_type = ' '
         101             :     character(len=80) :: file_interp_type = ' '
         102             :     character(len=256) :: tmp_string = ' '
         103             :     character(len=32) :: xchr = ' '
         104             :     real(r8) :: xdbl
         105             :     character(len=256) :: locfn
         106             :
         107             :     !-----------------------------------------------------------------------
         108             :
         109             :     call phys_getopts( &
         110             :          history_aerosol_out = history_aerosol, &
         111             :          history_chemistry_out = history_chemistry, &
         112        1536 :          history_cesm_forcing_out = history_cesm_forcing )
         113             :
         114             :     !-----------------------------------------------------------------------
         115             :     !   ... species has insitu forcing ?
         116             :     !-----------------------------------------------------------------------
         117             :
         118             :     !write(iulog,*) 'Species with insitu forcings'
         119        1536 :     mm = 0
         120      155136 :     indx(:) = 0
         121             :
         122       16896 :     count_emis: do n=1,size(extfrc_specifier)
         123             :
         124       16896 :        if ( len_trim(extfrc_specifier(n) ) == 0 ) then
         125             :           exit count_emis
         126             :        endif
         127             :
         128       15360 :        i = scan(extfrc_specifier(n),'->')
         129       15360 :        spc_name = trim(adjustl(extfrc_specifier(n)(:i-1)))
         130             :
         131             :        ! need to parse out scalefactor ...
         132       15360 :        tmp_string = adjustl(extfrc_specifier(n)(i+2:))
         133       15360 :        j = scan( tmp_string, '*' )
         134       15360 :        if (j>0) then
         135           0 :           xchr = tmp_string(1:j-1) ! get the multipler (left of the '*')
         136           0 :           read( xchr, * ) xdbl   ! convert the string to a real
         137           0 :           tmp_string = adjustl(tmp_string(j+1:)) ! get the filepath name (right of the '*')
         138             :        else
         139       15360 :           xdbl = 1._r8
         140             :        endif
         141       15360 :        filename = trim(tmp_string)
         142             :
         143       15360 :        m = get_extfrc_ndx( spc_name )
         144             :
         145       15360 :        if ( m < 1 ) then
         146           0 :           call endrun('extfrc_inti: '//trim(spc_name)// ' does not have an external source')
         147             :        endif
         148             :
         149       15360 :        if ( .not. frc_from_dataset(m) ) then
         150           0 :           call endrun('extfrc_inti: '//trim(spc_name)//' cannot have external forcing from additional dataset')
         151             :        endif
         152             :
         153       15360 :        mm = mm+1
         154       15360 :        frc_species(mm) = spc_name
         155       15360 :        frc_fnames(mm) = filename
         156       15360 :        frc_indexes(mm) = m
         157       15360 :        frc_scalefactor(mm) = xdbl
         158             :
         159       16896 :        indx(n)=n
         160             :
         161             :     enddo count_emis
         162             :
         163        1536 :     n_frc_files = mm
         164             :
         165        1536 :     if( n_frc_files < 1 ) then
         166           0 :        if (masterproc) write(iulog,*) 'There are no species with insitu forcings'
         167           0 :        return
         168             :     end if
         169             :
         170        1536 :     if (masterproc) write(iulog,*) ' '
         171             :
         172             :     !-----------------------------------------------------------------------
         173             :     !   ... allocate forcings type array
         174             :     !-----------------------------------------------------------------------
         175       26112 :     allocate( forcings(n_frc_files), stat=astat )
         176        1536 :     if( astat/= 0 ) then
         177           0 :        write(iulog,*) 'extfrc_inti: failed to allocate forcings array; error = ',astat
         178           0 :        call endrun('extfrc_inti: failed to allocate forcings array')
         179             :     end if
         180             :
         181             :     !-----------------------------------------------------------------------
         182             :     ! Sort the input files so that the emissions sources are summed in the
         183             :     ! same order regardless of the order of the input files in the namelist
         184             :     !-----------------------------------------------------------------------
         185        1536 :     if (n_frc_files > 0) then
         186        1536 :        call IndexSort(n_frc_files, indx, frc_fnames)
         187             :     end if
         188             :
         189             :     !-----------------------------------------------------------------------
         190             :     !   ... setup the forcing type array
         191             :     !-----------------------------------------------------------------------
         192       16896 :     do m=1,n_frc_files
         193       15360 :        forcings(m)%frc_ndx     = frc_indexes(indx(m))
         194       15360 :        forcings(m)%species     = frc_species(indx(m))
         195       15360 :        forcings(m)%filename    = frc_fnames(indx(m))
         196       16896 :        forcings(m)%scalefactor = frc_scalefactor(indx(m))
         197             :     enddo
         198             :
         199       18432 :     do n= 1,extcnt
         200       18432 :        if (frc_from_dataset(n)) then
         201       16896 :           spc_name = extfrc_lst(n)
         202             :           call addfld( trim(spc_name)//'_XFRC', (/ 'lev' /), 'A',  'molec/cm3/s', &
         203       33792 :                'external forcing for '//trim(spc_name) )
         204             :           call addfld( trim(spc_name)//'_CLXF', horiz_only,  'A',  'molec/cm2/s', &
         205       16896 :                'vertically intergrated external forcing for '//trim(spc_name) )
         206             :           call addfld( trim(spc_name)//'_CMXF', horiz_only,  'A',  'kg/m2/s', &
         207       16896 :                'vertically intergrated external forcing for '//trim(spc_name) )
         208       16896 :           if ( history_aerosol .or. history_chemistry ) then
         209       16896 :              call add_default( trim(spc_name)//'_CLXF', 1, ' ' )
         210       16896 :              call add_default( trim(spc_name)//'_CMXF', 1, ' ' )
         211             :           endif
         212       16896 :           if ( history_cesm_forcing .and. spc_name == 'NO2' ) then
         213           0 :              call add_default( trim(spc_name)//'_CLXF', 1, ' ' )
         214           0 :              call add_default( trim(spc_name)//'_CMXF', 1, ' ' )
         215             :           endif
         216             :        endif
         217             :     enddo
         218             :
         219        1536 :     if (masterproc) then
         220             :        !-----------------------------------------------------------------------
         221             :        !        ... diagnostics
         222             :        !-----------------------------------------------------------------------
         223           2 :        write(iulog,*) ' '
         224           2 :        write(iulog,*) 'extfrc_inti: diagnostics'
         225           2 :        write(iulog,*) ' '
         226           2 :        write(iulog,*) 'extfrc timing specs'
         227           2 :        write(iulog,*) 'type = ',extfrc_type
         228           2 :        if( extfrc_type == 'FIXED' ) then
         229           0 :           write(iulog,*) ' fixed date = ', extfrc_fixed_ymd
         230           0 :           write(iulog,*) ' fixed time = ', extfrc_fixed_tod
         231           2 :        else if( extfrc_type == 'CYCLICAL' ) then
         232           0 :           write(iulog,*) ' cycle year = ',extfrc_cycle_yr
         233             :        end if
         234           2 :        write(iulog,*) ' '
         235           2 :        write(iulog,*) 'there are ',n_frc_files,' species with external forcing files'
         236          22 :        do m = 1,n_frc_files
         237          20 :           write(iulog,*) ' '
         238          20 :           write(iulog,*) 'forcing type ',m
         239          20 :           write(iulog,*) 'species = ',trim(forcings(m)%species)
         240          20 :           write(iulog,*) 'frc ndx = ',forcings(m)%frc_ndx
         241          22 :           write(iulog,*) 'filename= ',trim(forcings(m)%filename)
         242             :        end do
         243           2 :        write(iulog,*) ' '
         244             :     endif
         245             :
         246             :     !-----------------------------------------------------------------------
         247             :     ! read emis files to determine number of sectors
         248             :     !-----------------------------------------------------------------------
         249       16896 :     frcing_loop: do m = 1, n_frc_files
         250             :
         251       15360 :        forcings(m)%nsectors = 0
         252             :
         253       15360 :        if (masterproc) then
         254          20 :           write(iulog,'(a,i3,a)') 'extfrc_inti m: ',m,' init file : '//trim(forcings(m)%filename)
         255             :        endif
         256             :
         257       15360 :        call getfil (forcings(m)%filename, locfn, 0)
         258       15360 :        call cam_pio_openfile ( ncid, trim(locfn), PIO_NOWRITE)
         259       15360 :        ierr = pio_inquire (ncid, nVariables=nvars)
         260             :
         261       15360 :        call pio_seterrorhandling(ncid, PIO_BCAST_ERROR)
         262       15360 :        ierr = pio_inq_dimid( ncid, 'ncol', dimid )
         263       15360 :        unstructured = ierr==PIO_NOERR
         264       15360 :        call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)
         265             :
         266       46080 :        allocate(is_sector(nvars))
         267      188928 :        is_sector(:) = .false.
         268             :
         269      188928 :        do vid = 1,nvars
         270             :
         271      173568 :           ierr = pio_inq_varndims (ncid, vid, ndims)
         272      173568 :           if (unstructured) then
         273             :              num_dims_xfrc = 3
         274             :           else
         275           0 :              num_dims_xfrc = 4
         276             :           endif
         277             :
         278      173568 :           if( ndims < num_dims_xfrc ) then
         279             :              cycle
         280       15360 :           elseif( ndims > num_dims_xfrc ) then
         281           0 :              ierr = pio_inq_varname (ncid, vid, varname)
         282           0 :              write(iulog,*) 'extfrc_inti: Skipping variable ', trim(varname),', ndims = ',ndims, &
         283           0 :                   ' , species=',trim(forcings(m)%species)
         284           0 :              cycle
         285             :           end if
         286             :
         287       15360 :           forcings(m)%nsectors = forcings(m)%nsectors+1
         288      188928 :           is_sector(vid)=.true.
         289             :
         290             :        enddo
         291             :
         292       46080 :        allocate( forcings(m)%sectors(forcings(m)%nsectors), stat=astat )
         293       15360 :        if( astat/= 0 ) then
         294           0 :          write(iulog,*) 'extfrc_inti: failed to allocate forcings(m)%sectors array; error = ',astat
         295           0 :          call endrun
         296             :        end if
         297             :
         298       15360 :        isec = 1
         299      188928 :        do vid = 1,nvars
         300      188928 :           if( is_sector(vid) ) then
         301       15360 :              ierr = pio_inq_varname(ncid, vid, forcings(m)%sectors(isec))
         302       15360 :              isec = isec+1
         303             :           endif
         304             :        enddo
         305       15360 :        deallocate(is_sector)
         306             :
         307             :        ! Global attribute 'input_method' overrides the ext_frc_type namelist setting on
         308             :        ! a file-by-file basis.  If the ext_frc file does not contain the 'input_method'
         309             :        ! attribute then the ext_frc_type namelist setting is used.
         310       15360 :        call pio_seterrorhandling(ncid, PIO_BCAST_ERROR)
         311       15360 :        ierr = pio_get_att(ncid, PIO_GLOBAL, 'input_method', file_interp_type)
         312       15360 :        call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)
         313       15360 :        if ( ierr == PIO_NOERR) then
         314           0 :           l = GLC(file_interp_type)
         315           0 :           extfrc_type(1:l) = file_interp_type(1:l)
         316           0 :           extfrc_type(l+1:) = ' '
         317             :        else
         318       15360 :           extfrc_type = trim(extfrc_type_in)
         319             :        endif
         320             :
         321       15360 :        call cam_pio_closefile (ncid)
         322             :
         323       46080 :        allocate(forcings(m)%file%in_pbuf(size(forcings(m)%sectors)))
         324       30720 :        forcings(m)%file%in_pbuf(:) = .false.
         325           0 :        call trcdata_init( forcings(m)%sectors, &
         326             :                           forcings(m)%filename, filelist, datapath, &
         327             :                           forcings(m)%fields,  &
         328             :                           forcings(m)%file, &
         329       16896 :                           rmv_file, extfrc_cycle_yr, extfrc_fixed_ymd, extfrc_fixed_tod, trim(extfrc_type) )
         330             :
         331             :     enddo frcing_loop
         332             :
         333             :
         334        3072 :   end subroutine extfrc_inti
         335             :
         336      741888 :   subroutine extfrc_timestep_init( pbuf2d, state )
         337             :     !-----------------------------------------------------------------------
         338             :     !       ... check serial case for time span
         339             :     !-----------------------------------------------------------------------
         340             :
         341        1536 :     use physics_types,only : physics_state
         342             :     use ppgrid,       only : begchunk, endchunk
         343             :     use tracer_data,  only : advance_trcdata
         344             :     use physics_buffer, only : physics_buffer_desc
         345             :
         346             :     implicit none
         347             :
         348             :     type(physics_state), intent(in):: state(begchunk:endchunk)
         349             :     type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         350             :
         351             :     !-----------------------------------------------------------------------
         352             :     !       ... local variables
         353             :     !-----------------------------------------------------------------------
         354             :     integer :: m
         355             :
         356     4080384 :     do m = 1,n_frc_files
         357     4080384 :        call advance_trcdata( forcings(m)%fields, forcings(m)%file, state, pbuf2d  )
         358             :     end do
         359             :
         360      370944 :   end subroutine extfrc_timestep_init
         361             :
         362     1489176 :   subroutine extfrc_set( lchnk, zint, frcing, ncol )
         363             :
         364             :     !--------------------------------------------------------
         365             :     !   ... form the external forcing
         366             :     !--------------------------------------------------------
         367      370944 :     use mo_chem_utls,  only : get_spc_ndx
         368             :
         369             :     implicit none
         370             :
         371             :     !--------------------------------------------------------
         372             :     !   ... dummy arguments
         373             :     !--------------------------------------------------------
         374             :     integer,  intent(in)    :: ncol                  ! columns in chunk
         375             :     integer,  intent(in)    :: lchnk                 ! chunk index
         376             :     real(r8), intent(in)    :: zint(ncol, pverp)                  ! interface geopot above surface (km)
         377             :     real(r8), intent(inout) :: frcing(ncol,pver,extcnt)   ! insitu forcings (molec/cm^3/s)
         378             :
         379             :     !--------------------------------------------------------
         380             :     !   ... local variables
         381             :     !--------------------------------------------------------
         382             :     integer  ::  m, n
         383             :     character(len=max_fieldname_len) :: xfcname
         384     1489176 :     real(r8) :: frcing_col(1:ncol), frcing_col_kg(1:ncol)
         385             :     integer  :: k, isec
         386             :     real(r8),parameter :: km_to_cm = 1.e5_r8
         387             :     real(r8),parameter :: cm2_to_m2 = 1.e4_r8
         388             :     real(r8),parameter :: kg_to_g = 1.e-3_r8
         389             :     real(r8) :: molec_to_kg
         390             :     integer  :: spc_ndx
         391             :
         392     1489176 :     if( n_frc_files < 1 .or. extcnt < 1 ) then
         393             :        return
         394             :     end if
         395             :
         396 25455558960 :     frcing(:,:,:) = 0._r8
         397             :
         398             :     !--------------------------------------------------------
         399             :     !   ... set non-zero forcings
         400             :     !--------------------------------------------------------
         401    16380936 :     file_loop : do m = 1,n_frc_files
         402             :
         403    14891760 :        n = forcings(m)%frc_ndx
         404             :
         405    31272696 :        do isec = 1,forcings(m)%nsectors
         406 23154955200 :           frcing(:ncol,:,n) = frcing(:ncol,:,n) + forcings(m)%scalefactor*forcings(m)%fields(isec)%data(:ncol,:,lchnk)
         407             :        enddo
         408             :
         409             :     enddo file_loop
         410             :
         411    17870112 :     frc_loop : do n = 1,extcnt
         412    17870112 :        if (frc_from_dataset(n)) then
         413             :
         414    16380936 :           xfcname = trim(extfrc_lst(n))//'_XFRC'
         415    16380936 :           call outfld( xfcname, frcing(:ncol,:,n), ncol, lchnk )
         416             :
         417    16380936 :           spc_ndx = get_spc_ndx( extfrc_lst(n) )
         418    16380936 :           molec_to_kg = adv_mass( spc_ndx ) / avogadro *cm2_to_m2 * kg_to_g
         419             :
         420   273523536 :           frcing_col(:ncol) = 0._r8
         421   273523536 :           frcing_col_kg(:ncol) = 0._r8
         422  1539807984 :           do k = 1,pver
         423 25437688848 :              frcing_col(:ncol) = frcing_col(:ncol) + frcing(:ncol,k,n)*(zint(:ncol,k)-zint(:ncol,k+1))*km_to_cm
         424 25454069784 :              frcing_col_kg(:ncol) = frcing_col_kg(:ncol) + frcing(:ncol,k,n)*(zint(:ncol,k)-zint(:ncol,k+1))*km_to_cm*molec_to_kg
         425             :           enddo
         426             :
         427    16380936 :           xfcname = trim(extfrc_lst(n))//'_CLXF'
         428    16380936 :           call outfld( xfcname, frcing_col(:ncol), ncol, lchnk )
         429    16380936 :           xfcname = trim(extfrc_lst(n))//'_CMXF'
         430    16380936 :           call outfld( xfcname, frcing_col_kg(:ncol), ncol, lchnk )
         431             :        endif
         432             :     end do frc_loop
         433             :
         434             :   end subroutine extfrc_set
         435             :
         436             :
         437           0 : end module mo_extfrc

Generated by: LCOV version 1.14