ctsm5.3.033: MIMICS updates

python/ctsm/site_and_regional/single_point_case.py

@@ -597,7 +597,7 @@ def extract_datm_at(self, file_in, file_out):
         f_out = f_out.expand_dims(["lat", "lon"])
 
         # specify dimension order
-        f_out = f_out.transpose("scalar", "time", "lat", "lon")
+        f_out = f_out.transpose("time", "lat", "lon")
 
         # update attributes
         self.update_metadata(f_out)
@@ -648,46 +648,42 @@ def create_datm_at_point(self, datm_tuple: DatmFiles, datm_syr, datm_eyr, datm_s
         tpqwfiles = []
         for year in range(datm_syr, datm_eyr + 1):
             ystr = str(year)
-            for month in range(FIRST_MONTH, LAST_MONTH + 1):
-                mstr = str(month)
-                if month < 10:
-                    mstr = "0" + mstr
 
-                dtag = ystr + "-" + mstr
+            dtag = ystr 
 
-                fsolar = os.path.join(
-                    datm_tuple.indir,
-                    datm_tuple.dir_solar,
-                    "{}{}.nc".format(datm_tuple.tag_solar, dtag),
-                )
-                fsolar2 = "{}{}.{}.nc".format(datm_tuple.tag_solar, self.tag, dtag)
-                fprecip = os.path.join(
-                    datm_tuple.indir,
-                    datm_tuple.dir_prec,
-                    "{}{}.nc".format(datm_tuple.tag_prec, dtag),
-                )
-                fprecip2 = "{}{}.{}.nc".format(datm_tuple.tag_prec, self.tag, dtag)
-                ftpqw = os.path.join(
-                    datm_tuple.indir,
-                    datm_tuple.dir_tpqw,
-                    "{}{}.nc".format(datm_tuple.tag_tpqw, dtag),
-                )
-                ftpqw2 = "{}{}.{}.nc".format(datm_tuple.tag_tpqw, self.tag, dtag)
-
-                outdir = os.path.join(self.out_dir, datm_tuple.outdir)
-                infile += [fsolar, fprecip, ftpqw]
-                outfile += [
-                    os.path.join(outdir, fsolar2),
-                    os.path.join(outdir, fprecip2),
-                    os.path.join(outdir, ftpqw2),
-                ]
-                solarfiles.append(
-                    os.path.join("${}".format(USRDAT_DIR), datm_tuple.outdir, fsolar2)
-                )
-                precfiles.append(
-                    os.path.join("${}".format(USRDAT_DIR), datm_tuple.outdir, fprecip2)
-                )
-                tpqwfiles.append(os.path.join("${}".format(USRDAT_DIR), datm_tuple.outdir, ftpqw2))
+            fsolar = os.path.join(
+                datm_tuple.indir,
+                datm_tuple.dir_solar,
+                "{}{}.nc".format(datm_tuple.tag_solar, dtag),
+            )
+            fsolar2 = "{}{}.{}.nc".format(datm_tuple.tag_solar, self.tag, dtag)
+            fprecip = os.path.join(
+                datm_tuple.indir,
+                datm_tuple.dir_prec,
+                "{}{}.nc".format(datm_tuple.tag_prec, dtag),
+            )
+            fprecip2 = "{}{}.{}.nc".format(datm_tuple.tag_prec, self.tag, dtag)
+            ftpqw = os.path.join(
+                datm_tuple.indir,
+                datm_tuple.dir_tpqw,
+                "{}{}.nc".format(datm_tuple.tag_tpqw, dtag),
+            )
+            ftpqw2 = "{}{}.{}.nc".format(datm_tuple.tag_tpqw, self.tag, dtag)
+
+            outdir = os.path.join(self.out_dir, datm_tuple.outdir)
+            infile += [fsolar, fprecip, ftpqw]
+            outfile += [
+                os.path.join(outdir, fsolar2),
+                os.path.join(outdir, fprecip2),
+                os.path.join(outdir, ftpqw2),
+            ]
+            solarfiles.append(
+                os.path.join("${}".format(USRDAT_DIR), datm_tuple.outdir, fsolar2)
+            )
+            precfiles.append(
+                os.path.join("${}".format(USRDAT_DIR), datm_tuple.outdir, fprecip2)
+            )
+            tpqwfiles.append(os.path.join("${}".format(USRDAT_DIR), datm_tuple.outdir, ftpqw2))
 
         for idx, out_f in enumerate(outfile):
             logger.debug(out_f)

src/biogeochem/CNCStateUpdate1Mod.F90

@@ -9,6 +9,7 @@ module CNCStateUpdate1Mod
   use clm_varpar                         , only : ndecomp_cascade_transitions, nlevdecomp
   use clm_time_manager                   , only : get_step_size_real
   use clm_varpar                         , only : i_litr_min, i_litr_max, i_cwd
+  use clm_varpar                         , only : i_met_lit, i_str_lit, i_phys_som, i_chem_som
   use pftconMod                          , only : npcropmin, nc3crop, pftcon
   use abortutils                         , only : endrun
   use decompMod                          , only : bounds_type
@@ -17,7 +18,7 @@ module CNCStateUpdate1Mod
   use CropType                           , only : crop_type
   use CropReprPoolsMod                   , only : nrepr, repr_grain_min, repr_grain_max
   use CropReprPoolsMod                   , only : repr_structure_min, repr_structure_max
-  use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con, use_soil_matrixcn
+  use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con, decomp_method, mimics_decomp, use_soil_matrixcn
   use SoilBiogeochemCarbonFluxType       , only : soilbiogeochem_carbonflux_type
   use SoilBiogeochemCarbonStateType      , only : soilbiogeochem_carbonstate_type
   use PatchType                          , only : patch
@@ -170,6 +171,7 @@ subroutine CStateUpdate1( num_soilc, filter_soilc, num_soilp, filter_soilp, &
     associate(                                                               & 
          ivt                   => patch%itype                                , & ! Input:  [integer  (:)     ]  patch vegetation type                                
 
+         mimics_fi             => pftcon%mimics_fi                         , & ! Input: MIMICS parameter fi
          woody                 => pftcon%woody                             , & ! Input:  binary flag for woody lifeform (1=woody, 0=not woody)
 
          cascade_donor_pool    => decomp_cascade_con%cascade_donor_pool    , & ! Input:  [integer  (:)     ]  which pool is C taken from for a given decomposition step
@@ -206,10 +208,21 @@ subroutine CStateUpdate1( num_soilc, filter_soilc, num_soilp, filter_soilp, &
                !
                if (.not. use_soil_matrixcn) then
                   ! phenology and dynamic land cover fluxes
-                  do i = i_litr_min, i_litr_max
-                     cf_soil%decomp_cpools_sourcesink_col(c,j,i) = &
-                          cf_veg%phenology_c_to_litr_c_col(c,j,i) * dt
-                  end do
+                  if (decomp_method == mimics_decomp) then
+                     do i = i_litr_min, i_litr_max  ! in MIMICS these are 1 and 2
+                        cf_soil%decomp_cpools_sourcesink_col(c,j,i) = (1 - mimics_fi(i)) * &
+                           cf_veg%phenology_c_to_litr_c_col(c,j,i) * dt
+                     end do
+                     cf_soil%decomp_cpools_sourcesink_col(c,j,i_phys_som) = mimics_fi(1) * &
+                        cf_veg%phenology_c_to_litr_c_col(c,j,i_met_lit) * dt
+                     cf_soil%decomp_cpools_sourcesink_col(c,j,i_chem_som) = mimics_fi(2) * &
+                        cf_veg%phenology_c_to_litr_c_col(c,j,i_str_lit) * dt
+                  else
+                     do i = i_litr_min, i_litr_max
+                        cf_soil%decomp_cpools_sourcesink_col(c,j,i) = &
+                             cf_veg%phenology_c_to_litr_c_col(c,j,i) * dt
+                     end do
+                  end if
 
                   ! NOTE(wjs, 2017-01-02) This used to be set to a non-zero value, but the
                   ! terms have been moved to CStateUpdateDynPatch. I think this is zeroed every

src/biogeochem/CNNStateUpdate1Mod.F90

@@ -12,8 +12,9 @@ module CNNStateUpdate1Mod
   use clm_time_manager                , only : get_step_size_real
   use clm_varpar                      , only : nlevdecomp
   use clm_varpar                      , only : i_litr_min, i_litr_max, i_cwd
+  use clm_varpar                      , only : i_met_lit, i_str_lit, i_phys_som, i_chem_som
   use clm_varctl                      , only : iulog, use_nitrif_denitrif
-  use SoilBiogeochemDecompCascadeConType, only : use_soil_matrixcn
+  use SoilBiogeochemDecompCascadeConType, only : decomp_method, mimics_decomp, use_soil_matrixcn
   use CNSharedParamsMod               , only : use_matrixcn
   use clm_varcon                      , only : nitrif_n2o_loss_frac
   use pftconMod                       , only : npcropmin, pftcon
@@ -129,6 +130,7 @@ subroutine NStateUpdate1(num_soilc, filter_soilc, num_soilp, filter_soilp, &
     associate(                                                                   & 
          ivt                   => patch%itype                                    , & ! Input:  [integer  (:)     ]  patch vegetation type                                
 
+         mimics_fi             => pftcon%mimics_fi                             , & ! Input:  MIMICS parameter fi
          woody                 => pftcon%woody                                 , & ! Input:  binary flag for woody lifeform (1=woody, 0=not woody)
 
          nf_veg                => cnveg_nitrogenflux_inst                      , & ! Input:
@@ -163,10 +165,21 @@ subroutine NStateUpdate1(num_soilc, filter_soilc, num_soilp, filter_soilp, &
             ! State update without the matrix solution
             !
             if (.not. use_soil_matrixcn) then ! to be consistent with C
-               do i = i_litr_min, i_litr_max
-                  nf_soil%decomp_npools_sourcesink_col(c,j,i) = &
-                     nf_veg%phenology_n_to_litr_n_col(c,j,i) * dt
-               end do
+               if (decomp_method == mimics_decomp) then
+                  do i = i_litr_min, i_litr_max  ! in MIMICS these are 1 and 2
+                     nf_soil%decomp_npools_sourcesink_col(c,j,i) = (1 - mimics_fi(i)) * &
+                        nf_veg%phenology_n_to_litr_n_col(c,j,i) * dt
+                  end do
+                  nf_soil%decomp_npools_sourcesink_col(c,j,i_phys_som) = mimics_fi(1) * &
+                     nf_veg%phenology_n_to_litr_n_col(c,j,i_met_lit) * dt
+                  nf_soil%decomp_npools_sourcesink_col(c,j,i_chem_som) = mimics_fi(2) * &
+                     nf_veg%phenology_n_to_litr_n_col(c,j,i_str_lit) * dt
+               else
+                  do i = i_litr_min, i_litr_max
+                     nf_soil%decomp_npools_sourcesink_col(c,j,i) = &
+                        nf_veg%phenology_n_to_litr_n_col(c,j,i) * dt
+                  end do
+               end if
 
                ! NOTE(wjs, 2017-01-02) This used to be set to a non-zero value, but the
                ! terms have been moved to CStateUpdateDynPatch. I think this is zeroed every

src/main/clm_varpar.F90

@@ -82,6 +82,9 @@ module clm_varpar
   integer, public            :: i_litr_min    = -9    ! min index of litter pools; overwritten in SoilBiogeochemDecompCascade*Mod
   integer, public            :: i_litr_max    = -9    ! max index of litter pools; overwritten in SoilBiogeochemDecompCascade*Mod
   integer, public            :: i_met_lit     = -9    ! index of metabolic litter pool; overwritten in SoilBiogeochemDecompCascade*Mod
+  integer, public            :: i_str_lit     = -9    ! index of structural litter pool; overwritten in SoilBiogeochemDecompCascade*Mod
+  integer, public            :: i_phys_som    = -9    ! index of physically protected Soil Organic Matter (SOM); overwritten in SoilBiogeochemDecompCascade*Mod
+  integer, public            :: i_chem_som    = -9    ! index of chemically protected Soil Organic Matter (SOM); overwritten in SoilBiogeochemDecompCascade*Mod
   integer, public            :: i_cop_mic     = -9    ! index of copiotrophic microbial pool; overwritten in SoilBiogeochemDecompCascade*Mod
   integer, public            :: i_oli_mic     = -9    ! index of oligotrophic microbial pool; overwritten in SoilBiogeochemDecompCascade*Mod
   integer, public            :: i_cwd         = -9    ! index of cwd pool; overwritten in SoilBiogeochemDecompCascade*Mod

src/main/pftconMod.F90

@@ -165,6 +165,9 @@ module pftconMod
      real(r8), allocatable :: crit_onset_gdd_sf(:)! scale factor for crit_onset_gdd
      real(r8), allocatable :: ndays_on(:)         ! number of days to complete leaf onset
 
+     ! MIMICS
+     real(r8), allocatable :: mimics_fi(:)
+
      !  crop
 
      ! These arrays give information about the merge of unused crop types to the types CLM
@@ -504,6 +507,7 @@ subroutine InitAllocate (this)
     allocate( this%taper         (0:mxpft) )
     allocate( this%rstem_per_dbh (0:mxpft) )
     allocate( this%wood_density  (0:mxpft) )
+    allocate( this%mimics_fi(2) )
     allocate( this%crit_onset_gdd_sf (0:mxpft) )
     allocate( this%ndays_on      (0:mxpft) )
 
@@ -1095,6 +1099,9 @@ subroutine InitRead(this)
     !
     ! clm 5 nitrogen variables
     !
+    call ncd_io('mimics_fi',this%mimics_fi, 'read', ncid, readvar=readv)
+    if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(sourcefile, __LINE__))
+
     if (use_flexibleCN) then
        call ncd_io('i_vcad', this%i_vcad, 'read', ncid, readvar=readv) 
        if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(sourcefile, __LINE__)) 
@@ -1132,6 +1139,7 @@ subroutine InitRead(this)
        if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(sourcefile, __LINE__))
        call ncd_io('wood_density',this%wood_density, 'read', ncid, readvar=readv)
        if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(sourcefile, __LINE__))
+
     else
        this%dbh = 0.0_r8
        this%fbw = 0.0_r8
@@ -1594,6 +1602,7 @@ subroutine Clean(this)
     deallocate( this%rstem_per_dbh)
     deallocate( this%wood_density)
     deallocate( this%taper)
+    deallocate( this%mimics_fi)
     deallocate( this%crit_onset_gdd_sf)
     deallocate( this%ndays_on)
   end subroutine Clean

src/soilbiogeochem/SoilBiogeochemDecompCascadeMIMICSMod.F90

@@ -10,7 +10,7 @@ module SoilBiogeochemDecompCascadeMIMICSMod
   use shr_const_mod                      , only : SHR_CONST_TKFRZ
   use shr_log_mod                        , only : errMsg => shr_log_errMsg
   use clm_varpar                         , only : nlevdecomp, ndecomp_pools_max
-  use clm_varpar                         , only : i_met_lit, i_cop_mic, i_oli_mic, i_cwd
+  use clm_varpar                         , only : i_phys_som, i_chem_som, i_str_lit, i_met_lit, i_cop_mic, i_oli_mic, i_cwd
   use clm_varpar                         , only : i_litr_min, i_litr_max, i_cwdl2
   use clm_varctl                         , only : iulog, spinup_state, anoxia, use_lch4, use_fates
   use clm_varcon                         , only : zsoi
@@ -48,10 +48,7 @@ module SoilBiogeochemDecompCascadeMIMICSMod
   real(r8), private, allocatable :: fphys_m1(:,:)
   real(r8), private, allocatable :: fphys_m2(:,:)
   real(r8), private, allocatable :: p_scalar(:,:)
-  integer, private :: i_phys_som  ! index of physically protected Soil Organic Matter (SOM)
-  integer, private :: i_chem_som  ! index of chemically protected SOM
   integer, private :: i_avl_som  ! index of available (aka active) SOM
-  integer, private :: i_str_lit  ! index of structural litter pool
   integer, private :: i_l1m1  ! indices of transitions, eg l1m1: litter 1 -> first microbial pool
   integer, private :: i_l1m2
   integer, private :: i_l2m1
@@ -1293,14 +1290,14 @@ subroutine decomp_rates_mimics(bounds, num_bgc_soilc, filter_bgc_soilc, &
             decomp_k(c,j,i_chem_som) = (term_1 + term_2) * w_d_o_scalars
 
             ! Currently, mimics_densdep = 1 so as to have no effect
-            decomp_k(c,j,i_cop_mic) = tau_m1 * &
-                   m1_conc**(mimics_densdep - 1.0_r8) * w_d_o_scalars
+            decomp_k(c,j,i_cop_mic) = tau_m1 * m1_conc**(mimics_densdep) 
+
             favl = min(1.0_r8, max(0.0_r8, 1.0_r8 - fphys_m1(c,j) - fchem_m1))
             pathfrac_decomp_cascade(c,j,i_m1s1) = favl
             pathfrac_decomp_cascade(c,j,i_m1s2) = fchem_m1
 
-            decomp_k(c,j,i_oli_mic) = tau_m2 * &
-                   m2_conc**(mimics_densdep - 1.0_r8) * w_d_o_scalars
+            decomp_k(c,j,i_oli_mic) = tau_m2 * m2_conc**(mimics_densdep)
+
             favl = min(1.0_r8, max(0.0_r8, 1.0_r8 - fphys_m2(c,j) - fchem_m2))
             pathfrac_decomp_cascade(c,j,i_m2s1) = favl
             pathfrac_decomp_cascade(c,j,i_m2s2) = fchem_m2