Computing and sending drydep velocities from CTSM to CAM

[mvertens]

drv_flds_in

• the namelist file drv_flds_in contains the namelist group drydep_inparm which is read by CTSM to determine if drydep velocities will be sent by CTSM to the mediator
• When running a case with CAM and CTSM,
  • CAM determines the drydep_inparm setting in drv_flds_in and both cam and clm use the file drv_flds_in and the namelist group drydep_inparm to determine if drydep velocities will be sent from clm -> cmeps -> cam
• When running a case with DATM and CTSM
  • to turn on having CTSM send drydep velocities to the mediator simply invoke the following command after case creation
    ./xmlchange --append CLM_BLDNML_OPTS='-drydep'

---

shr_drydep_mod

• clm and cam read in the drydep information using the routine
  ~/$SRCROOT/components/cmeps/cesm/nuopc_cap_share/shr_drydep_mod.F90
• shr_drydep_mod.F90 reads in from the file drv_flds_in the namelist group dry_dep_inparm - which for the n1850 case is

&drydep_inparm
  dep_data_file = '/cluster/shared/noresm/inputdata/atm/cam/chem/trop_mozart/dvel/dep_data_c20221208.nc'
  drydep_list = 'DMS', 'H2O2', 'H2SO4', 'SO2'
/ 

• dep_data_c20221208.nc contains the following fields:

netcdf dep_data_c20221208 {
dimensions:
        n_species_table = 193 ;
        NHen = 6 ;
        chars = 16 ;
variables:
        char species_name_table(n_species_table, chars) ;
                species_name_table:long_name = "names of species in dry deposition tables" ;
        double dfoxd(n_species_table) ;
                dfoxd:long_name = "data for foxd (reactivity factor for oxidation)" ;
        double dheff(n_species_table, NHen) ;
                dheff:long_name = "data for effective Henry\'s Law coefficient" ;
        double mol_wghts(n_species_table) ;
                mol_wghts:long_name = "species molecular mass" ;
                mol_wghts:units = "grams/mole" ;

Note that 'DMS', 'H2O2', 'H2SO4', 'SO2' are all entries in species_name_table

• shr_drydep_mod.F90 contains the following hard-code module parameters

  integer, private, parameter :: NSeas = 5                 ! Number of seasons
  integer, public,  parameter :: NLUse = 11                ! Number of land-use types

and the following module arrays are used by CTSM in DryDepVelocity.F90 and depend on the season and land type

  real(r8), public, protected, dimension(NSeas,NLUse) :: ri, rlu, rac, rgss, rgso, rcls, rclo
  !--- ri:   Richardson number                      (dimensionless)
  !--- rlu:  Resistance of leaves in upper canopy   (s.m-1)
  !--- rac:  Aerodynamic resistance to lower canopy (s.m-1)
  !--- rgss: Ground surface resistance for SO2      (s.m-1)
  !--- rgso: Ground surface resistance for O3       (s.m-1)
  !--- rcls: Lower canopy resistance for SO2        (s.m-1)
  !--- rclo: Lower canopy resistance for O3         (s.m-1)

where the indices in the above arrays correspond to the following

  ! index #1 : season
  !           1 -> midsummer with lush vegetation
  !           2 -> autumn with unharvested cropland
  !           3 -> late autumn after frost, no snow
  !           4 -> winter, snow on ground, and subfreezing
  !           5 -> transitional spring with partially green short annuals
  !
  ! index #2 : landuse type (Wesley veg type)
  !           1 -> urban land
  !           2 -> agricultural land
  !           3 -> range land
  !           4 -> deciduous forest
  !           5 -> coniferous forest
  !           6 -> mixed forest including wetland
  !           7 -> water, both salt and fresh
  !           8 -> barren land, mostly desert
  !           9 -> nonforested wetland
  !           10 -> mixed agricultural and range land
  !           11 -> rocky open areas with low growing shrubs

The season index is calculated as follows:

     pft_loop: do pi = bounds%begp,bounds%endp
        l = patch%landunit(pi)
        active: if (patch%active(pi)) then
           .......

           !mlaidiff=jan-feb
           minlai=minval(annlai(:,pi))
           maxlai=maxval(annlai(:,pi))

           index_season = -1

           if ( lun%itype(l) /= istsoil )then
              if ( lun%itype(l) == istice ) then
                 wesveg       = 8
                 index_season = 4
              elseif ( lun%itype(l) == istdlak ) then
                 wesveg       = 7
                 index_season = 4
              elseif ( lun%itype(l) == istwet ) then
                 wesveg       = 9
                 index_season = 2
              elseif ( lun%urbpoi(l) ) then
                 wesveg       = 1
                 index_season = 2
              end if
           else if ( snow_depth(c) > 0 ) then
              index_season = 4
           else if(elai(pi) > 0.5_r8*maxlai) then
              index_season = 1
           endif
           if (index_season<0) then
              if (elai(pi) < (minlai+0.05_r8*(maxlai-minlai))) then
                 index_season = 3
              endif
           endif

           if (index_season<0) then
              if (mlaidiff(pi) > 0.0_r8) then
                 index_season = 2
              elseif (mlaidiff(pi) < 0.0_r8) then
                 index_season = 5
     	       elseif (mlaidiff(pi).eq.0.0_r8) then
                 index_season = 3
              endif
           endif

---

CTSM DryDepVelocity.F90

The drydep velocity is computed as follows:

• a loop over patches is at the top level and for each patch:
  • the patch vegetation type is mapped to the landuse type in shr_drydep_mod (i.e. the Wesley veg type above).
  • the patch seasonal index is determined from mlaidiff (input argument) that is obtained by reading in the monthly lai data
  • the surface resistance is computed which depends on both land type and species
    • as an example, ground resistance is computed as follows

rgsx(ispec) = 1._r8/((heff(ispec)/(1.e5_r8*(rgss(index_season,wesveg)+cts))) + (foxd(ispec)/(rgso(index_season,wesveg)+cts)))

where heff and foxd are read in from dep_data_c20221208.nc

Suggestion from @mvdebolskiy

• NOCOMP mode in FATES: In this mode, all processes are active except for light competition between PFTs. Each PFT is given a total patch area to grow on, but unlike FATES-SP mode, disturbance can occur on each patch and thus the space allocated to each PFT may be split into one or more patches based on disturbance history. Each patch has a PFT label, and only that PFT is allowed to grow on the patch. Cohorts of a given PFT compete against each other for canopy access and thus light.
• Just make a new map for nocomp. and initialize fates pft numver to equal hlm_pfts same way as in sp mode.

[mvertens]

@rosiealice @mvdebolskiy @kjetilaas @maritsandstad - I've added a new section at the top called drv_flds_in that clarifies the recommendation that @mvdebolskiy made at our meeting this morning.