From bc6d32127f8b2aa484b7b543e0d4a5c1e8731a6d Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Mon, 8 Sep 2025 18:21:06 -0400 Subject: [PATCH 01/29] mix of fixes for catchCN CLM 5.1 --- .../CLM51/CNBalanceCheckMod.F90 | 31 ++++- .../CLM51/CNCLM_DriverMod.F90 | 3 + .../CLM51/CNCLM_init_mod.F90 | 2 +- .../CLM51/CNProductsMod.F90 | 24 +++- .../CLM51/CNVegCarbonFluxType.F90 | 20 ++- .../CLM51/CNVegCarbonStateType.F90 | 31 ++++- .../CLM51/CNVegNitrogenStateType.F90 | 69 ++++++---- .../CLM51/CNVegStateType.F90 | 1 + .../CLM51/CanopyStateType.F90 | 3 +- .../CLM51/PhotosynthesisMod.F90 | 16 +-- .../CLM51/SoilBiogeochemNitrogenStateType.F90 | 71 +++++++--- .../CLM51/SurfaceAlbedoType.F90 | 28 ++-- .../Utils/mk_restarts/CatchmentCNRst.F90 | 124 +++++++++++++++--- 13 files changed, 318 insertions(+), 105 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNBalanceCheckMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNBalanceCheckMod.F90 index c10a48c9c..d64e96d02 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNBalanceCheckMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNBalanceCheckMod.F90 @@ -523,9 +523,9 @@ subroutine NBalanceCheck(this, bounds, num_soilc, filter_soilc, & end if if (abs(col_errnb(c)) > this%nwarning) then - ! write(iulog,*) 'nbalance warning at c =', c, col_errnb(c), col_endnb(c) - ! write(iulog,*)'inputs,ffix,nfix,ndep = ',ffix_to_sminn(c)*dt,nfix_to_sminn(c)*dt,ndep_to_sminn(c)*dt - ! write(iulog,*)'outputs,lch,roff,dnit = ',smin_no3_leached(c)*dt, smin_no3_runoff(c)*dt,f_n2o_nit(c)*dt + write(iulog,*) 'nbalance warning at c =', c, col_errnb(c), col_endnb(c) + write(iulog,*)'inputs,ffix,nfix,ndep = ',ffix_to_sminn(c)*dt,nfix_to_sminn(c)*dt,ndep_to_sminn(c)*dt + write(iulog,*)'outputs,lch,roff,dnit = ',smin_no3_leached(c)*dt, smin_no3_runoff(c)*dt,f_n2o_nit(c)*dt end if end do ! end of columns loop @@ -540,8 +540,26 @@ subroutine NBalanceCheck(this, bounds, num_soilc, filter_soilc, & write(iulog,*)'input mass = ',col_ninputs(c)*dt write(iulog,*)'output mass = ',col_noutputs(c)*dt write(iulog,*)'net flux = ',(col_ninputs(c)-col_noutputs(c))*dt - write(iulog,*)'inputs,ffix,nfix,ndep = ',ffix_to_sminn(c)*dt,nfix_to_sminn(c)*dt,ndep_to_sminn(c)*dt - write(iulog,*)'outputs,ffix,nfix,ndep = ',smin_no3_leached(c)*dt, smin_no3_runoff(c)*dt,f_n2o_nit(c)*dt + write(iulog,*)'column nbalance error = ', col_errnb(c), c + write(iulog,*)'Latdeg,Londeg = ', grc%latdeg(col%gridcell(c)), grc%londeg(col%gridcell(c)) + write(iulog,*)'begnb = ', col_begnb(c) + write(iulog,*)'endnb = ', col_endnb(c) + write(iulog,*)'delta store = ', col_endnb(c)-col_begnb(c) + write(iulog,*)'input mass (dt) = ', col_ninputs(c)*dt + write(iulog,*)'output mass (dt) = ', col_noutputs(c)*dt + write(iulog,*)'net flux (dt) = ', (col_ninputs(c)-col_noutputs(c))*dt + + ! Full IN terms (rates) + write(iulog,*) ' IN rates: ndep=', ndep_to_sminn(c), ' nfix=', nfix_to_sminn(c), & + ' ffix=', ffix_to_sminn(c), ' supplement=', supplement_to_sminn(c), & + ' fert=', fert_to_sminn(c), ' soyfix=', soyfixn_to_sminn(c) + + ! Full OUT terms (rates) exactly as used in col_noutputs(c) + write(iulog,*) 'OUT rates: denit=', denit(c), ' fireN=', col_fire_nloss(c), & + ' wood=', wood_harvestn(c), ' grain=', grainn_to_cropprodn(c), & + ' f_n2o=', f_n2o_nit(c), ' NO3_leach=', smin_no3_leached(c), & + ' NO3_runoff=', smin_no3_runoff(c), ' somN_leach(NEG)=', som_n_leached(c) + @@ -597,10 +615,9 @@ subroutine NBalanceCheck(this, bounds, num_soilc, filter_soilc, & end if if (abs(grc_errnb(g)) > this%nwarning) then - !write(iulog,*) 'nbalance warning at g =', g, grc_errnb(g), grc_endnb(g) + write(iulog,*) 'nbalance warning at g =', g, grc_errnb(g), grc_endnb(g) end if end do - if (err_found) then g = err_index write(iulog,*) 'gridcell nbalance error =', grc_errnb(g), g diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 index 3536ef2f5..0d5aef3dd 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 @@ -548,6 +548,9 @@ subroutine CN_exit(nch,ityp,fveg,cncol,cnpft) integer, dimension(8) :: decomp_cpool_cncol_index = (/ 3, 4, 5, 2, 10, 11, 12, 13 /) integer, dimension(8) :: decomp_npool_cncol_index = (/ 18, 19, 20, 17,25, 26, 27, 28 /) !---------------- + ! Clean slate: zero arrays so any unassigned slots are safe in the restart + cncol(:,:,:) = 0.0 + cnpft(:,:,:,:) = 0.0 n = 0 np = 0 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_init_mod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_init_mod.F90 index 0ae6280e0..1bc97bb76 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_init_mod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_init_mod.F90 @@ -201,7 +201,7 @@ subroutine CN_init(NLFilename, nch,ityp,fveg,cncol,cnpft,lats,lons,dtcn,paramfil call solarabs_inst%Init (bounds) - call surfalb_inst%Init (bounds, nch, cncol, cnpft) + call surfalb_inst%Init (bounds, nch, ityp, fveg, cncol, cnpft) call ozone_inst%Init (bounds) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNProductsMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNProductsMod.F90 index 416316c82..96ac1813c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNProductsMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNProductsMod.F90 @@ -74,7 +74,21 @@ module CNProductsMod __FILE__ contains - + pure elemental logical function valid(v) + use shr_kind_mod, only : r8 => shr_kind_r8 + use clm_varcon , only : spval + real(r8), intent(in) :: v + real(r8), parameter :: CF_FILL = 9.96921e36_r8 + real(r8), parameter :: BADHUGE = 1.0e20_r8 ! tighten from 1.0e30 + valid = (v == v) .and. (abs(v) < BADHUGE) .and. & + (abs(v - spval) > 1.0e-6_r8*max(1._r8,abs(spval))) .and. & + (abs(v - CF_FILL) > 1.0e-6_r8*max(1._r8,abs(CF_FILL))) + end function valid + + pure elemental real(r8) function safe(v) + real(r8), intent(in) :: v + safe = merge(v, 0._r8, valid(v)) + end function safe !-------------------------------------------------------------- subroutine Init(this, bounds, nch, cncol, species, rc) @@ -147,11 +161,11 @@ subroutine Init(this, bounds, nch, cncol, species, rc) do nz = 1,num_zon ! CN zone loop if (trim(species) == 'C') then - this%prod100_grc(nc) = this%prod100_grc(nc) + cncol(nc,nz,7)*CN_zone_weight(nz) - this%prod10_grc(nc) = this%prod10_grc(nc) + cncol(nc,nz,8)*CN_zone_weight(nz) + this%prod100_grc(nc) = this%prod100_grc(nc) + safe( real(cncol(nc,nz, 7), r8) ) + this%prod10_grc (nc) = this%prod10_grc (nc) + safe( real(cncol(nc,nz, 8), r8) ) elseif (trim(species) == 'N') then - this%prod100_grc(nc) = this%prod100_grc(nc) + cncol(nc,nz,21)*CN_zone_weight(nz) - this%prod10_grc(nc) = this%prod10_grc(nc) + cncol(nc,nz,22)*CN_zone_weight(nz) + this%prod100_grc(nc) = this%prod100_grc(nc) + safe( real(cncol(nc,nz,21), r8) ) + this%prod10_grc (nc) = this%prod10_grc (nc) + safe( real(cncol(nc,nz,22), r8) ) else _ASSERT(.FALSE.,'unknown species') end if diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonFluxType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonFluxType.F90 index 642c87c1d..9244928b7 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonFluxType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonFluxType.F90 @@ -500,7 +500,21 @@ module CNVegCarbonFluxType __FILE__ contains - + pure elemental logical function valid(v) + use shr_kind_mod, only : r8 => shr_kind_r8 + use clm_varcon , only : spval + real(r8), intent(in) :: v + real(r8), parameter :: CF_FILL = 9.96921e36_r8 + real(r8), parameter :: BADHUGE = 1.0e20_r8 ! tighten from 1.0e30 + valid = (v == v) .and. (abs(v) < BADHUGE) .and. & + (abs(v - spval) > 1.0e-6_r8*max(1._r8,abs(spval))) .and. & + (abs(v - CF_FILL) > 1.0e-6_r8*max(1._r8,abs(CF_FILL))) + end function valid + + pure elemental real(r8) function safe(v) + real(r8), intent(in) :: v + safe = merge(v, 0._r8, valid(v)) + end function safe !--------------------------------------- subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft, carbon_type, cn5_cold_start, rc) @@ -1143,8 +1157,8 @@ subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft, carbon_type, cn5_co ! "new" variables: introduced in CNCLM50 if (cold_start.eqv..false.) then - this%annsum_litfall_patch(np) = cnpft(nc,nz,nv, 82) - this%tempsum_litfall_patch(np) = cnpft(nc,nz,nv, 83) + this%annsum_litfall_patch(np) = safe( real(cnpft(nc,nz,nv, 80), r8) ) + this%tempsum_litfall_patch(np) = safe( real(cnpft(nc,nz,nv, 81), r8) ) elseif (cold_start) then this%annsum_litfall_patch(np) = 0._r8 this%tempsum_litfall_patch(np) = 0._r8 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 index 0a06ce606..fe208bd7f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 @@ -224,6 +224,21 @@ module CNVegCarbonStateType __FILE__ contains + pure elemental logical function valid(v) + use shr_kind_mod, only : r8 => shr_kind_r8 + use clm_varcon , only : spval + real(r8), intent(in) :: v + real(r8), parameter :: CF_FILL = 9.96921e36_r8 + real(r8), parameter :: BADHUGE = 1.0e20_r8 ! tighten from 1.0e30 + valid = (v == v) .and. (abs(v) < BADHUGE) .and. & + (abs(v - spval) > 1.0e-6_r8*max(1._r8,abs(spval))) .and. & + (abs(v - CF_FILL) > 1.0e-6_r8*max(1._r8,abs(CF_FILL))) + end function valid + + pure elemental real(r8) function safe(v) + real(r8), intent(in) :: v + safe = merge(v, 0._r8, valid(v)) + end function safe !---------------------------------------------- subroutine Init(this, bounds, NLFilename, nch, ityp, fveg, cncol, cnpft) @@ -449,8 +464,8 @@ subroutine Init(this, bounds, NLFilename, nch, ityp, fveg, cncol, cnpft) allocate(this%deadstemc_col (begc:endc)) ; this%deadstemc_col (:) = nan allocate(this%fuelc_col (begc:endc)) ; this%fuelc_col (:) = nan allocate(this%fuelc_crop_col (begc:endc)) ; this%fuelc_crop_col (:) = nan - - allocate(this%totvegc_patch (begp:endp)) ; this%totvegc_patch (:) = spval + ! totvegc_patch set to zero so non visited patches don’t print CF_FILL. + allocate(this%totvegc_patch (begp:endp)) ; this%totvegc_patch (:) = 0._r8 allocate(this%totvegc_col (begc:endc)) ; this%totvegc_col (:) = nan allocate(this%totc_p2c_col (begc:endc)) ; this%totc_p2c_col (:) = nan @@ -469,7 +484,8 @@ subroutine Init(this, bounds, NLFilename, nch, ityp, fveg, cncol, cnpft) n = n + 1 this%totvegc_col(n) = cncol(nc,nz, 6) - this%seedc_grc (nc) = this%seedc_grc(nc) + cncol(nc,nz,9)*CN_zone_weight(nz) + this%seedc_grc(nc) = this%seedc_grc(nc) + safe( real(cncol(nc,nz, 9), r8) ) + this%totc_col (n) = cncol(nc,nz,14) do p = 0,numpft ! PFT index loop @@ -532,6 +548,8 @@ subroutine Init(this, bounds, NLFilename, nch, ityp, fveg, cncol, cnpft) end do ! nc call this%InitReadNML ( NLFilename ) + + ! --------------------------------------------------------------------------- end subroutine Init @@ -632,7 +650,12 @@ subroutine Summary_carbonstate(this, bounds, num_allc, filter_allc, & SHR_ASSERT_ALL_FL((ubound(soilbiogeochem_totsomc_col) == (/bounds%endc/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(soilbiogeochem_ctrunc_col) == (/bounds%endc/)), sourcefile, __LINE__) - ! calculate patch -level summary of carbon state + ! --- make non visited patches safe + this%dispvegc_patch(bounds%begp:bounds%endp) = 0._r8 + this%storvegc_patch(bounds%begp:bounds%endp) = 0._r8 + this%totvegc_patch(bounds%begp:bounds%endp) = 0._r8 + this%totc_patch (bounds%begp:bounds%endp) = 0._r8 + this%woodc_patch (bounds%begp:bounds%endp) = 0._r8 do fp = 1,num_soilp p = filter_soilp(fp) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 index 6354b3163..a513578ce 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 @@ -4,6 +4,7 @@ module CNVegNitrogenStateType use MAPL_ExceptionHandling use clm_varctl , only : use_matrixcn, use_crop use clm_varctl , only : use_nitrif_denitrif, use_vertsoilc, use_century_decomp + use clm_varctl , only : iulog use clm_varpar , only : NUM_ZON, NUM_VEG, VAR_COL, VAR_PFT, & numpft, CN_zone_weight use clm_varpar , only : ndecomp_cascade_transitions, ndecomp_pools, nlevcan @@ -213,6 +214,21 @@ module CNVegNitrogenStateType type(cnveg_nitrogenstate_type), public, target, save :: cnveg_nitrogenstate_inst contains + pure elemental logical function valid(v) + use shr_kind_mod, only : r8 => shr_kind_r8 + use clm_varcon , only : spval + real(r8), intent(in) :: v + real(r8), parameter :: CF_FILL = 9.96921e36_r8 + real(r8), parameter :: BADHUGE = 1.0e20_r8 ! tighten from 1.0e30 + valid = (v == v) .and. (abs(v) < BADHUGE) .and. & + (abs(v - spval) > 1.0e-6_r8*max(1._r8,abs(spval))) .and. & + (abs(v - CF_FILL) > 1.0e-6_r8*max(1._r8,abs(CF_FILL))) + end function valid + + pure elemental real(r8) function safe(v) + real(r8), intent(in) :: v + safe = merge(v, 0._r8, valid(v)) + end function safe !------------------------------------------------------------- subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft) @@ -445,35 +461,40 @@ subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft) do nz = 1,NUM_ZON ! CN zone loop n = n + 1 - this%seedn_grc(nc) = this%seedn_grc(nc) + cncol(nc,nz,23)*CN_zone_weight(nz) - this%totn_col(n) = cncol(nc,nz,29) + this%seedn_grc(nc) = this%seedn_grc(nc) + safe( real(cncol(nc,nz,23), r8) ) + + ! Slot 29 in CNCOL is an aggregate diagnostic - total column N and is where all junk shows up this avoids hard coding + this%totn_col(n) = spval + ! this%totn_col(n) = safe( real(cncol(nc,nz,29), r8) ) ! produces large values and if we do not wish to hard code + ! solution is NOT to ingest the restart aggregate; we'll anyway recompute in subroutine Summary_nitrogenstate at line: + ! this%totn_col(c) = this%totn_p2c_col(c) + cwdn + totlitn + totsomn + sminn + ntrunc do p = 0,numpft ! PFT index loop np = np + 1 do nv = 1,NUM_VEG ! defined veg loop if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>1.e-4) then - this%deadcrootn_patch (np) = cnpft(nc,nz,nv, 48) - this%deadcrootn_storage_patch (np) = cnpft(nc,nz,nv, 49) - this%deadcrootn_xfer_patch (np) = cnpft(nc,nz,nv, 50) - this%deadstemn_patch (np) = cnpft(nc,nz,nv, 51) - this%deadstemn_storage_patch (np) = cnpft(nc,nz,nv, 52) - this%deadstemn_xfer_patch (np) = cnpft(nc,nz,nv, 53) - this%frootn_patch (np) = cnpft(nc,nz,nv, 54) - this%frootn_storage_patch (np) = cnpft(nc,nz,nv, 55) - this%frootn_xfer_patch (np) = cnpft(nc,nz,nv, 56) - this%leafn_patch (np) = cnpft(nc,nz,nv, 57) - this%leafn_storage_patch (np) = cnpft(nc,nz,nv, 58) - this%leafn_xfer_patch (np) = cnpft(nc,nz,nv, 59) - this%livecrootn_patch (np) = cnpft(nc,nz,nv, 60) - this%livecrootn_storage_patch (np) = cnpft(nc,nz,nv, 61) - this%livecrootn_xfer_patch (np) = cnpft(nc,nz,nv, 62) - this%livestemn_patch (np) = cnpft(nc,nz,nv, 63) - this%livestemn_storage_patch (np) = cnpft(nc,nz,nv, 64) - this%livestemn_xfer_patch (np) = cnpft(nc,nz,nv, 65) - this%npool_patch (np) = cnpft(nc,nz,nv, 66) - this%ntrunc_patch (np) = cnpft(nc,nz,nv, 67) - this%retransn_patch (np) = cnpft(nc,nz,nv, 68) + this%deadcrootn_patch (np) = safe( real(cnpft(nc,nz,nv, 48), r8) ) + this%deadcrootn_storage_patch (np) = safe( real(cnpft(nc,nz,nv, 49), r8) ) + this%deadcrootn_xfer_patch (np) = safe( real(cnpft(nc,nz,nv, 50), r8) ) + this%deadstemn_patch (np) = safe( real(cnpft(nc,nz,nv, 51), r8) ) + this%deadstemn_storage_patch (np) = safe( real(cnpft(nc,nz,nv, 52), r8) ) + this%deadstemn_xfer_patch (np) = safe( real(cnpft(nc,nz,nv, 53), r8) ) + this%frootn_patch (np) = safe( real(cnpft(nc,nz,nv, 54), r8) ) + this%frootn_storage_patch (np) = safe( real(cnpft(nc,nz,nv, 55), r8) ) + this%frootn_xfer_patch (np) = safe( real(cnpft(nc,nz,nv, 56), r8) ) + this%leafn_patch (np) = safe( real(cnpft(nc,nz,nv, 57), r8) ) + this%leafn_storage_patch (np) = safe( real(cnpft(nc,nz,nv, 58), r8) ) + this%leafn_xfer_patch (np) = safe( real(cnpft(nc,nz,nv, 59), r8) ) + this%livecrootn_patch (np) = safe( real(cnpft(nc,nz,nv, 60), r8) ) + this%livecrootn_storage_patch (np) = safe( real(cnpft(nc,nz,nv, 61), r8) ) + this%livecrootn_xfer_patch (np) = safe( real(cnpft(nc,nz,nv, 62), r8) ) + this%livestemn_patch (np) = safe( real(cnpft(nc,nz,nv, 63), r8) ) + this%livestemn_storage_patch (np) = safe( real(cnpft(nc,nz,nv, 64), r8) ) + this%livestemn_xfer_patch (np) = safe( real(cnpft(nc,nz,nv, 65), r8) ) + this%npool_patch (np) = safe( real(cnpft(nc,nz,nv, 66), r8) ) + this%ntrunc_patch (np) = safe( real(cnpft(nc,nz,nv, 67), r8) ) + this%retransn_patch (np) = safe( real(cnpft(nc,nz,nv, 68), r8) ) end if end do !nv @@ -591,14 +612,12 @@ subroutine Summary_nitrogenstate(this, bounds, num_allc, filter_allc, & this%totvegn_col(c) ! total column nitrogen, including patch (TOTCOLN) - this%totn_col(c) = this%totn_p2c_col(c) + & soilbiogeochem_nitrogenstate_inst%cwdn_col(c) + & soilbiogeochem_nitrogenstate_inst%totlitn_col(c) + & soilbiogeochem_nitrogenstate_inst%totsomn_col(c) + & soilbiogeochem_nitrogenstate_inst%sminn_col(c) + & soilbiogeochem_nitrogenstate_inst%ntrunc_col(c) - end do end subroutine Summary_nitrogenstate diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegStateType.F90 index 6b911b12c..0fc447c9d 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegStateType.F90 @@ -5,6 +5,7 @@ module CNVegStateType use clm_varpar , only : nlevsno, nlevgrnd, nlevlak, nlevsoi, & num_zon, num_veg, var_col, var_pft, numpft use clm_varcon , only : spval, ispval + use clm_varctl , only : iulog use decompMod , only : bounds_type use AnnualFluxDribbler, only : annual_flux_dribbler_type, annual_flux_dribbler_patch diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 index b9b4c616c..a93a9c575 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 @@ -6,6 +6,7 @@ module CanopyStateType use clm_varpar , only : nlevcan, nvegwcs, numpft, num_zon, num_veg, & var_col, var_pft use clm_varcon , only : spval + use clm_varctl , only : iulog use nanMod , only : nan use decompMod , only : bounds_type use MAPL_ExceptionHandling @@ -167,7 +168,7 @@ subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft, cn5_cold_start, rc) ! "new" variables: introduced in CNCLM50 if (cold_start.eqv..false.) then do nw = 1,nvegwcs - this%vegwp_patch(np,nw) = cnpft(nc,nz,nv, 78+(nw-1)) + this%vegwp_patch(np,nw) = cnpft(nc,nz,nv, 76+(nw-1)) end do elseif (cold_start) then this%vegwp_patch(np,1:nvegwcs) = -2.5e4_r8 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PhotosynthesisMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PhotosynthesisMod.F90 index b1cc0f0c8..e7e7c4bc3 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PhotosynthesisMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PhotosynthesisMod.F90 @@ -379,17 +379,11 @@ subroutine Init(this,bounds,nch,ityp,fveg,cncol,cnpft,cn5_cold_start,rc) do nz = 1,num_zon ! CN zone loop do p = 0,numpft ! PFT index loop np = np + 1 - do nv = 1,num_veg ! defined veg loop - if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>1.e-4) then - if (cold_start) then - this%alphapsnsun_patch(np) = spval - this%alphapsnsha_patch(np) = spval - else if (cold_start.eqv..false.) then - this%alphapsnsun_patch(np) = cnpft(nc,nz,nv, 76) - this%alphapsnsha_patch(np) = cnpft(nc,nz,nv, 77) - end if - end if ! ityp =p - end do !nv + ! Using 0._r8 is safer than spval, because spval value can propagate if any path uses alpha before it’s overwritten. + ! If someday we do want alpha to persist across restarts, we need to: (1) pick free indices, (2) write them in CN_exit, and (3) read those same indices here. + ! Until then, don’t read them from cnpft ... we will use what we calculate. + this%alphapsnsun_patch(np) = 0._r8 + this%alphapsnsha_patch(np) = 0._r8 end do ! p end do ! nz end do ! nc diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 index 4d1f3057e..231ad6159 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 @@ -8,7 +8,7 @@ module SoilBiogeochemNitrogenStateType use clm_varpar , only : nlevdecomp_full, nlevdecomp, nlevsoi, & NUM_ZON, VAR_COL use clm_varcon , only : spval, dzsoi_decomp, zisoi - use clm_varctl , only : use_nitrif_denitrif, use_vertsoilc, use_century_decomp, use_soil_matrixcn + use clm_varctl , only : use_nitrif_denitrif, use_vertsoilc, use_century_decomp, use_soil_matrixcn, iulog use decompMod , only : bounds_type use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con use LandunitType , only : lun @@ -23,6 +23,9 @@ module SoilBiogeochemNitrogenStateType ! ! !PUBLIC MEMBER FUNCTIONS: + ! limit against corrupted restart values + real(r8), parameter :: N_MAX = 1.0e5_r8 ! any plausible pool, sanitaze monster values + type, public :: soilbiogeochem_nitrogenstate_type real(r8), pointer :: decomp_npools_vr_col (:,:,:) ! col (gN/m3) vertically-resolved decomposing (litter, cwd, soil) N pools @@ -71,6 +74,7 @@ module SoilBiogeochemNitrogenStateType ! type(sparse_matrix_type) :: AKXnacc ! col (gN/m3/yr) accumulated N transfers from j to i (col,i,j) per year in dimension(col,nlev*npools,nlev*npools) in sparse matrix type ! type(vector_type) :: matrix_Ninter ! col (gN/m3) vertically-resolved decomposing (litter, cwd, soil) N pools in dimension(col,nlev*npools) in vector type + contains procedure , public :: Summary @@ -85,6 +89,31 @@ module SoilBiogeochemNitrogenStateType contains + pure elemental logical function valid(v) + use shr_kind_mod, only : r8 => shr_kind_r8 + use clm_varcon , only : spval + real(r8), intent(in) :: v + real(r8), parameter :: CF_FILL = 9.96921e36_r8 + real(r8), parameter :: BADHUGE = 1.0e20_r8 ! tighten from 1.0e30 + valid = (v == v) .and. (abs(v) < BADHUGE) .and. & + (abs(v - spval) > 1.0e-6_r8*max(1._r8,abs(spval))) .and. & + (abs(v - CF_FILL) > 1.0e-6_r8*max(1._r8,abs(CF_FILL))) + end function valid + + pure elemental real(r8) function safe(v) + real(r8), intent(in) :: v + safe = merge(v, 0._r8, valid(v)) + end function safe + + pure elemental real(r8) function clamp_nonneg(v) result(out) + real(r8), intent(in) :: v + real(r8) :: x + x = safe(v) ! NaN/CF-fill/SPVAL -> 0 (via safe), else pass through + if (x < 0._r8) x = 0._r8 ! pools are nonnegative + if (x > N_MAX) x = 0._r8 ! treat absurd values as missing + out = x + end function clamp_nonneg + !------------------------------------------- subroutine Init(this, bounds, nch, cncol) @@ -162,31 +191,38 @@ subroutine Init(this, bounds, nch, cncol) ! initialize variables from restart file or set to cold start value n = 0 do nc = 1,nch ! catchment tile loop - do nz = 1,num_zon ! CN zone loop + do nz = 1,NUM_ZON ! CN zone loop n = n + 1 - - this%ntrunc_vr_col (n,1:nlevdecomp_full) = cncol(nc,nz,16) + + ! 1) ntrunc can be ± (sink), just scrub NaN/fill + this%ntrunc_vr_col(n,1:nlevdecomp_full) = safe( real(cncol(nc,nz,16), r8) ) + + ! 2) total mineral N ! jkolassa May 2022: for now nlevdecomp_full = 1; will need to add loop if we introduce more soil layers - this%sminn_vr_col (n,1:nlevdecomp_full) = cncol(nc,nz,24) - this%sminn_col (n) = this%sminn_vr_col(n,1) - + this%sminn_vr_col(n,1:nlevdecomp_full) = clamp_nonneg( real(cncol(nc,nz,24), r8) ) + this%sminn_col(n) = this%sminn_vr_col(n,1) + + ! 3) split into NO3/NH4 if (no_cn51_rst) then ! jkolassa Nov 2024: when no CN51 restart file is available compute NO3 and NH4 from N - this%smin_no3_col(n) = (1.25/2.25)*this%sminn_col(n) - this%smin_nh4_col(n) = this%sminn_col(n)/2.25 + ! derive split if no CN51 restart present + this%smin_no3_col(n) = clamp_nonneg( (1.25_r8/2.25_r8)*this%sminn_col(n) ) + this%smin_nh4_col(n) = clamp_nonneg( this%sminn_col(n)/2.25_r8 ) else - this%smin_no3_col(n) = cncol(nc,nz,36); - this%smin_nh4_col(n) = cncol(nc,nz,37); + this%smin_no3_col(n) = clamp_nonneg( real(cncol(nc,nz,36), r8) ) + this%smin_nh4_col(n) = clamp_nonneg( real(cncol(nc,nz,37), r8) ) end if - + + ! 4) mirror scalars into vertically-resolved arrays (nlevdecomp_full is 1 today, but this stays valid if >1 later) this%smin_no3_vr_col(n,1:nlevdecomp_full) = this%smin_no3_col(n) - this%smin_nh4_vr_col(n,1:nlevdecomp_full) = this%smin_nh4_col(n) + this%smin_nh4_vr_col(n,1:nlevdecomp_full) = this%smin_nh4_col(n) do np = 1,ndecomp_pools ! jkolassa May 2022: accounting for fact that pool order in CNCOL is different from CTSM - this%decomp_npools_col (n,np) = cncol(nc,nz,decomp_npool_cncol_index(np)) - this%decomp_npools_1m_col (n,np) = cncol(nc,nz,decomp_npool_cncol_index(np)) + this%decomp_npools_col (n,np) = clamp_nonneg( real(cncol(nc,nz,decomp_npool_cncol_index(np)), r8) ) + this%decomp_npools_1m_col(n,np) = this%decomp_npools_col(n,np) ! jkolassa May 2022: loop has to be added below if we add more biogeochemical (or soil) layers - this%decomp_npools_vr_col (n,1,np) = cncol(nc,nz,decomp_npool_cncol_index(np)) + this%decomp_npools_vr_col(n,1,np) = this%decomp_npools_col(n,np) + end do !np end do !nz end do @@ -205,9 +241,6 @@ subroutine Init(this, bounds, nch, cncol) end if end do - - - end subroutine Init !----------------------------------------------------------------------- diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SurfaceAlbedoType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SurfaceAlbedoType.F90 index 94531140d..1228e7cf4 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SurfaceAlbedoType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SurfaceAlbedoType.F90 @@ -72,7 +72,7 @@ module SurfaceAlbedoType contains !--------------------------------------------------- - subroutine Init(this, bounds, nch, cncol, cnpft) + subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft) ! !DESCRIPTION: ! Initialize CTSM surface albedo needed for calling CTSM routines @@ -86,12 +86,15 @@ subroutine Init(this, bounds, nch, cncol, cnpft) integer, intent(in) :: nch ! number of Catchment tiles real, dimension(nch,num_zon,var_col), intent(in) :: cncol ! column-level restart variable array real, dimension(nch,num_zon,num_veg,var_pft), intent(in) :: cnpft ! pft-level (patch-level) restart variable array + integer, dimension(nch,num_veg,num_zon), intent(in) :: ityp + real, dimension(nch,num_veg,num_zon), intent(in) :: fveg class(surfalb_type) :: this ! LOCAL integer :: begp, endp integer :: begc, endc integer :: np, nc, nz, p, nv, n + real, parameter :: fmin = 1.0e-4 ! Ignore vegetation fractions at or below this value !------------------------------- begp = bounds%begp ; endp = bounds%endp @@ -153,18 +156,17 @@ subroutine Init(this, bounds, nch, cncol, cnpft) this%nrad_patch(np) = 1 do nv = 1,num_veg ! defined veg loop - do n = 1,nlevcan - this%tlai_z_patch(np,n) = cnpft(nc,nz,nv, 73) - if (isnan(this%tlai_z_patch(np,n))) then - this%tlai_z_patch(np,n) = 0. - end if - this%tsai_z_patch(np,n) = cnpft(nc,nz,nv, 74) - if (isnan(this%tsai_z_patch(np,n))) then - this%tsai_z_patch(np,n) = 0. - end if - this%vcmaxcintsha_patch(np) = 1._r8 - this%vcmaxcintsun_patch(np) = 1._r8 - end do + if (ityp(nc,nv,nz) == p .and. fveg(nc,nv,nz) > fmin) then + ! without this if loop we are overwriting each patch multiple times because we loop nv=1:num_veg without + ! checking which nv corresponds to the active PFT for patch p. + do n = 1, nlevcan + this%tlai_z_patch(np,n) = cnpft(nc,nz,nv,73); if (isnan(this%tlai_z_patch(np,n))) this%tlai_z_patch(np,n)=0._r8 + this%tsai_z_patch(np,n) = cnpft(nc,nz,nv,74); if (isnan(this%tsai_z_patch(np,n))) this%tsai_z_patch(np,n)=0._r8 + end do + this%vcmaxcintsun_patch(np) = 1._r8 + this%vcmaxcintsha_patch(np) = 1._r8 + exit ! <- stop after filling the matching nv for this patch + end if ! ityp =p end do !nv end do ! p end do ! nz diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 index 1d9c6409d..e6dcd08a4 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 @@ -98,7 +98,6 @@ module CatchmentCNRstMod real, allocatable :: runsurf(:) contains - procedure :: write_nc4 procedure :: allocate_cn procedure :: add_bcs_to_cnrst @@ -126,6 +125,11 @@ function CatchmentCNRst_create(filename, cnclm, time, rc) result (catch) character(len=:), pointer :: dname type(FileMetadata) :: meta character(len=256) :: Iam = "CatchmentCNRst_create" + integer :: nt_col, blk_col, off + real :: N_MAX + integer :: joff + integer :: n_fixed + call MAPL_NCIOGetFileType(filename, filetype, __RC__) if (filetype /= 0) then @@ -224,11 +228,61 @@ function CatchmentCNRst_create(filename, cnclm, time, rc) result (catch) call MAPL_VarRead( formatter, "TPREC10D", catch%TPREC10D, __RC__) call MAPL_VarRead( formatter, "TPREC60D", catch%TPREC60D, __RC__) call MAPL_VarRead( formatter, "ET365D", catch%ET365D , __RC__) ! rreichle, bug fix 4 Aug 2025: added for CLM51 - call MAPL_VarRead( formatter, "RUNSURF", catch%ET365D , __RC__) ! rreichle, bug fix 4 Aug 2025: added for CLM51 + call MAPL_VarRead( formatter, "RUNSURF", catch%RUNSURF , __RC__) ! rreichle, bug fix 4 Aug 2025: added for CLM51 - endif + endif !! <— close CLM40/CLM51 branch - endif + ! strip NaNs + where (.not.(catch%CNCOL == catch%CNCOL)) catch%CNCOL = 0.0 + where (.not.(catch%CNPFT == catch%CNPFT)) catch%CNPFT = 0.0 + + ! Targeted sanitize for N in CNCOL: 24=sminn, 36=NO3, 37=NH4 + ! Without this section model can't run values are junk! + + N_MAX = 1.0e5 ! conservative guardrail + + ! Slot 24: sminn + joff = (24-1)*nzone + where (catch%CNCOL(:, joff+1:joff+nzone) /= catch%CNCOL(:, joff+1:joff+nzone)) + catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + end where + where (catch%CNCOL(:, joff+1:joff+nzone) < 0.0) + catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + end where + where (catch%CNCOL(:, joff+1:joff+nzone) > N_MAX) + catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + end where + + ! Slot 36: smin_no3 + joff = (36-1)*nzone + where (catch%CNCOL(:, joff+1:joff+nzone) /= catch%CNCOL(:, joff+1:joff+nzone)) + catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + end where + where (catch%CNCOL(:, joff+1:joff+nzone) < 0.0) + catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + end where + where (catch%CNCOL(:, joff+1:joff+nzone) > N_MAX) + catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + end where + + ! Slot 37: smin_nh4 + joff = (37-1)*nzone + where (catch%CNCOL(:, joff+1:joff+nzone) /= catch%CNCOL(:, joff+1:joff+nzone)) + catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + end where + where (catch%CNCOL(:, joff+1:joff+nzone) < 0.0) + catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + end where + where (catch%CNCOL(:, joff+1:joff+nzone) > N_MAX) + catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + end where + + n_fixed = count(.not.(catch%CNCOL(:, (24-1)*nzone+1:(24)*nzone) == catch%CNCOL(:, (24-1)*nzone+1:(24)*nzone))) & + + count(catch%CNCOL(:, (24-1)*nzone+1:(24)*nzone) < 0.0) & + + count(catch%CNCOL(:, (24-1)*nzone+1:(24)*nzone) > N_MAX) + if (n_fixed > 0) write(*,*) 'RSTCHK sanitize CNCOL(24): fixed=', n_fixed + + endif !! <— closes myid call formatter%close() @@ -355,7 +409,8 @@ subroutine write_nc4(this, filename, rc) call MAPL_VarWrite( formatter, "TPREC10D", this%TPREC10D, __RC__) call MAPL_VarWrite( formatter, "TPREC60D", this%TPREC60D, __RC__) call MAPL_VarWrite( formatter, "ET365D", this%ET365D , __RC__) - call MAPL_VarWrite( formatter, "RUNSURF", this%ET365D , __RC__) + call MAPL_VarWrite( formatter, "RUNSURF", this%RUNSURF , __RC__) + endif @@ -1114,8 +1169,7 @@ subroutine re_tile(this, InTileFile, OutBcsDir, OutTileFile, surflay, rc) end do end do - where(isnan(var_off_pft)) var_off_pft = 0. - where(var_off_pft /= var_off_pft) var_off_pft = 0. ! should do same as previous line, but keep for now just in case + print *, 'calculating regridded carbn' @@ -1169,7 +1223,19 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & integer :: N, STATUS, nv, nx, offl_cell, ityp_new, i, j, nz, iv, jj real :: fveg_new character(256) :: Iam = "write_regridded_carbon" - + ! Column slots (CNCOL) that must be >= 0 + integer, parameter :: nonneg_col(*) = (/ 1,2,3,4,5, 6, 10,11,12,13, 14,15,16, 21,22,23,24,25,26,27,28,29, 30,31,32,33,34,35,36,37 /) + ! PFT slots (CNPFT) that must be >= 0 (DO NOT include 76–79: veg water potentials) + integer, parameter :: nonneg_pft(*) = (/ & + 1,2,3,4,5, 6,7,8,9,10, 11,12,13,14,15, 16,17,18,19,20, 21,22,23, & + 24,25,26,27, & + 41,42, & ! prev_*_to_litter + 45,47, & ! tempsum_npp, xsmrpool_recover + 69,70,71,72,73,74, & ! canopy geometry + 75, & ! plant_ndemand + 80,81 & ! litfall sums + /) + allocate (CLMC_pf1(NTILES)) allocate (CLMC_pf2(NTILES)) @@ -1205,7 +1271,8 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & allocate (var_pft_out (1: NTILES, 1 : nzone,1 : nveg, 1 : var_pft)) var_col_out = 0. - var_pft_out = NaN + !var_pft_out = NaN + var_pft_out = 0.0 OUT_TILE : DO N = 1, NTILES @@ -1397,27 +1464,47 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & var_pft_out(n, nz,nv,68) = max(var_pft_out(n, nz,nv,68),0.) var_pft_out(n, nz,nv,69) = max(var_pft_out(n, nz,nv,69),0.) var_pft_out(n, nz,nv,70) = max(var_pft_out(n, nz,nv,70),0.) + var_pft_out(n, nz,nv,71) = max(var_pft_out(n, nz,nv,71),0.) ! HBOT are heights + var_pft_out(n, nz,nv,72) = max(var_pft_out(n, nz,nv,72),0.) ! HTOP are heights var_pft_out(n, nz,nv,73) = max(var_pft_out(n, nz,nv,73),0.) var_pft_out(n, nz,nv,74) = max(var_pft_out(n, nz,nv,74),0.) !if(this%isCLM45) var_pft_out(n, nz,nv,75) = max(var_pft_out(n, nz,nv,75),0.) if(this%isCLM51) then var_pft_out(n, nz,nv,75) = max(var_pft_out(n, nz,nv,75),0.) - var_pft_out(n, nz,nv,76) = max(var_pft_out(n, nz,nv,76),0.) - var_pft_out(n, nz,nv,77) = max(var_pft_out(n, nz,nv,77),0.) - var_pft_out(n, nz,nv,78) = max(var_pft_out(n, nz,nv,78),0.) - var_pft_out(n, nz,nv,79) = max(var_pft_out(n, nz,nv,79),0.) + !! These indices are veg water potentials and can be negative by design.. !cnpft(:, :, :, 76..79) => vegwp_patch(1..4) + !These are pressures, not masses. + !Plants pull water under tension, so the numbers are negative ( example −0.3 to −3 MPa, or negative Pascals). + !var_pft_out(n, nz,nv,76) = max(var_pft_out(n, nz,nv,76),0.) + !var_pft_out(n, nz,nv,77) = max(var_pft_out(n, nz,nv,77),0.) + !var_pft_out(n, nz,nv,78) = max(var_pft_out(n, nz,nv,78),0.) + !var_pft_out(n, nz,nv,79) = max(var_pft_out(n, nz,nv,79),0.) var_pft_out(n, nz,nv,80) = max(var_pft_out(n, nz,nv,80),0.) var_pft_out(n, nz,nv,81) = max(var_pft_out(n, nz,nv,81),0.) - var_pft_out(n, nz,nv,82) = max(var_pft_out(n, nz,nv,82),0.) - var_pft_out(n, nz,nv,83) = max(var_pft_out(n, nz,nv,83),0.) + ! imaginary... 82/83 !! + !var_pft_out(n, nz,nv,82) = max(var_pft_out(n, nz,nv,82),0.) + !var_pft_out(n, nz,nv,83) = max(var_pft_out(n, nz,nv,83),0.) end if endif end do NVLOOP3 ! end veg loop endif ! end carbon check end do NZLOOP ! end zone loop - ! Update dayx variable var_pft_out (:,:,28) + ! --- Enforce non-negativity on known non-negative fields (all tiles) --- + do jj = 1, size(nonneg_col) + do nz = 1, nzone + var_col_out(:,nz,nonneg_col(jj)) = max( var_col_out(:,nz,nonneg_col(jj)), 0.0 ) + end do + end do + + do j = 1, size(nonneg_pft) + do nz = 1, nzone + do nv = 1, nveg + var_pft_out(:,nz,nv,nonneg_pft(j)) = max( var_pft_out(:,nz,nv,nonneg_pft(j)), 0.0 ) + end do + end do + end do + ! Update dayx variable var_pft_out (:,:,28) do j = 28, 28 ! 1,VAR_PFT var_pft_out (:,:,:,28) do nv = 1,nveg do nz = 1,nzone @@ -1557,6 +1644,11 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & end do OUT_TILE + ! --- Final scrub: only finite, plausible magnitudes make it to the restart + where (var_col_out /= var_col_out .or. abs(var_col_out) > 1.0e20) var_col_out = 0.0 + where (var_pft_out /= var_pft_out .or. abs(var_pft_out) > 1.0e20) var_pft_out = 0.0 + + i = 1 deallocate(this%cncol) allocate(this%cncol(NTILES, nzone*VAR_COL)) From 4b321a21d929df74a9d3ab648c85a80900047a72 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 10 Sep 2025 17:40:45 -0400 Subject: [PATCH 02/29] for reference, added commented version of Biljana's init of CNCOL, CNPFT to zero, which got lost in merge of 0-diff changes from Jana's branch (CNCLM_DriverMod.F90) --- .../GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 index 53aff7c07..0367f3131 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 @@ -62,8 +62,9 @@ module CNCLM_DriverMod public :: get_CN_LAI contains - -!--------------------------------- + + !--------------------------------- + subroutine CN_Driver(istep,nch,ityp,fveg,ndep,tp1,tairm,psis,bee,dayl,btran_fire,car1m,& rzm,sfm,rhm,windm,rainfm,snowfm,prec10d,prec60d,et365d,gdp,& abm,peatf,hdm,lnfm,poros,rh30,totwat,bflow,runsrf,sndzn,& @@ -592,6 +593,10 @@ subroutine CN_exit(nch,ityp,fveg,cncol,cnpft) !---------------- + !rrXbo10Sep2025 ! Clean slate: zero arrays so any unassigned slots are safe in the restart + !rrXbo10Sep2025 cncol(:,:,:) = 0.0 + !rrXbo10Sep2025 cnpft(:,:,:,:) = 0.0 + n = 0 np = 0 From 4d0ef802929def4f8d12846ed989d0f7afa861b1 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 10 Sep 2025 17:50:35 -0400 Subject: [PATCH 03/29] turned Biljana's "safe" function into identity operator (CNProductsMod.F90, CNVegCarbonFluxType.F90, CNVegCarbonStateType.F90, CNVegNitrogenStateType.F90, SoilBiogeochemNitrogenStateType.F90) --- .../GEOScatchCNCLM51_GridComp/CLM51/CNProductsMod.F90 | 2 +- .../GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonFluxType.F90 | 2 +- .../GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 | 2 +- .../GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 | 2 +- .../CLM51/SoilBiogeochemNitrogenStateType.F90 | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNProductsMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNProductsMod.F90 index 96ac1813c..ae6fe7f0c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNProductsMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNProductsMod.F90 @@ -87,7 +87,7 @@ end function valid pure elemental real(r8) function safe(v) real(r8), intent(in) :: v - safe = merge(v, 0._r8, valid(v)) + safe = merge(v, 0._r8, .true.) ! valid(v)) !rrXbo 10Sep2025 end function safe !-------------------------------------------------------------- subroutine Init(this, bounds, nch, cncol, species, rc) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonFluxType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonFluxType.F90 index 9244928b7..e86e9b0c3 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonFluxType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonFluxType.F90 @@ -513,7 +513,7 @@ end function valid pure elemental real(r8) function safe(v) real(r8), intent(in) :: v - safe = merge(v, 0._r8, valid(v)) + safe = merge(v, 0._r8, .true.) ! valid(v)) !rrXbo 10Sep2025 end function safe !--------------------------------------- subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft, carbon_type, cn5_cold_start, rc) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 index fe208bd7f..ce25a1e4b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 @@ -237,7 +237,7 @@ end function valid pure elemental real(r8) function safe(v) real(r8), intent(in) :: v - safe = merge(v, 0._r8, valid(v)) + safe = merge(v, 0._r8, .true.) ! valid(v)) !rrXbo 10Sep2025 end function safe !---------------------------------------------- diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 index a513578ce..30f571b8e 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 @@ -227,7 +227,7 @@ end function valid pure elemental real(r8) function safe(v) real(r8), intent(in) :: v - safe = merge(v, 0._r8, valid(v)) + safe = merge(v, 0._r8, .true.) ! valid(v)) !rrXbo 10Sep2025 end function safe !------------------------------------------------------------- diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 index 231ad6159..c41a9d3f0 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 @@ -102,7 +102,7 @@ end function valid pure elemental real(r8) function safe(v) real(r8), intent(in) :: v - safe = merge(v, 0._r8, valid(v)) + safe = merge(v, 0._r8, .true.) ! valid(v)) !rrXbo 10Sep2025 end function safe pure elemental real(r8) function clamp_nonneg(v) result(out) From d803329e931bb9fe3b4f07e6610438e94e94f4b1 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 10 Sep 2025 18:05:37 -0400 Subject: [PATCH 04/29] reverted unwanted changes from Biljana's commit (CNVegCarbonStateType.F90, CNVegNitrogenStateType.F90) --- .../CLM51/CNVegCarbonStateType.F90 | 4 ++-- .../CLM51/CNVegNitrogenStateType.F90 | 12 +++++++----- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 index ce25a1e4b..980ab59a6 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 @@ -464,8 +464,8 @@ subroutine Init(this, bounds, NLFilename, nch, ityp, fveg, cncol, cnpft) allocate(this%deadstemc_col (begc:endc)) ; this%deadstemc_col (:) = nan allocate(this%fuelc_col (begc:endc)) ; this%fuelc_col (:) = nan allocate(this%fuelc_crop_col (begc:endc)) ; this%fuelc_crop_col (:) = nan - ! totvegc_patch set to zero so non visited patches don’t print CF_FILL. - allocate(this%totvegc_patch (begp:endp)) ; this%totvegc_patch (:) = 0._r8 + ! totvegc_patch set to zero so non visited patches don’t print CF_FILL. !rrXbo 10Sep2025 + allocate(this%totvegc_patch (begp:endp)) ; this%totvegc_patch (:) = spval ! 0._r8 !rrXbo 10Sep2025 allocate(this%totvegc_col (begc:endc)) ; this%totvegc_col (:) = nan allocate(this%totc_p2c_col (begc:endc)) ; this%totc_p2c_col (:) = nan diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 index 30f571b8e..52225cc75 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 @@ -463,11 +463,13 @@ subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft) this%seedn_grc(nc) = this%seedn_grc(nc) + safe( real(cncol(nc,nz,23), r8) ) - ! Slot 29 in CNCOL is an aggregate diagnostic - total column N and is where all junk shows up this avoids hard coding - this%totn_col(n) = spval - ! this%totn_col(n) = safe( real(cncol(nc,nz,29), r8) ) ! produces large values and if we do not wish to hard code - ! solution is NOT to ingest the restart aggregate; we'll anyway recompute in subroutine Summary_nitrogenstate at line: - ! this%totn_col(c) = this%totn_p2c_col(c) + cwdn + totlitn + totsomn + sminn + ntrunc + this%totn_col(n) = cncol(nc,nz,29) ! borescan, reichle, 10 Sep 2025: may not be needed in restart file, diagnosed in subroutine Summary_nitrogenstate() (?) + + !rrXbo 10Sep2025 ! Slot 29 in CNCOL is an aggregate diagnostic - total column N and is where all junk shows up this avoids hard coding + !rrXbo 10Sep2025 this%totn_col(n) = spval + !rrXbo 10Sep2025 ! this%totn_col(n) = safe( real(cncol(nc,nz,29), r8) ) ! produces large values and if we do not wish to hard code + !rrXbo 10Sep2025 ! solution is NOT to ingest the restart aggregate; we'll anyway recompute in subroutine Summary_nitrogenstate at line: + !rrXbo 10Sep2025 ! this%totn_col(c) = this%totn_p2c_col(c) + cwdn + totlitn + totsomn + sminn + ntrunc do p = 0,numpft ! PFT index loop np = np + 1 From 8c01247046ac0c7b6062822c5ec72aed6e15d76d Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 10 Sep 2025 18:09:47 -0400 Subject: [PATCH 05/29] reverted more of Biljana's changes (CNVegCarbonStateType.F90) --- .../CLM51/CNVegCarbonStateType.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 index 980ab59a6..414c4be76 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 @@ -650,12 +650,12 @@ subroutine Summary_carbonstate(this, bounds, num_allc, filter_allc, & SHR_ASSERT_ALL_FL((ubound(soilbiogeochem_totsomc_col) == (/bounds%endc/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(soilbiogeochem_ctrunc_col) == (/bounds%endc/)), sourcefile, __LINE__) - ! --- make non visited patches safe - this%dispvegc_patch(bounds%begp:bounds%endp) = 0._r8 - this%storvegc_patch(bounds%begp:bounds%endp) = 0._r8 - this%totvegc_patch(bounds%begp:bounds%endp) = 0._r8 - this%totc_patch (bounds%begp:bounds%endp) = 0._r8 - this%woodc_patch (bounds%begp:bounds%endp) = 0._r8 + !rrXbo 10Sep2025 ! --- make non visited patches safe + !rrXbo 10Sep2025 this%dispvegc_patch(bounds%begp:bounds%endp) = 0._r8 + !rrXbo 10Sep2025 this%storvegc_patch(bounds%begp:bounds%endp) = 0._r8 + !rrXbo 10Sep2025 this%totvegc_patch(bounds%begp:bounds%endp) = 0._r8 + !rrXbo 10Sep2025 this%totc_patch (bounds%begp:bounds%endp) = 0._r8 + !rrXbo 10Sep2025 this%woodc_patch (bounds%begp:bounds%endp) = 0._r8 do fp = 1,num_soilp p = filter_soilp(fp) From 3feb4cf55012dc80d66ad20ef1a529c6831b8beb Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 10 Sep 2025 18:13:05 -0400 Subject: [PATCH 06/29] removed obsolete "use iulog" statements (CNVegNitrogenStateType.F90, CNVegStateType.F90, CanopyStateType.F90, SoilBiogeochemNitrogenStateType.F90) --- .../GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 | 2 +- .../GEOScatchCNCLM51_GridComp/CLM51/CNVegStateType.F90 | 2 +- .../GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 | 2 +- .../CLM51/SoilBiogeochemNitrogenStateType.F90 | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 index 52225cc75..f314b5563 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 @@ -4,7 +4,7 @@ module CNVegNitrogenStateType use MAPL_ExceptionHandling use clm_varctl , only : use_matrixcn, use_crop use clm_varctl , only : use_nitrif_denitrif, use_vertsoilc, use_century_decomp - use clm_varctl , only : iulog + !use clm_varctl , only : iulog !rrXbo 10Sep2025 use clm_varpar , only : NUM_ZON, NUM_VEG, VAR_COL, VAR_PFT, & numpft, CN_zone_weight use clm_varpar , only : ndecomp_cascade_transitions, ndecomp_pools, nlevcan diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegStateType.F90 index 0fc447c9d..52e5144af 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegStateType.F90 @@ -5,7 +5,7 @@ module CNVegStateType use clm_varpar , only : nlevsno, nlevgrnd, nlevlak, nlevsoi, & num_zon, num_veg, var_col, var_pft, numpft use clm_varcon , only : spval, ispval - use clm_varctl , only : iulog + !use clm_varctl , only : iulog !rrXbo 10Sep2025 use decompMod , only : bounds_type use AnnualFluxDribbler, only : annual_flux_dribbler_type, annual_flux_dribbler_patch diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 index a93a9c575..b075186bd 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 @@ -6,7 +6,7 @@ module CanopyStateType use clm_varpar , only : nlevcan, nvegwcs, numpft, num_zon, num_veg, & var_col, var_pft use clm_varcon , only : spval - use clm_varctl , only : iulog + !use clm_varctl , only : iulog !rrXbo 10Sep2025 use nanMod , only : nan use decompMod , only : bounds_type use MAPL_ExceptionHandling diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 index c41a9d3f0..f31b2ff93 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 @@ -8,7 +8,7 @@ module SoilBiogeochemNitrogenStateType use clm_varpar , only : nlevdecomp_full, nlevdecomp, nlevsoi, & NUM_ZON, VAR_COL use clm_varcon , only : spval, dzsoi_decomp, zisoi - use clm_varctl , only : use_nitrif_denitrif, use_vertsoilc, use_century_decomp, use_soil_matrixcn, iulog + use clm_varctl , only : use_nitrif_denitrif, use_vertsoilc, use_century_decomp, use_soil_matrixcn !, iulog !rrXbo 10Sep2025 use decompMod , only : bounds_type use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con use LandunitType , only : lun From 0d586920fb35849e342d2160dbfbaa695a1a25e5 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 11 Sep 2025 09:51:13 -0400 Subject: [PATCH 07/29] fix comment about PFTs in CatchCN51 (GEOS_CatchCNCLM51GridComp.F90) --- .../GEOS_CatchCNCLM51GridComp.F90 | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 index 18a8a58eb..fd13ae4df 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 @@ -132,7 +132,10 @@ module GEOS_CatchCNCLM51GridCompMod ! ROOTL import from GEOS_VegdynGridComp was disabled and brought the look up table ! in order to obtain ROOTL for primary and secondary types. + + ! map catchment type into PFT +! NOTE: Unlike CatchCN40/45, CatchCN51 no longer uses split PFTs (depending on moisture stress) ! --------------------------- !PFT Description ! 0 bare @@ -145,16 +148,12 @@ module GEOS_CatchCNCLM51GridCompMod ! 7 broadleaf deciduous temperate tree ! 8 broadleaf deciduous boreal tree ! 9 broadleaf evergreen temperate shrub -! 10 broadleaf deciduous temperate shrub [moisture + deciduous] -! 11 broadleaf deciduous temperate shrub [moisture stress only] -! 12 broadleaf deciduous boreal shrub -! 13 arctic c3 grass -! 14 cool c3 grass [moisture + deciduous] -! 15 cool c3 grass [moisture stress only] -! 16 warm c4 grass [moisture + deciduous] -! 17 warm c4 grass [moisture stress only] -! 18 crop [moisture + deciduous] -! 19 crop [moisture stress only] +! 10 broadleaf deciduous temperate shrub +! 11 broadleaf deciduous boreal shrub +! 12 arctic c3 grass +! 13 cool c3 grass +! 14 warm c4 grass +! 15 crop ! Catchment types and PFT mapping: ! From 69f46b5263f4d6ceacade3729ecce563e34e8175 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 11 Sep 2025 12:58:19 -0400 Subject: [PATCH 08/29] use FVEG_MIN parameter instead of hard-coded 1.e-4 --- .../CLM51/CNCLM_DriverMod.F90 | 12 +++--- .../CLM51/CNCLM_Photosynthesis.F90 | 6 +-- .../CLM51/CNVegCarbonFluxType.F90 | 5 +-- .../CLM51/CNVegCarbonStateType.F90 | 4 +- .../CLM51/CNVegNitrogenFluxType.F90 | 7 ++-- .../CLM51/CNVegNitrogenStateType.F90 | 6 +-- .../CLM51/CNVegStateType.F90 | 9 ++--- .../CLM51/CanopyStateType.F90 | 7 ++-- .../CLM51/PatchType.F90 | 6 +-- .../CLM51/SoilBiogeochemStateType.F90 | 7 ++-- .../CLM51/SurfaceAlbedoType.F90 | 5 +-- .../CLM51/clm_varpar.F90 | 13 ++++--- .../CLM51/filterMod.F90 | 39 +++++++++---------- .../GEOS_CatchCNCLM51GridComp.F90 | 7 ++-- .../Shared/clm_varpar_shared.F90 | 30 +++++++------- 15 files changed, 81 insertions(+), 82 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 index 0367f3131..fac8b19ca 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 @@ -4,8 +4,8 @@ module CNCLM_DriverMod use nanMod , only : nan use CNVegetationFacade - use clm_varpar , only : nlevsno, nlevmaxurbgrnd, num_veg, num_zon, CN_zone_weight,& - var_col, var_pft, nlevgrnd, numpft, ndecomp_pools + use clm_varpar , only : nlevsno, nlevmaxurbgrnd, num_veg, num_zon, CN_zone_weight, & + var_col, var_pft, nlevgrnd, numpft, ndecomp_pools, FVEG_MIN use clm_varcon , only : grav, denh2o use clm_time_manager , only : is_first_step, get_nstep use decompMod @@ -319,7 +319,7 @@ subroutine CN_Driver(istep,nch,ityp,fveg,ndep,tp1,tairm,psis,bee,dayl,btran_fire do nv = 1,num_veg ! defined veg loop - if(ityp(nc,nv,nz)==np .and. fveg(nc,nv,nz)>1.e-4) then + if(ityp(nc,nv,nz)==np .and. fveg(nc,nv,nz)>FVEG_MIN) then photosyns_inst%psnsun_patch(p) = psnsunm(nc,nv,nz) photosyns_inst%psnsha_patch(p) = psnsham(nc,nv,nz) @@ -507,7 +507,7 @@ subroutine CN_Driver(istep,nch,ityp,fveg,ndep,tp1,tairm,psis,bee,dayl,btran_fire do nv = 1,num_veg ! defined veg loop - if(ityp(nc,nv,nz)==np .and. fveg(nc,nv,nz)>1.e-4) then + if(ityp(nc,nv,nz)==np .and. fveg(nc,nv,nz)>FVEG_MIN) then zlai(nc,nv,nz) = canopystate_inst%elai_patch(p) @@ -649,7 +649,7 @@ subroutine CN_exit(nch,ityp,fveg,cncol,cnpft) do nv = 1,num_veg ! defined veg loop - if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>1.e-4) then + if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>FVEG_MIN) then cnpft(nc,nz,nv, 1) = bgc_vegetation_inst%cnveg_carbonstate_inst%cpool_patch (np ) cnpft(nc,nz,nv, 2) = bgc_vegetation_inst%cnveg_carbonstate_inst%deadcrootc_patch (np ) @@ -791,7 +791,7 @@ subroutine get_CN_LAI(nch,ityp,fveg,elai,esai,tlai,tsai) ! extract LAI & SAI from CN clmtype ! --------------------------------- - if(ityp(nc,nv,nz)==p .and. ityp(nc,nv,nz)>0 .and. fveg(nc,nv,nz)>1.e-4) then + if(ityp(nc,nv,nz)==p .and. ityp(nc,nv,nz)>0 .and. fveg(nc,nv,nz)>FVEG_MIN) then elai(nc,nv,nz) = elai_clm(np) if(present(esai)) esai(nc,nv,nz) = esai_clm(np) if(present(tlai)) tlai(nc,nv,nz) = tlai_clm(np) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_Photosynthesis.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_Photosynthesis.F90 index 36e752e47..e50502a37 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_Photosynthesis.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_Photosynthesis.F90 @@ -1,7 +1,7 @@ module CNCLM_Photosynthesis use MAPL_ConstantsMod - use clm_varpar, only : numpft, numrad, num_veg, num_zon, & + use clm_varpar, only : numpft, numrad, num_veg, num_zon, FVEG_MIN, & nlevcan use decompMod use PatchType @@ -23,7 +23,7 @@ module CNCLM_Photosynthesis use WaterStateType use WaterType use CNVegetationFacade - + implicit none private @@ -459,7 +459,7 @@ subroutine catchcn_calc_rc(nch,fveg,tc,qa,pbot,co2v,dayl_factor, & do p = 0,numpft ! PFT index loop np = np + 1 do nv = 1,num_veg ! defined veg loop - if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>1.e-4) then + if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>FVEG_MIN) then ! stomatal resistances rs = laisun(np)/max(rssun(np), 1.e-06_r8 ) + laisha(np)/max(rssha(np), 1.e-06_r8 ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonFluxType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonFluxType.F90 index e86e9b0c3..0884c6c3d 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonFluxType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonFluxType.F90 @@ -17,7 +17,7 @@ module CNVegCarbonFluxType ideadcroot,ideadcroot_st,ideadcroot_xf,& igrain,igrain_st,igrain_xf,ioutc use clm_varpar , only : numpft, num_zon, num_veg, & - var_col, var_pft, CN_zone_weight + var_col, var_pft, CN_zone_weight, FVEG_MIN use clm_varctl , only : use_crop, use_matrixcn, use_cndv, use_grainproduct, iulog use clm_varcon , only : dzsoi_decomp use pftconMod , only : npcropmin @@ -30,7 +30,6 @@ module CNVegCarbonFluxType use shr_log_mod , only : errMsg => shr_log_errMsg ! use SPMMod , only : sparse_matrix_type, diag_matrix_type, vector_type - ! !PUBLIC TYPES: implicit none save @@ -1140,7 +1139,7 @@ subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft, carbon_type, cn5_co do p = 0,numpft ! PFT index loop np = np + 1 do nv = 1,num_veg ! defined veg loop - if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>1.e-4) then + if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>FVEG_MIN) then this%gpp_patch(p) = 0._r8 this%excess_cflux_patch(p) = 0._r8 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 index 414c4be76..699cc88a8 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegCarbonStateType.F90 @@ -5,7 +5,7 @@ module CNVegCarbonStateType use shr_kind_mod , only : r8 => shr_kind_r8 use clm_varctl , only : iulog, use_cndv, use_crop, use_matrixcn use clm_varpar , only : numpft, num_zon, num_veg, & - var_col, var_pft, CN_zone_weight + var_col, var_pft, CN_zone_weight, FVEG_MIN use clm_varcon , only : spval use nanMod , only : nan use decompMod , only : bounds_type @@ -491,7 +491,7 @@ subroutine Init(this, bounds, NLFilename, nch, ityp, fveg, cncol, cnpft) do p = 0,numpft ! PFT index loop np = np + 1 do nv = 1,num_veg ! defined veg loop - if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>1.e-4) then + if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>FVEG_MIN) then ! "old" variables: CNCLM45 and before this%cpool_patch (np) = cnpft(nc,nz,nv, 1) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenFluxType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenFluxType.F90 index 3710e26b1..23c40bb36 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenFluxType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenFluxType.F90 @@ -15,15 +15,14 @@ module CNVegNitrogenFluxType ideadcroot,ideadcroot_st,ideadcroot_xf,& igrain,igrain_st,igrain_xf,iretransn,ioutn use clm_varpar , only : numpft, num_zon, num_veg, & - var_col, var_pft, CN_zone_weight + var_col, var_pft, CN_zone_weight, FVEG_MIN use clm_varcon , only : spval, ispval, dzsoi_decomp use clm_varctl , only : use_nitrif_denitrif, use_vertsoilc, use_crop, use_matrixcn use PatchType , only : patch use CNSharedParamsMod , only : use_fun use LandunitType , only : lun use landunit_varcon , only : istsoil, istcrop - - + ! !PUBLIC TYPES: implicit none save @@ -983,7 +982,7 @@ subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft) do p = 0,numpft ! PFT index loop np = np + 1 do nv = 1,num_veg ! defined veg loop - if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>1.e-4) then + if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>FVEG_MIN) then this%plant_ndemand_patch (np) = cnpft(nc,nz,nv, 75) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 index f314b5563..9db1476be 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegNitrogenStateType.F90 @@ -6,7 +6,7 @@ module CNVegNitrogenStateType use clm_varctl , only : use_nitrif_denitrif, use_vertsoilc, use_century_decomp !use clm_varctl , only : iulog !rrXbo 10Sep2025 use clm_varpar , only : NUM_ZON, NUM_VEG, VAR_COL, VAR_PFT, & - numpft, CN_zone_weight + numpft, CN_zone_weight, FVEG_MIN use clm_varpar , only : ndecomp_cascade_transitions, ndecomp_pools, nlevcan use clm_varpar , only : nlevdecomp_full, nlevdecomp use clm_varcon , only : spval, ispval, dzsoi_decomp, zisoi @@ -14,7 +14,7 @@ module CNVegNitrogenStateType use decompMod , only : bounds_type use pftconMod , only : npcropmin use PatchType , only : patch - + ! !PUBLIC TYPES: implicit none save @@ -474,7 +474,7 @@ subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft) do p = 0,numpft ! PFT index loop np = np + 1 do nv = 1,NUM_VEG ! defined veg loop - if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>1.e-4) then + if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>FVEG_MIN) then this%deadcrootn_patch (np) = safe( real(cnpft(nc,nz,nv, 48), r8) ) this%deadcrootn_storage_patch (np) = safe( real(cnpft(nc,nz,nv, 49), r8) ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegStateType.F90 index 52e5144af..d46915f2b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNVegStateType.F90 @@ -2,14 +2,13 @@ module CNVegStateType use shr_kind_mod , only : r8 => shr_kind_r8 use nanMod , only : nan - use clm_varpar , only : nlevsno, nlevgrnd, nlevlak, nlevsoi, & - num_zon, num_veg, var_col, var_pft, numpft + use clm_varpar , only : nlevsno, nlevgrnd, nlevlak, nlevsoi, & + num_zon, num_veg, var_col, var_pft, numpft, FVEG_MIN use clm_varcon , only : spval, ispval !use clm_varctl , only : iulog !rrXbo 10Sep2025 use decompMod , only : bounds_type use AnnualFluxDribbler, only : annual_flux_dribbler_type, annual_flux_dribbler_patch - - + ! !PUBLIC TYPES: implicit none save @@ -224,7 +223,7 @@ subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft) do p = 0,numpft ! PFT index loop np = np + 1 do nv = 1,num_veg ! defined veg loop - if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>1.e-4) then + if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>FVEG_MIN) then this%annavg_t2m_patch (np) = cnpft(nc,nz,nv, 24) this%annmax_retransn_patch (np) = cnpft(nc,nz,nv, 25) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 index b075186bd..2cd668d3e 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 @@ -4,13 +4,14 @@ module CanopyStateType use shr_kind_mod , only : r8 => shr_kind_r8 use clm_varpar , only : nlevcan, nvegwcs, numpft, num_zon, num_veg, & - var_col, var_pft + var_col, var_pft, FVEG_MIN use clm_varcon , only : spval !use clm_varctl , only : iulog !rrXbo 10Sep2025 use nanMod , only : nan use decompMod , only : bounds_type - use MAPL_ExceptionHandling + use MAPL_ExceptionHandling + ! !PUBLIC TYPES: implicit none save @@ -155,7 +156,7 @@ subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft, cn5_cold_start, rc) do p = 0,numpft ! PFT index loop np = np + 1 do nv = 1,num_veg ! defined veg loop - if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>1.e-4) then + if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>FVEG_MIN) then ! "old" variables: CNCLM45 and before this%elai_patch (np) = cnpft(nc,nz,nv, 69) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PatchType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PatchType.F90 index 7781c2e59..2d290fa21 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PatchType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PatchType.F90 @@ -5,8 +5,8 @@ module PatchType use decompMod , only : bounds_type use clm_varcon , only : ispval use clm_varctl , only : use_fates - use clm_varpar , only : numpft, NUM_ZON, NUM_VEG, CN_zone_weight - + use clm_varpar , only : numpft, NUM_ZON, NUM_VEG, CN_zone_weight, FVEG_MIN + !----------------------------------------------------------------------- ! !DESCRIPTION: ! Patch data type allocation @@ -138,7 +138,7 @@ subroutine Init(this, bounds, nch, ityp, fveg) this%landunit(np) = nc this%wtlunit(np) = 0. do nv = 1,num_veg ! defined veg loop - if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>1.e-4) then + if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>FVEG_MIN) then this%active(np) = .true. this%wtcol(np) = this%wtcol(np) + fveg(nc,nv,nz) this%wtgcell(np) = this%wtgcell(np) + (fveg(nc,nv,nz)*CN_zone_weight(nz)) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemStateType.F90 index 20fc62c3a..8b67f5a8c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemStateType.F90 @@ -7,12 +7,13 @@ module SoilBiogeochemStateType use clm_varpar , only : ndecomp_cascade_transitions, ndecomp_pools, nlevcan, & nlevsno, nlevgrnd, nlevlak use clm_varpar , only : nlevdecomp_full, nlevdecomp, nlevsoi, & - VAR_COL, VAR_PFT, num_zon, num_veg, numpft + VAR_COL, VAR_PFT, num_zon, num_veg, numpft, FVEG_MIN use clm_varctl , only : use_cn use clm_varcon , only : spval use decompMod , only : bounds_type - use MAPL_ExceptionHandling + use MAPL_ExceptionHandling + ! !PUBLIC TYPES: implicit none save @@ -109,7 +110,7 @@ subroutine Init(this, bounds, nch, cncol, cnpft, ityp, fveg, rc) this%plant_ndemand_col(n) = 0._r8 do p = 0,numpft ! PFT index loop do nv = 1,num_veg ! defined veg loop - if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>1.e-4) then + if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>FVEG_MIN) then this%plant_ndemand_col(n) = this%plant_ndemand_col(n) + cnpft(nc,nz,nv, 75) end if end do ! nv diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SurfaceAlbedoType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SurfaceAlbedoType.F90 index 1228e7cf4..559278dfe 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SurfaceAlbedoType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SurfaceAlbedoType.F90 @@ -3,7 +3,7 @@ module SurfaceAlbedoType use shr_kind_mod , only : r8 => shr_kind_r8 use nanMod , only : nan use clm_varpar , only : numrad, nlevcan, nlevsno, numpft, num_zon, num_veg, & - var_col, var_pft + var_col, var_pft, FVEG_MIN use clm_varcon , only : spval, ispval use decompMod , only : bounds_type @@ -94,7 +94,6 @@ subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft) integer :: begp, endp integer :: begc, endc integer :: np, nc, nz, p, nv, n - real, parameter :: fmin = 1.0e-4 ! Ignore vegetation fractions at or below this value !------------------------------- begp = bounds%begp ; endp = bounds%endp @@ -156,7 +155,7 @@ subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft) this%nrad_patch(np) = 1 do nv = 1,num_veg ! defined veg loop - if (ityp(nc,nv,nz) == p .and. fveg(nc,nv,nz) > fmin) then + if (ityp(nc,nv,nz) == p .and. fveg(nc,nv,nz) > FVEG_MIN) then ! without this if loop we are overwriting each patch multiple times because we loop nv=1:num_veg without ! checking which nv corresponds to the active PFT for patch p. do n = 1, nlevcan diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/clm_varpar.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/clm_varpar.F90 index 21ecad8bc..b566db6e2 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/clm_varpar.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/clm_varpar.F90 @@ -10,12 +10,13 @@ module clm_varpar ! ! !USES: ! - use clm_varpar_shared, only : & - VAR_COL => VAR_COL_51, & - VAR_PFT => VAR_PFT_51, & - numpft => NUM_PFT_CN_51, & - NUM_ZON => NUM_ZON_CN, & - NUM_VEG => NUM_VEG_CN_51 + use clm_varpar_shared, only : & + VAR_COL => VAR_COL_51, & + VAR_PFT => VAR_PFT_51, & + numpft => NUM_PFT_CN_51, & + NUM_ZON => NUM_ZON_CN, & + NUM_VEG => NUM_VEG_CN_51, & + FVEG_MIN ! !PUBLIC TYPES: implicit none diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/filterMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/filterMod.F90 index ad15933f5..c30022349 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/filterMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/filterMod.F90 @@ -3,9 +3,9 @@ module filterMod use shr_kind_mod, only: r8 => shr_kind_r8 use nanMod , only : nan use decompMod , only : bounds_type - use clm_varpar , only : NUM_ZON, NUM_VEG, numpft + use clm_varpar , only : NUM_ZON, NUM_VEG, numpft, FVEG_MIN use pftconMod , only : npcropmin - + ! !PUBLIC TYPES: implicit none save @@ -224,28 +224,27 @@ subroutine init_filter_type(bounds, nch, ityp, fveg, this_filter) do nv = 1,num_veg ! defined veg loop if(ityp(nc,nv,nz)==p) then - if (fveg(nc,nv,nz)>1.e-4) then + if (fveg(nc,nv,nz)>FVEG_MIN) then - this_filter(1)%num_nourbanp = this_filter(1)%num_nourbanp + 1 - this_filter(1)%nourbanp(this_filter(1)%num_nourbanp) = np - - this_filter(1)%num_soilp = this_filter(1)%num_soilp + 1 - this_filter(1)%soilp(this_filter(1)%num_soilp) = np - - ! jkolassa: not sure this is needed, since we do not use prognostic crop information - if(ityp(nc,nv,nz) >= npcropmin) then - this_filter(1)%num_pcropp = this_filter(1)%num_pcropp + 1 - this_filter(1)%pcropp(this_filter(1)%num_pcropp) = np - endif - - -! if (fveg(nc,nv,nz)>1.e-4) then + this_filter(1)%num_nourbanp = this_filter(1)%num_nourbanp + 1 + this_filter(1)%nourbanp(this_filter(1)%num_nourbanp) = np + + this_filter(1)%num_soilp = this_filter(1)%num_soilp + 1 + this_filter(1)%soilp(this_filter(1)%num_soilp) = np + + ! jkolassa: not sure this is needed, since we do not use prognostic crop information + if(ityp(nc,nv,nz) >= npcropmin) then + this_filter(1)%num_pcropp = this_filter(1)%num_pcropp + 1 + this_filter(1)%pcropp(this_filter(1)%num_pcropp) = np + endif + ! if (fveg(nc,nv,nz)>FVEG_MIN) then + this_filter(1)%num_exposedvegp = this_filter(1)%num_exposedvegp + 1 this_filter(1)%exposedvegp(this_filter(1)%num_exposedvegp) = np - - elseif (fveg(nc,nv,nz)<=1.e-4) then - + + elseif (fveg(nc,nv,nz)<=FVEG_MIN) then + this_filter(1)%num_noexposedvegp = this_filter(1)%num_noexposedvegp + 1 this_filter(1)%noexposedvegp(this_filter(1)%num_noexposedvegp) = np diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 index fd13ae4df..9643ab116 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 @@ -62,7 +62,8 @@ module GEOS_CatchCNCLM51GridCompMod USE clm_varpar, ONLY : & NUM_ZON, NUM_VEG, VAR_COL, VAR_PFT, & - CN_zone_weight, map_cat, numpft + CN_zone_weight, map_cat, numpft, & + FVEG_MIN USE MAPL use MAPL_ConstantsMod,only: Tzero => MAPL_TICE, pi => MAPL_PI, MAPL_RHOWTR @@ -132,8 +133,6 @@ module GEOS_CatchCNCLM51GridCompMod ! ROOTL import from GEOS_VegdynGridComp was disabled and brought the look up table ! in order to obtain ROOTL for primary and secondary types. - - ! map catchment type into PFT ! NOTE: Unlike CatchCN40/45, CatchCN51 no longer uses split PFTs (depending on moisture stress) ! --------------------------- @@ -7355,7 +7354,7 @@ subroutine Driver ( RC ) do nz = 1,nzone do nv = 1,nveg do n = 1,ntiles - if (fveg(n,nv,nz)>1.e-4) then ! account for fact that parzone is undefined if fveg = 0 + if (fveg(n,nv,nz)>FVEG_MIN) then ! account for fact that parzone is undefined if fveg = 0 para(n) = para(n) + parzone(n,nv,nz)*wtzone(n,nz)*fveg(n,nv,nz) if(associated(BTRANT)) then btrant(n) = btrant(n) + btran(n,nv,nz)*fveg(n,nv,nz)*wtzone(n,nz) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/clm_varpar_shared.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/clm_varpar_shared.F90 index a5dc27f3a..41c149e9c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/clm_varpar_shared.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/clm_varpar_shared.F90 @@ -16,20 +16,22 @@ module clm_varpar_shared ! Define number of levels - integer, parameter, PUBLIC :: NUM_PFT_CN_40 = 19 ! actual # of pfts (without bare) for Catchment-CN4.0 - integer, parameter, PUBLIC :: NUM_PFT_CN_51 = 15 ! actual # of pfts (without bare) for Catchment-CN5.1 - - integer, parameter, PUBLIC :: NUM_ZON_CN = 3 ! number of CN hydrology zones per tile - - integer, parameter, PUBLIC :: NUM_VEG_CN_40 = 4 ! number of CN PFTs per zone for Catchment-CN4.0 - integer, parameter, PUBLIC :: NUM_VEG_CN_51 = 2 ! number of CN PFTs per zone for Catchment-CN5.1 - - integer, parameter, PUBLIC :: VAR_COL_40 = 40 ! number of CN column restart variables - integer, parameter, PUBLIC :: VAR_PFT_40 = 74 ! number of CN PFT variables per column - - integer, parameter, PUBLIC :: VAR_COL_51 = 37 ! number of CN column restart variables - integer, parameter, PUBLIC :: VAR_PFT_51 = 83 ! number of CN PFT restart variables - + integer, parameter, PUBLIC :: NUM_PFT_CN_40 = 19 ! actual # of pfts (without bare) for Catchment-CN4.0 + integer, parameter, PUBLIC :: NUM_PFT_CN_51 = 15 ! actual # of pfts (without bare) for Catchment-CN5.1 + + integer, parameter, PUBLIC :: NUM_ZON_CN = 3 ! number of CN hydrology zones per tile + + integer, parameter, PUBLIC :: NUM_VEG_CN_40 = 4 ! number of CN PFTs per zone for Catchment-CN4.0 + integer, parameter, PUBLIC :: NUM_VEG_CN_51 = 2 ! number of CN PFTs per zone for Catchment-CN5.1 + + integer, parameter, PUBLIC :: VAR_COL_40 = 40 ! number of CN column restart variables + integer, parameter, PUBLIC :: VAR_PFT_40 = 74 ! number of CN PFT variables per column + + integer, parameter, PUBLIC :: VAR_COL_51 = 37 ! number of CN column restart variables + integer, parameter, PUBLIC :: VAR_PFT_51 = 83 ! number of CN PFT restart variables + + real, parameter, PUBLIC :: FVEG_MIN = 1.e-4 ! ignore vegetation when vegetation fraction is at or below this value + end module clm_varpar_shared ! ============================ EOF =========================================================== From 7454ddac18b3adea084629c03c1fe4745dc369eb Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 11 Sep 2025 14:57:57 -0400 Subject: [PATCH 09/29] fixed indent, improved vertical alignment, clarified/edited comments (SurfaceAlbedoType.F90) --- .../CLM51/SurfaceAlbedoType.F90 | 88 ++++++++++--------- 1 file changed, 48 insertions(+), 40 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SurfaceAlbedoType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SurfaceAlbedoType.F90 index 559278dfe..a6e59b89b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SurfaceAlbedoType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SurfaceAlbedoType.F90 @@ -10,10 +10,10 @@ module SurfaceAlbedoType ! !PUBLIC TYPES: implicit none save -! -! !PUBLIC MEMBER FUNCTIONS: + ! + ! !PUBLIC MEMBER FUNCTIONS: - ! + ! type, public :: surfalb_type real(r8), pointer :: coszen_col (:) ! col cosine of solar zenith angle @@ -61,34 +61,37 @@ module SurfaceAlbedoType integer , pointer :: nrad_patch (:) ! patch number of canopy layers, above snow for radiative transfer real(r8) , pointer :: vcmaxcintsun_patch (:) ! patch leaf to canopy scaling coefficient, sunlit leaf vcmax real(r8) , pointer :: vcmaxcintsha_patch (:) ! patch leaf to canopy scaling coefficient, shaded leaf vcmax - + contains - + procedure, public :: Init - -end type surfalb_type -type(surfalb_type), public, target, save :: surfalb_inst + + end type surfalb_type + + type(surfalb_type), public, target, save :: surfalb_inst contains -!--------------------------------------------------- + !--------------------------------------------------- + subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft) - ! !DESCRIPTION: -! Initialize CTSM surface albedo needed for calling CTSM routines -! jk Apr 2021: type is allocated and initialized to NaN; values are assigned from Catchment states before calls to CLM subroutines are made -! this type is only used to be able to pass Catchment states and fluxes to CLM subroutines in the format they expect -! -! !ARGUMENTS: + ! !DESCRIPTION: + ! Initialize CTSM surface albedo needed for calling CTSM routines + ! jk Apr 2021: type is allocated and initialized to NaN; values are assigned from Catchment states before calls to CLM subroutines are made + ! this type is only used to be able to pass Catchment states and fluxes to CLM subroutines in the format they expect + ! + ! !ARGUMENTS: implicit none + !INPUT/OUTPUT - type(bounds_type), intent(in) :: bounds - integer, intent(in) :: nch ! number of Catchment tiles - real, dimension(nch,num_zon,var_col), intent(in) :: cncol ! column-level restart variable array - real, dimension(nch,num_zon,num_veg,var_pft), intent(in) :: cnpft ! pft-level (patch-level) restart variable array - integer, dimension(nch,num_veg,num_zon), intent(in) :: ityp - real, dimension(nch,num_veg,num_zon), intent(in) :: fveg - class(surfalb_type) :: this + type(bounds_type), intent(in) :: bounds + integer, intent(in) :: nch ! number of Catchment tiles + real, dimension(nch,num_zon,var_col), intent(in) :: cncol ! column-level restart variable array + real, dimension(nch,num_zon,num_veg,var_pft), intent(in) :: cnpft ! pft-level (patch-level) restart variable array + integer, dimension(nch,num_veg,num_zon), intent(in) :: ityp + real, dimension(nch,num_veg,num_zon), intent(in) :: fveg + class(surfalb_type) :: this ! LOCAL integer :: begp, endp @@ -147,30 +150,35 @@ subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft) ! initialize variables from restart files np = 0 - do nc = 1,nch ! catchment tile loop + + do nc = 1,nch ! catchment tile loop do nz = 1,num_zon ! CN zone loop - do p = 0,numpft ! PFT index loop - np = np + 1 + do p = 0,numpft ! PFT index loop + np = np + 1 + this%nrad_patch(np) = 1 do nv = 1,num_veg ! defined veg loop - if (ityp(nc,nv,nz) == p .and. fveg(nc,nv,nz) > FVEG_MIN) then - ! without this if loop we are overwriting each patch multiple times because we loop nv=1:num_veg without - ! checking which nv corresponds to the active PFT for patch p. - do n = 1, nlevcan - this%tlai_z_patch(np,n) = cnpft(nc,nz,nv,73); if (isnan(this%tlai_z_patch(np,n))) this%tlai_z_patch(np,n)=0._r8 - this%tsai_z_patch(np,n) = cnpft(nc,nz,nv,74); if (isnan(this%tsai_z_patch(np,n))) this%tsai_z_patch(np,n)=0._r8 - end do - this%vcmaxcintsun_patch(np) = 1._r8 - this%vcmaxcintsha_patch(np) = 1._r8 - exit ! <- stop after filling the matching nv for this patch - end if ! ityp =p - end do !nv - end do ! p - end do ! nz - end do ! nc + if (ityp(nc,nv,nz) == p .and. fveg(nc,nv,nz) > FVEG_MIN) then ! "if" condition was missing in original implementation !rrXbo 10Sep2025 + + do n = 1, nlevcan + this%tlai_z_patch(np,n) = cnpft(nc,nz,nv,73); if (isnan(this%tlai_z_patch(np,n))) this%tlai_z_patch(np,n)=0._r8 + this%tsai_z_patch(np,n) = cnpft(nc,nz,nv,74); if (isnan(this%tsai_z_patch(np,n))) this%tsai_z_patch(np,n)=0._r8 + end do + + this%vcmaxcintsun_patch(np) = 1._r8 + this%vcmaxcintsha_patch(np) = 1._r8 + + exit ! stop after filling the matching nv for this patch + + end if ! ityp =p + end do ! nv + end do ! p + end do ! nz + end do ! nc + end subroutine Init end module SurfaceAlbedoType From 7f5950400dbdc16bdc5e49f9e32a7a8beb2abbf5 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 11 Sep 2025 15:05:36 -0400 Subject: [PATCH 10/29] fixed indentation, improved vertical alignment (CNCLM_init_mod.F90) --- .../CLM51/CNCLM_init_mod.F90 | 240 +++++++++--------- 1 file changed, 121 insertions(+), 119 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_init_mod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_init_mod.F90 index 1bc97bb76..ceabfe4a4 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_init_mod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_init_mod.F90 @@ -36,14 +36,14 @@ module CNCLM_initMod use PatchType use ColumnType use ch4Mod - use SoilBiogeochemDecompCascadeConType, only : init_decomp_cascade_constants + use SoilBiogeochemDecompCascadeConType , only : init_decomp_cascade_constants use ActiveLayerMod use CropType use CNDVType - use LandunitType , only : lun + use LandunitType , only : lun use RootBiophysMod - use CNMRespMod , only : readCNMRespParams => readParams, CNMRespReadNML - use CNSharedParamsMod , only : CNParamsReadShared + use CNMRespMod , only : readCNMRespParams => readParams, CNMRespReadNML + use CNSharedParamsMod , only : CNParamsReadShared use spmdMod use Wateratm2lndBulkType use WaterDiagnosticBulkType @@ -57,7 +57,7 @@ module CNCLM_initMod use CNVegetationFacade use initSubgridMod use CN2CLMType - use WaterType , only : water_type + use WaterType , only : water_type use CNBalanceCheckMod use SoilBiogeochemDecompCascadeBGCMod , only : init_decompcascade_bgc, DecompCascadeBGCreadNML @@ -86,60 +86,61 @@ module CNCLM_initMod use SoilBiogeochemNitrifDenitrifMod , only : readSoilBiogeochemNitrifDenitrifParams => readParams use SoilStateInitTimeConstMod , only : readParams_SoilStateInitTimeConst => readParams, & SoilStateInitTimeConst - use clm_varpar , only : numpft, num_zon, num_veg, var_pft, var_col, & - nlevgrnd, nlevsoi + use clm_varpar , only : numpft, num_zon, num_veg, var_pft, var_col, & + nlevgrnd, nlevsoi use MAPL ! , only : NetCDF4_FileFormatter, pFIO_READ use ESMF - implicit none - private + implicit none + private - class(nutrient_competition_method_type), public, allocatable :: nutrient_competition_method - class(fire_method_type), allocatable :: cnfire_method - type(saturated_excess_runoff_type), public :: saturated_excess_runoff_inst - type(waterflux_type), public :: waterflux_inst - type(waterfluxbulk_type), public :: waterfluxbulk_inst + class(nutrient_competition_method_type), public, allocatable :: nutrient_competition_method + class(fire_method_type), allocatable :: cnfire_method + type(saturated_excess_runoff_type), public :: saturated_excess_runoff_inst + type(waterflux_type), public :: waterflux_inst + type(waterfluxbulk_type), public :: waterfluxbulk_inst -! !PUBLIC MEMBER FUNCTIONS: - public :: CN_init + ! !PUBLIC MEMBER FUNCTIONS: + public :: CN_init - contains +contains -!------------------------------------------------------ - subroutine CN_init(NLFilename, nch,ityp,fveg,cncol,cnpft,lats,lons,dtcn,paramfile,water_inst,bgc_vegetation_inst,cn5_cold_start) + !------------------------------------------------------ - !ARGUMENTS - implicit none - !INPUT/OUTPUT - character(256), intent(in) :: NLFilename ! len=256 for consistency with SHR_KIND_CL of CLM51/shr_kind_mod.F90 - integer, intent(in) :: nch ! number of tiles - integer, dimension(nch,num_veg,num_zon), intent(in) :: ityp ! PFT index - real, dimension(nch,num_veg,num_zon), intent(in) :: fveg ! PFT fraction - real, dimension(nch,num_zon,var_col), intent(in) :: cncol ! column-level CN restart variables - real, dimension(nch,num_zon,num_veg,var_pft), intent(in) :: cnpft ! patch/pft-level CN restart variables - real, dimension(nch), intent(in) :: lats ! Catchment tile latitudes [rad] - real, dimension(nch), intent(in) :: lons ! Catchment tile longitudes [rad] - real, intent(in) :: dtcn ! Catchment-CN step size - character(len=ESMF_MAXSTR), intent(in) :: paramfile - logical, optional, intent(in) :: cn5_cold_start ! cold start for the CLM variables that are new in Catchment-CN5.0 - type(water_type), intent(out) :: water_inst - type(cn_vegetation_type), intent(out) :: bgc_vegetation_inst - !LOCAL - - type(Netcdf4_fileformatter) :: ncid - integer :: rc, status, ndt + subroutine CN_init(NLFilename, nch,ityp,fveg,cncol,cnpft,lats,lons,dtcn,paramfile,water_inst,bgc_vegetation_inst,cn5_cold_start) - !----------------------------------------- + !ARGUMENTS + implicit none + !INPUT/OUTPUT + character(256), intent(in) :: NLFilename ! len=256 for consistency with SHR_KIND_CL of CLM51/shr_kind_mod.F90 + integer, intent(in) :: nch ! number of tiles + integer, dimension(nch,num_veg,num_zon), intent(in) :: ityp ! PFT index + real, dimension(nch,num_veg,num_zon), intent(in) :: fveg ! PFT fraction + real, dimension(nch,num_zon,var_col), intent(in) :: cncol ! column-level CN restart variables + real, dimension(nch,num_zon,num_veg,var_pft), intent(in) :: cnpft ! patch/pft-level CN restart variables + real, dimension(nch), intent(in) :: lats ! Catchment tile latitudes [rad] + real, dimension(nch), intent(in) :: lons ! Catchment tile longitudes [rad] + real, intent(in) :: dtcn ! Catchment-CN step size + character(len=ESMF_MAXSTR), intent(in) :: paramfile + logical, optional, intent(in) :: cn5_cold_start ! cold start for the CLM variables that are new in Catchment-CN5.0 + type(water_type), intent(out) :: water_inst + type(cn_vegetation_type), intent(out) :: bgc_vegetation_inst + !LOCAL -! initialize CN step size + type(Netcdf4_fileformatter) :: ncid + integer :: rc, status, ndt - ndt = get_step_size( nint(dtcn) ) + !----------------------------------------- -! initialize CN model -! ------------------- + ! initialize CN step size + + ndt = get_step_size( nint(dtcn) ) + + ! initialize CN model + ! ------------------- call spmd_init() @@ -149,30 +150,30 @@ subroutine CN_init(NLFilename, nch,ityp,fveg,cncol,cnpft,lats,lons,dtcn,paramfil call init_clm_varctl() - call bounds%Init (nch) + call bounds%Init (nch) ! initialize subrgid types - call patch%Init (bounds, nch, ityp, fveg) + call patch%Init (bounds, nch, ityp, fveg) - call col%Init (bounds, nch) + call col%Init (bounds, nch) - call lun%Init (bounds, nch) + call lun%Init (bounds, nch) - call grc%Init (bounds, nch, cnpft, lats, lons) + call grc%Init (bounds, nch, cnpft, lats, lons) ! create subgrid structure - call clm_ptrs_compdown (bounds) + call clm_ptrs_compdown (bounds) ! initialize filters - call allocFilters (bounds, nch, ityp, fveg) + call allocFilters (bounds, nch, ityp, fveg) ! read parameters and configurations from namelist file - + call CNPhenologyReadNML ( NLFilename ) - call dynSubgridControl_init ( ) + call dynSubgridControl_init ( ) call CNFireReadNML ( NLFilename ) call CNNDynamicsReadNML ( NLFilename ) call photosyns_inst%ReadNML ( NLFilename ) @@ -183,57 +184,57 @@ subroutine CN_init(NLFilename, nch,ityp,fveg,cncol,cnpft,lats,lons,dtcn,paramfil ! initialize states and fluxes - call pftcon%init_pftcon_type (paramfile) + call pftcon%init_pftcon_type (paramfile) - call bgc_vegetation_inst%Init(bounds, NLFilename, nch, ityp, fveg, cncol, cnpft, paramfile, cn5_cold_start) + call bgc_vegetation_inst%Init (bounds, NLFilename, nch, ityp, fveg, cncol, cnpft, paramfile, cn5_cold_start) - call atm2lnd_inst%Init (bounds) + call atm2lnd_inst%Init (bounds) - call temperature_inst%Init (bounds) + call temperature_inst%Init (bounds) - call soilstate_inst%Init (bounds) + call soilstate_inst%Init (bounds) - call SoilStateInitTimeConst (bounds, soilstate_inst, NLFilename) ! sets hydraulic and thermal soil properties + call SoilStateInitTimeConst (bounds, soilstate_inst, NLFilename) ! sets hydraulic and thermal soil properties - call water_inst%Init (bounds) + call water_inst%Init (bounds) - call canopystate_inst%Init (bounds, nch, ityp, fveg, cncol, cnpft, cn5_cold_start) + call canopystate_inst%Init (bounds, nch, ityp, fveg, cncol, cnpft, cn5_cold_start) - call solarabs_inst%Init (bounds) + call solarabs_inst%Init (bounds) - call surfalb_inst%Init (bounds, nch, ityp, fveg, cncol, cnpft) + call surfalb_inst%Init (bounds, nch, ityp, fveg, cncol, cnpft) - call ozone_inst%Init (bounds) + call ozone_inst%Init (bounds) - call photosyns_inst%Init (bounds, nch, ityp, fveg, cncol, cnpft, cn5_cold_start) + call photosyns_inst%Init (bounds, nch, ityp, fveg, cncol, cnpft, cn5_cold_start) - call soilbiogeochem_carbonstate_inst%Init(bounds, nch, cncol) + call soilbiogeochem_carbonstate_inst%Init (bounds, nch, cncol) call soilbiogeochem_nitrogenstate_inst%Init(bounds, nch, cncol) - call soilbiogeochem_state_inst%Init (bounds, nch, cncol, cnpft, ityp, fveg) + call soilbiogeochem_state_inst%Init (bounds, nch, cncol, cnpft, ityp, fveg) - call soilbiogeochem_carbonflux_inst%Init (bounds) + call soilbiogeochem_carbonflux_inst%Init (bounds) - call soilbiogeochem_nitrogenflux_inst%Init(bounds) + call soilbiogeochem_nitrogenflux_inst%Init (bounds) - call ch4_inst%Init (bounds) + call ch4_inst%Init (bounds) - call init_decomp_cascade_constants (use_century_decomp) - - call active_layer_inst%Init (bounds) + call init_decomp_cascade_constants (use_century_decomp) - call crop_inst%Init (bounds) + call active_layer_inst%Init (bounds) - call dgvs_inst%Init (bounds) + call crop_inst%Init (bounds) - call saturated_excess_runoff_inst%Init(bounds) + call dgvs_inst%Init (bounds) - call energyflux_inst%Init (bounds) + call saturated_excess_runoff_inst%Init (bounds) - call frictionvel_inst%Init (bounds) + call energyflux_inst%Init (bounds) - call cn_balance_inst%Init (bounds) + call frictionvel_inst%Init (bounds) + + call cn_balance_inst%Init (bounds) ! initialize rooting profile parameters from namelist @@ -242,67 +243,68 @@ subroutine CN_init(NLFilename, nch,ityp,fveg,cncol,cnpft,lats,lons,dtcn,paramfil ! initialize root fractions call init_vegrootfr(bounds, nlevsoi, nlevgrnd, & - soilstate_inst%rootfr_patch(bounds%begp:bounds%endp,1:nlevgrnd),'water') + soilstate_inst%rootfr_patch(bounds%begp:bounds%endp,1:nlevgrnd),'water') call init_vegrootfr(bounds, nlevsoi, nlevgrnd, & - soilstate_inst%crootfr_patch(bounds%begp:bounds%endp,1:nlevgrnd),'carbon') + soilstate_inst%crootfr_patch(bounds%begp:bounds%endp,1:nlevgrnd),'carbon') - ! allocate CLM arrays that are not allocated in their modules + ! allocate CLM arrays that are not allocated in their modules allocate(nutrient_competition_method, & source=create_nutrient_competition_method(bounds)) ! jkolassa: this allocates and initializes the nutrient_competition_method_type - ! initialize CLM parameters from parameter file + ! initialize CLM parameters from parameter file - call ncid%open(trim(paramfile),pFIO_READ, RC=status) + call ncid%open(trim(paramfile),pFIO_READ, RC=status) - call readCNMRespParams(ncid) - call CNParamsReadShared(ncid, NLFilename) ! this is called CN params but really is for the soil biogeochem parameters - call readSoilBiogeochemDecompCnParams(ncid) - call readSoilBiogeochemDecompBgcParams(ncid) - call nutrient_competition_method%readParams(ncid) - call readSoilBiogeochemDecompParams(ncid) - call readCNPhenolParams(ncid) - call readSoilBiogeochemLittVertTranspParams(ncid) - call photosyns_inst%ReadParams( ncid ) - call readSoilBiogeochemNLeachingParams(ncid) - call readSoilBiogeochemCompetitionParams(ncid) - call readSoilBiogeochemPotentialParams(ncid) - call readCNGapMortalityParams(ncid) - call readCNFUNParams(ncid) - call readSoilBiogeochemNitrifDenitrifParams(ncid) - !call readParams_SoilStateInitTimeConst(ncid) + call readCNMRespParams( ncid) + call CNParamsReadShared( ncid, NLFilename) ! this is called CN params but really is for the soil biogeochem parameters + call readSoilBiogeochemDecompCnParams( ncid) + call readSoilBiogeochemDecompBgcParams( ncid) + call nutrient_competition_method%readParams(ncid) + call readSoilBiogeochemDecompParams( ncid) + call readCNPhenolParams( ncid) + call readSoilBiogeochemLittVertTranspParams(ncid) + call photosyns_inst%ReadParams( ncid) + call readSoilBiogeochemNLeachingParams( ncid) + call readSoilBiogeochemCompetitionParams( ncid) + call readSoilBiogeochemPotentialParams( ncid) + call readCNGapMortalityParams( ncid) + call readCNFUNParams( ncid) + call readSoilBiogeochemNitrifDenitrifParams(ncid) + !call readParams_SoilStateInitTimeConst( ncid) - call ncid%close(rc=status) + call ncid%close(rc=status) - ! initialize types that depend on parameters + ! initialize types that depend on parameters - call CNPhenologyInit (bounds) - call SoilBiogeochemCompetitionInit (bounds) - call CNFUNInit(bounds,bgc_vegetation_inst%cnveg_state_inst,bgc_vegetation_inst%cnveg_carbonstate_inst,bgc_vegetation_inst%cnveg_nitrogenstate_inst) + call CNPhenologyInit (bounds) + call SoilBiogeochemCompetitionInit (bounds) + call CNFUNInit(bounds,bgc_vegetation_inst%cnveg_state_inst,bgc_vegetation_inst%cnveg_carbonstate_inst,bgc_vegetation_inst%cnveg_nitrogenstate_inst) - ! Initialize precision control for soil biogeochemistry (use soilbiogeochem_carbonstate three times, since we do not currently use isotopes) + ! Initialize precision control for soil biogeochemistry (use soilbiogeochem_carbonstate three times, since we do not currently use isotopes) call SoilBiogeochemPrecisionControlInit( soilbiogeochem_carbonstate_inst, soilbiogeochem_carbonstate_inst, & - soilbiogeochem_carbonstate_inst, soilbiogeochem_nitrogenstate_inst) + soilbiogeochem_carbonstate_inst, soilbiogeochem_nitrogenstate_inst) + + ! call FireMethodInit(bounds,paramfile) - ! call FireMethodInit(bounds,paramfile) + if (use_century_decomp) then + call init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, & + soilstate_inst ) + else + call init_decompcascade_cn(bounds, soilbiogeochem_state_inst) + end if - if (use_century_decomp) then - call init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, & - soilstate_inst ) - else - call init_decompcascade_cn(bounds, soilbiogeochem_state_inst) - end if + ! initialize custom type used to pass Catchment information to nested CLM fire types - ! initialize custom type used to pass Catchment information to nested CLM fire types + call cn2clm_inst%Init (bounds) - call cn2clm_inst%Init (bounds) + ! initialize radiation time - ! initialize radiation time + call update_rad_dtime(.true.) - call update_rad_dtime(.true.) + end subroutine CN_init - end subroutine CN_init end module CNCLM_initMod From 58e67706448b9503aa0c4be39743d4db9ddbbe01 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 11 Sep 2025 15:08:32 -0400 Subject: [PATCH 11/29] bug fix: corrected assignment of restart variables cncol(:,:,[8,22]) (CNCLM_DriverMod.F90) --- .../GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 index fac8b19ca..1cae78e58 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 @@ -619,7 +619,7 @@ subroutine CN_exit(nch,ityp,fveg,cncol,cnpft) ! jkolassa: variables below transitioned from being column-level to being gridcell-level in CLM; ! assuming here that quantities are spread over zones according to zone weight cncol(nc,nz, 7) = bgc_vegetation_inst%c_products_inst%prod100_grc( nc) * CN_zone_weight(nz) - cncol(nc,nz, 8) = bgc_vegetation_inst%c_products_inst%prod100_grc( nc) * CN_zone_weight(nz) + cncol(nc,nz, 8) = bgc_vegetation_inst%c_products_inst%prod10_grc( nc) * CN_zone_weight(nz) cncol(nc,nz, 9) = bgc_vegetation_inst%cnveg_carbonstate_inst%seedc_grc( nc) * CN_zone_weight(nz) cncol(nc,nz,14) = bgc_vegetation_inst%cnveg_carbonstate_inst%totc_col( n ) @@ -628,7 +628,7 @@ subroutine CN_exit(nch,ityp,fveg,cncol,cnpft) cncol(nc,nz,21) = bgc_vegetation_inst%n_products_inst%prod100_grc( nc) * CN_zone_weight(nz) - cncol(nc,nz,22) = bgc_vegetation_inst%n_products_inst%prod100_grc( nc) * CN_zone_weight(nz) + cncol(nc,nz,22) = bgc_vegetation_inst%n_products_inst%prod10_grc( nc) * CN_zone_weight(nz) cncol(nc,nz,23) = bgc_vegetation_inst%cnveg_nitrogenstate_inst%seedn_grc( nc) * CN_zone_weight(nz) cncol(nc,nz,24) = soilbiogeochem_nitrogenstate_inst%sminn_vr_col( n,1) From 88875762270297a2e25f205a6bb930327b755960 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 11 Sep 2025 15:30:02 -0400 Subject: [PATCH 12/29] fixed "if" condition involving "optional" argument (CanopyStateType.F90) --- .../GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 index 2cd668d3e..1cf4817ff 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CanopyStateType.F90 @@ -103,8 +103,10 @@ subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft, cn5_cold_start, rc) begg = bounds%begg ; endg = bounds%endg ! check whether a cn5_cold_start option was set and change cold_start accordingly - if (present(cn5_cold_start) .and. (cn5_cold_start.eqv..true.)) then - cold_start = .true. + if (present(cn5_cold_start)) then + if (cn5_cold_start .eqv. .true.) then + cold_start = .true. + end if end if ! jkolassa: if cold_start is false, check that both CNCOL and CNPFT have the expected size for CNCLM50, else abort From 766560f9a106526542cb9cd68da08bee5d623749 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 11 Sep 2025 15:57:23 -0400 Subject: [PATCH 13/29] reverted init of "alphapsn" back to "spval", added comments (PhotosynthesisMod.F90) --- .../CLM51/PhotosynthesisMod.F90 | 57 +++++++++++-------- 1 file changed, 33 insertions(+), 24 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PhotosynthesisMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PhotosynthesisMod.F90 index e7e7c4bc3..f06f780e9 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PhotosynthesisMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PhotosynthesisMod.F90 @@ -264,13 +264,12 @@ subroutine Init(this,bounds,nch,ityp,fveg,cncol,cnpft,cn5_cold_start,rc) cold_start = .true. end if - ! jkolassa: if cold_start is false, check that both CNCOL and CNPFT have the expected size for CNCLM50, else abort + ! jkolassa: if cold_start is false, check that both CNCOL and CNPFT have the expected size for CNCLM51, else abort if ((cold_start.eqv..false.) .and. ((size(cncol,3).ne.var_col) .or. & (size(cnpft,4).ne.var_pft))) then _ASSERT(.FALSE.,'option CNCLM50_cold_start = .FALSE. requires a CNCLM50 restart file') end if - allocate(this%c3flag_patch (begp:endp)) ; this%c3flag_patch (:) =.false. allocate(this%ac_phs_patch (begp:endp,2,1:nlevcan)) ; this%ac_phs_patch (:,:,:) = nan allocate(this%aj_phs_patch (begp:endp,2,1:nlevcan)) ; this%aj_phs_patch (:,:,:) = nan @@ -346,9 +345,10 @@ subroutine Init(this,bounds,nch,ityp,fveg,cncol,cnpft,cn5_cold_start,rc) allocate(this%luvcmax25top_patch(begp:endp)) ; this%luvcmax25top_patch(:) = nan allocate(this%lujmax25top_patch (begp:endp)) ; this%lujmax25top_patch(:) = nan allocate(this%lutpu25top_patch (begp:endp)) ; this%lutpu25top_patch(:) = nan -!! -! allocate(this%psncanopy_patch (begp:endp)) ; this%psncanopy_patch (:) = nan -! allocate(this%lmrcanopy_patch (begp:endp)) ; this%lmrcanopy_patch (:) = nan + !! + ! allocate(this%psncanopy_patch (begp:endp)) ; this%psncanopy_patch (:) = nan + ! allocate(this%lmrcanopy_patch (begp:endp)) ; this%lmrcanopy_patch (:) = nan + if(use_luna)then ! NOTE(bja, 2015-09) because these variables are only allocated ! when luna is turned on, they can not be placed into associate @@ -372,21 +372,30 @@ subroutine Init(this,bounds,nch,ityp,fveg,cncol,cnpft,cn5_cold_start,rc) ! ! this%modifyphoto_and_lmr_forcrop = .true. ! jkolassa, Feb 2022: Default for CLM50 and up - ! initialize types from restart file or through cold start values - - np = 0 - do nc = 1,nch ! catchment tile loop - do nz = 1,num_zon ! CN zone loop - do p = 0,numpft ! PFT index loop - np = np + 1 - ! Using 0._r8 is safer than spval, because spval value can propagate if any path uses alpha before it’s overwritten. - ! If someday we do want alpha to persist across restarts, we need to: (1) pick free indices, (2) write them in CN_exit, and (3) read those same indices here. - ! Until then, don’t read them from cnpft ... we will use what we calculate. - this%alphapsnsun_patch(np) = 0._r8 - this%alphapsnsha_patch(np) = 0._r8 - end do ! p - end do ! nz - end do ! nc + + ! The original CatchCN51 implementation was filling "alphapsn" from the CatchCN51 restart (CNPFT(:,:,[76,77]), + ! but the variables indices [76,77] used to do so are associated with vegwp_patch(:,[1,2]) in subroutine CN_exit(). + ! It is not clear if "alphapsn" is used anywhere in CatchCN51. For now, change the generic NaN + ! initialization from above to "spval" for the "defined" PFTs. - reichle, 11 Sep 2025 + + np = 0 + do nc = 1,nch ! catchment tile loop + do nz = 1,num_zon ! CN zone loop + do p = 0,numpft ! PFT index loop + np = np + 1 + do nv = 1,num_veg ! defined veg loop + + if(ityp(nc,nv,nz)==p .and. fveg(nc,nv,nz)>FVEG_MIN) then + + this%alphapsnsun_patch(np) = spval ! spval instead of NaN from generic init above !rrXbo, 10Sep2025 + this%alphapsnsha_patch(np) = spval ! spval instead of NaN from generic init above !rrXbo, 10Sep2025 + + end if ! ityp==p .and. fveg>FVEG_MIN + + end do ! nv + end do ! p + end do ! nz + end do ! nc end subroutine Init @@ -750,8 +759,8 @@ subroutine TimeStepInit (this, bounds) this%fpsn_wp_patch(p) = 0._r8 if ( use_c13 ) then - this%alphapsnsun_patch(p) = 0._r8 - this%alphapsnsha_patch(p) = 0._r8 + this%alphapsnsun_patch(p) = 0._r8 ! should this be 1.? see subroutine PhotosynthesisHydraulicStress(), [irrelevant as long as use_c13=.false. remains hard-coded], reichle, 11 Sep 2025] + this%alphapsnsha_patch(p) = 0._r8 ! should this be 1.? see subroutine PhotosynthesisHydraulicStress(), [irrelevant as long as use_c13=.false. remains hard-coded], reichle, 11 Sep 2025] this%c13_psnsun_patch(p) = 0._r8 this%c13_psnsha_patch(p) = 0._r8 endif @@ -1961,8 +1970,8 @@ subroutine PhotosynthesisHydraulicStress ( bounds, fn, filterp, & kp_z(p,sha,iv) = 0._r8 if ( use_c13 ) then - alphapsn_sun(p) = 1._r8 - alphapsn_sha(p) = 1._r8 + alphapsn_sun(p) = 1._r8 ! should this be 0.? see subroutine TimeStepInit(), [irrelevant as long as use_c13=.false. remains hard-coded], reichle, 11 Sep 2025] + alphapsn_sha(p) = 1._r8 ! should this be 0.? see subroutine TimeStepInit(), [irrelevant as long as use_c13=.false. remains hard-coded], reichle, 11 Sep 2025] end if else ! day time From cef80d917409918a07fc29fb566fa44a4e2ef615 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 11 Sep 2025 16:19:20 -0400 Subject: [PATCH 14/29] added comment (CNCLM_DriverMod.F90) --- .../GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 index 1cae78e58..ea1bc9975 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 @@ -307,15 +307,15 @@ subroutine CN_Driver(istep,nch,ityp,fveg,ndep,tp1,tairm,psis,bee,dayl,btran_fire temperature_inst%t_ref2m_patch(p) = tairm(nc) temperature_inst%soila10_patch(p) = tg10d(nc) temperature_inst%t_a5min_patch(p) = t2m5d(nc) - cn2clm_inst%btran2_patch_cn2clm(p) = btran_fire(nc,nz) - + + cn2clm_inst%btran2_patch_cn2clm(p) = btran_fire(nc,nz) ! add root zone soil wetness ("btran_fire") into CN2CLM structure; fire modules read it from there ! cnfire_li2016_inst%cnfire_base_type%btran2_patch(p) = btran_fire(nc,nz) ! cnfire_li2021_inst%cnfire_base_type%btran2_patch(p) = btran_fire(nc,nz) water_inst%wateratm2lndbulk_inst%prec60_patch(p) = prec60d(nc) water_inst%wateratm2lndbulk_inst%prec10_patch(p) = prec10d(nc) water_inst%wateratm2lndbulk_inst%rh30_patch(p) = rh30(nc) - frictionvel_inst%forc_hgt_u_patch(p) = 30. ! following CNCLM45 implementation, but this should be available from the GridComp + frictionvel_inst%forc_hgt_u_patch(p) = 30. ! following CNCLM45 implementation, but this should be available from the GridComp do nv = 1,num_veg ! defined veg loop From c3bb4a3323dbed564b665c3361ed57ab51b37b59 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 11 Sep 2025 16:37:27 -0400 Subject: [PATCH 15/29] enforce non-negative values without also setting upper bound (SoilBiogeochemNitrogenStateType.F90) --- .../CLM51/SoilBiogeochemNitrogenStateType.F90 | 29 ++++++++++--------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 index f31b2ff93..3f7b48c53 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 @@ -105,14 +105,15 @@ pure elemental real(r8) function safe(v) safe = merge(v, 0._r8, .true.) ! valid(v)) !rrXbo 10Sep2025 end function safe - pure elemental real(r8) function clamp_nonneg(v) result(out) - real(r8), intent(in) :: v - real(r8) :: x - x = safe(v) ! NaN/CF-fill/SPVAL -> 0 (via safe), else pass through - if (x < 0._r8) x = 0._r8 ! pools are nonnegative - if (x > N_MAX) x = 0._r8 ! treat absurd values as missing - out = x - end function clamp_nonneg +!rrXbo, 10Sep2025 +! pure elemental real(r8) function clamp_nonneg(v) result(out) +! real(r8), intent(in) :: v +! real(r8) :: x +! x = safe(v) ! NaN/CF-fill/SPVAL -> 0 (via safe), else pass through +! if (x < 0._r8) x = 0._r8 ! pools are nonnegative +! if (x > N_MAX) x = 0._r8 ! treat absurd values as missing +! out = x +! end function clamp_nonneg !------------------------------------------- subroutine Init(this, bounds, nch, cncol) @@ -199,17 +200,17 @@ subroutine Init(this, bounds, nch, cncol) ! 2) total mineral N ! jkolassa May 2022: for now nlevdecomp_full = 1; will need to add loop if we introduce more soil layers - this%sminn_vr_col(n,1:nlevdecomp_full) = clamp_nonneg( real(cncol(nc,nz,24), r8) ) + this%sminn_vr_col(n,1:nlevdecomp_full) = max(cncol(nc,nz,24), 0._r8 ) !rrXbo, 10Sep2025 this%sminn_col(n) = this%sminn_vr_col(n,1) ! 3) split into NO3/NH4 if (no_cn51_rst) then ! jkolassa Nov 2024: when no CN51 restart file is available compute NO3 and NH4 from N ! derive split if no CN51 restart present - this%smin_no3_col(n) = clamp_nonneg( (1.25_r8/2.25_r8)*this%sminn_col(n) ) - this%smin_nh4_col(n) = clamp_nonneg( this%sminn_col(n)/2.25_r8 ) + this%smin_no3_col(n) = max( (1.25_r8/2.25_r8)*this%sminn_col(n), 0._r8 ) !rrXbo, 10Sep2025 + this%smin_nh4_col(n) = max( this%sminn_col(n)/2.25_r8 , 0._r8 ) !rrXbo, 10Sep2025 else - this%smin_no3_col(n) = clamp_nonneg( real(cncol(nc,nz,36), r8) ) - this%smin_nh4_col(n) = clamp_nonneg( real(cncol(nc,nz,37), r8) ) + this%smin_no3_col(n) = max( cncol(nc,nz,36), 0._r8 ) !rrXbo, 10Sep2025 + this%smin_nh4_col(n) = max( cncol(nc,nz,37), 0._r8 ) !rrXbo, 10Sep2025 end if ! 4) mirror scalars into vertically-resolved arrays (nlevdecomp_full is 1 today, but this stays valid if >1 later) @@ -218,7 +219,7 @@ subroutine Init(this, bounds, nch, cncol) do np = 1,ndecomp_pools ! jkolassa May 2022: accounting for fact that pool order in CNCOL is different from CTSM - this%decomp_npools_col (n,np) = clamp_nonneg( real(cncol(nc,nz,decomp_npool_cncol_index(np)), r8) ) + this%decomp_npools_col (n,np) = max( cncol(nc,nz,decomp_npool_cncol_index(np)), 0._r8 ) !rrXbo, 10Sep2025 this%decomp_npools_1m_col(n,np) = this%decomp_npools_col(n,np) ! jkolassa May 2022: loop has to be added below if we add more biogeochemical (or soil) layers this%decomp_npools_vr_col(n,1,np) = this%decomp_npools_col(n,np) From 67cb83959ccc3226dc1b9f04772135517283bf58 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 11 Sep 2025 16:50:30 -0400 Subject: [PATCH 16/29] use FVEG_MIN from clm_varpar_shared module (CatchmentCNRst.F90) --- .../Utils/mk_restarts/CatchmentCNRst.F90 | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 index 1edb91106..f9aabeac1 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 @@ -15,14 +15,13 @@ module CatchmentCNRstMod VAR_COL_40, VAR_PFT_40, & VAR_COL_51, VAR_PFT_51, & npft_40 => NUM_PFT_CN_40, & - npft_51 => NUM_PFT_CN_51 + npft_51 => NUM_PFT_CN_51, & + FVEG_MIN use nanMod, only : nan implicit none - real, parameter :: fmin = 1.e-4 ! ignore vegetation fractions at or below this value - integer :: iclass_40(npft_40) = (/1,1,2,3,3,4,5,5,6,7,8,9,10,11,12,11,12,11,12/) integer :: iclass_51(npft_51) = (/1,1,2,3,3,4,5,5,6,7, 9,10,11, 11, 11 /) @@ -1312,17 +1311,17 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & end if - if (fveg_new > fmin) then + if (fveg_new > FVEG_MIN) then offl_cell = Id_glb(n,nv) - if(ityp_new == ityp_offl (offl_cell,nv) .and. fveg_offl (offl_cell,nv)> fmin) then + if(ityp_new == ityp_offl (offl_cell,nv) .and. fveg_offl (offl_cell,nv)> FVEG_MIN) then iv = nv ! same type fraction (primary of secondary) - else if(ityp_new == ityp_offl (offl_cell,nx) .and. fveg_offl (offl_cell,nx)> fmin) then + else if(ityp_new == ityp_offl (offl_cell,nx) .and. fveg_offl (offl_cell,nx)> FVEG_MIN) then iv = nx ! not same fraction - else if(iclass_in(ityp_new)==iclass_in(ityp_offl(offl_cell,nv)) .and. fveg_offl (offl_cell,nv)> fmin) then + else if(iclass_in(ityp_new)==iclass_in(ityp_offl(offl_cell,nv)) .and. fveg_offl (offl_cell,nv)> FVEG_MIN) then iv = nv ! primary, other type (same class) - else if(fveg_offl (offl_cell,nx)> fmin) then + else if(fveg_offl (offl_cell,nx)> FVEG_MIN) then iv = nx ! secondary, other type (same class) endif @@ -1395,7 +1394,7 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & if (nv == 2) fveg_new = CLMC_sf1(n) end if - if(fveg_new > fmin) then + if(fveg_new > FVEG_MIN) then var_pft_out(n, nz,nv, 1) = max(var_pft_out(n, nz,nv, 1),0.) var_pft_out(n, nz,nv, 2) = max(var_pft_out(n, nz,nv, 2),0.) var_pft_out(n, nz,nv, 3) = max(var_pft_out(n, nz,nv, 3),0.) From aef24b8fbbb45ba39bd6689f4e4e3f541b500d68 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 12 Sep 2025 09:12:33 -0400 Subject: [PATCH 17/29] added missing use statement for FVEG_MIN (PhotosynthesisMod.F90) --- .../GEOScatchCNCLM51_GridComp/CLM51/PhotosynthesisMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PhotosynthesisMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PhotosynthesisMod.F90 index f06f780e9..2e13df731 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PhotosynthesisMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PhotosynthesisMod.F90 @@ -17,7 +17,7 @@ module PhotosynthesisMod use clm_varctl , only : use_c13, use_c14, use_cn, use_cndv, use_fates, use_luna, use_hydrstress use clm_varctl , only : iulog use clm_varpar , only : nlevcan, nvegwcs, mxpft - use clm_varpar , only : numpft, NUM_VEG, NUM_ZON, VAR_COL, VAR_PFT + use clm_varpar , only : numpft, NUM_VEG, NUM_ZON, VAR_COL, VAR_PFT, FVEG_MIN use clm_varcon , only : namep, spval, isecspday use decompMod , only : bounds_type use QuadraticMod , only : quadratic From 6d2be4737a5181962b9e026e3d7c7e48be63e557 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 12 Sep 2025 14:21:14 -0400 Subject: [PATCH 18/29] fixed T2MMIN5D exponential moving average (GEOS_CatchCNCLM51GridComp.F90) --- .../GEOS_CatchCNCLM51GridComp.F90 | 153 +++++++++++++----- 1 file changed, 115 insertions(+), 38 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 index 9643ab116..51fc42bce 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 @@ -5530,10 +5530,11 @@ subroutine Driver ( RC ) integer, save :: istep ! model time step index integer, save :: istep_365 ! model time step index integer :: accper ! number of time steps accumulated in a period of XX days, increases from 1 to nXXd in the first XX days, - ! and remains as nXXd thereafter - integer, allocatable, dimension(:) :: ta_count - real, allocatable, dimension(:) :: TA_MIN - + ! and remains as nXXd thereafter + integer :: accper_days + real :: tmplon1 tmplon2 + logical, allocatable, dimension(:) :: mask_5amLT + integer :: AGCM_YY, AGCM_MM, AGCM_DD, AGCM_MI, AGCM_S, AGCM_HH, dofyr, AGCM_S_ofday logical, save :: first = .true. integer(INT64), save :: istep_cn = 1 ! gkw: legacy variable from offline @@ -6231,9 +6232,9 @@ subroutine Driver ( RC ) allocate(TOTDEPOS (NTILES,N_constit)) allocate(RMELT (NTILES,N_constit)) - allocate(TA_MIN (NTILES)) - allocate(ta_count (NTILES)) + allocate(mask_5amLT(NTILES)) + call ESMF_VMGetCurrent ( VM, RC=STATUS ) debugzth = .false. @@ -7122,14 +7123,54 @@ subroutine Driver ( RC ) istep = istep + 1 istep_365 = istep_365 + 1 - TA_MIN(:) = 1000. ! running mean - reset accumulation period until greater than nstep ! fzeng & gkw: may not be exactly 2m, but it is consistent with t_ref2m in CN model ! T2M10D (T10 in CLM4.5): 10-day running mean of 2-m temperature (K) ! TPREC10D (PREC10 in CLM4.5): 10-day running mean of total precipitation (mm H2O/s) ! TPREC60D (PREC60 in CLM4.5): 60-day running mean of total precipitation (mm H2O/s) - ! --------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------- + + ! For T2MMIN5D, always use surface air temperature at ~5am local time to avoid need for additional + ! restart variables. Update every full hour. -- reichle, 12 Sep 2025 + + if ( AGCM_MI==0 .and. AGCM_S==0 ) then + + tmplon1 = LT2lon( AGCM_HH, AGCM_MI, AGCM_S, 4.5 ) ! get longitude [radians] where it is 4:30am local time + tmplon2 = LT2lon( AGCM_HH, AGCM_MI, AGCM_S, 5.5 ) ! get longitude [radians] where it is 5:30am local time + + ! create mask for tiles within longitude band (4:30-5:30am local time); + ! assumes -pi <= LONS <= pi + + if (tmplon1= 24.) ) then + + err_msg = 'input local_hour falls outside allowed range of [0,24)' + + _ASSERT( .FALSE., err_msg) + + end if + + ! determine fractional UTC hour and time difference with local_hour + + UTC_hour = ( UTC_HH*3600 + UTC_MI*60 + UTC_S )/3600. ! 0 <= UTC_hour < 24 + + time_diff = local_hour - UTC_hour + + ! enforce -12. < time_diff <= 12. + + if (time_diff <= -12.) then + + time_diff = time_diff + 24. + + elseif (time_diff > 12.) then + + time_diff = time_diff - 24. + + end if + + ! determine longitude + + longitude = time_diff/24.*360. ! [degree] -180. < longitude <= +180. + + LT2lon = longitude*pi/180. ! [radians] + + end function LT2lon + +! ***************************************************************** + end module GEOS_CatchCNCLM51GridCompMod subroutine SetServices(gc, rc) From e9ae56f6c8afe5afd52d1aef91f320452bda715a Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 12 Sep 2025 14:38:31 -0400 Subject: [PATCH 19/29] refine calculation for previous commit (GEOS_CatchCNCLM51GridComp.F90) --- .../GEOS_CatchCNCLM51GridComp.F90 | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 index 51fc42bce..4107f1c04 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 @@ -7154,18 +7154,24 @@ subroutine Driver ( RC ) end if - ! compute 5-day exponential moving average (executed at daily time step from the perspective of a given tile) + ! get (inverse) weight for exponential moving average if (init_accum) then - - accper_days = min( istep/n1d, 5 ) - + + accper_days = min( (istep-1)/n1d + 1, 5 ) + else accper_days = 5 end if + ! during init_accum only, verify that from day 2 we have sensible values in T2MMIN5D (Kelvin) + + if (init_accum .and. accper_days>1) _ASSERT( (100. < T2MMIN5D) .and. (T2MMIN5D < 400.), 'error in T2MMIN5D calculation' ) + + ! compute 5-day exponential moving average (executed at daily time step from the perspective of a given tile) + T2MMIN5D(mask_5amLT) = ( (accper_days-1)*T2MMIN5D(mask_5amLT) + TA(mask_5amLT) )/accper_days end if @@ -7187,7 +7193,6 @@ subroutine Driver ( RC ) accper = min(istep,n30d) RH30D = ((accper-1)*RH30D + Qair_relative) / accper - ! (2) 60-day exponential moving average of total precipitation (mm H2O/s) accper = min(istep,n60d) TPREC60D = ((accper-1)*TPREC60D + PCU + PLS + SNO) / accper From c41c399930bf0d122de357db15f577fc5dd21423 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 12 Sep 2025 15:12:41 -0400 Subject: [PATCH 20/29] fixed syntax error in previous commit (GEOS_CatchCNCLM51GridComp.F90) --- .../GEOS_CatchCNCLM51GridComp.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 index 4107f1c04..279b31730 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 @@ -5532,7 +5532,7 @@ subroutine Driver ( RC ) integer :: accper ! number of time steps accumulated in a period of XX days, increases from 1 to nXXd in the first XX days, ! and remains as nXXd thereafter integer :: accper_days - real :: tmplon1 tmplon2 + real :: tmplon1, tmplon2 logical, allocatable, dimension(:) :: mask_5amLT integer :: AGCM_YY, AGCM_MM, AGCM_DD, AGCM_MI, AGCM_S, AGCM_HH, dofyr, AGCM_S_ofday @@ -7168,8 +7168,10 @@ subroutine Driver ( RC ) ! during init_accum only, verify that from day 2 we have sensible values in T2MMIN5D (Kelvin) - if (init_accum .and. accper_days>1) _ASSERT( (100. < T2MMIN5D) .and. (T2MMIN5D < 400.), 'error in T2MMIN5D calculation' ) - + if ( init_accum .and. (accper_days>1) ) then + _ASSERT( (100. < T2MMIN5D) .and. (T2MMIN5D < 400.), 'error in T2MMIN5D calculation' ) + end if + ! compute 5-day exponential moving average (executed at daily time step from the perspective of a given tile) T2MMIN5D(mask_5amLT) = ( (accper_days-1)*T2MMIN5D(mask_5amLT) + TA(mask_5amLT) )/accper_days From 7e32903011fe81d3eb3a50f00065298bcfd101b0 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 12 Sep 2025 15:42:21 -0400 Subject: [PATCH 21/29] fixed _ASSERT error in previous commit (GEOS_CatchCNCLM51GridComp.F90) --- .../GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 index 279b31730..b3b6a635f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 @@ -7169,7 +7169,7 @@ subroutine Driver ( RC ) ! during init_accum only, verify that from day 2 we have sensible values in T2MMIN5D (Kelvin) if ( init_accum .and. (accper_days>1) ) then - _ASSERT( (100. < T2MMIN5D) .and. (T2MMIN5D < 400.), 'error in T2MMIN5D calculation' ) + _ASSERT( all( (100. < T2MMIN5D) .and. (T2MMIN5D < 400.) ), 'error in T2MMIN5D calculation' ) end if ! compute 5-day exponential moving average (executed at daily time step from the perspective of a given tile) From cd634db08e18ffb5090aed458d6e932ad83cf443 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 12 Sep 2025 16:34:16 -0400 Subject: [PATCH 22/29] convert function LT2lon() into subroutine to make MAPL _ASSERT() work (GEOS_CatchCNCLM51GridComp.F90) --- .../GEOS_CatchCNCLM51GridComp.F90 | 22 ++++++++++--------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 index b3b6a635f..f8d45d174 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 @@ -7136,8 +7136,8 @@ subroutine Driver ( RC ) if ( AGCM_MI==0 .and. AGCM_S==0 ) then - tmplon1 = LT2lon( AGCM_HH, AGCM_MI, AGCM_S, 4.5 ) ! get longitude [radians] where it is 4:30am local time - tmplon2 = LT2lon( AGCM_HH, AGCM_MI, AGCM_S, 5.5 ) ! get longitude [radians] where it is 5:30am local time + call LT2lon( AGCM_HH, AGCM_MI, AGCM_S, 4.5, tmplon1, rc ) ! get longitude [radians] where it is 4:30am local time + call LT2lon( AGCM_HH, AGCM_MI, AGCM_S, 5.5, tmplon2, rc ) ! get longitude [radians] where it is 5:30am local time ! create mask for tiles within longitude band (4:30-5:30am local time); ! assumes -pi <= LONS <= pi @@ -9188,7 +9188,7 @@ end subroutine RUN0 ! ***************************************************************** - real function LT2lon( UTC_HH, UTC_MI, UTC_S, local_hour ) + subroutine LT2lon( UTC_HH, UTC_MI, UTC_S, local_hour, longitude, rc ) ! for UTC hour/min/sec, find the longitude where the local time is (fractional) local_hour; ! when UTC and local time are exactly 12 hours offset, always return longitude=+PI. @@ -9197,12 +9197,14 @@ real function LT2lon( UTC_HH, UTC_MI, UTC_S, local_hour ) implicit none - integer, intent(in) :: UTC_HH, UTC_MI, UTC_S ! UTC hour/min/sec - - real, intent(in) :: local_hour ! fractional local hour (e.g., 1:30pm is 13.5); range=[0,24) - - real :: longitude ! [degree] -pi < longitude <= +pi + integer, intent(in) :: UTC_HH, UTC_MI, UTC_S ! UTC hour/min/sec + real, intent(in) :: local_hour ! fractional local hour (e.g., 1:30pm is 13.5); range=[0,24) + + real, intent(out) :: longitude ! [degree] -pi < longitude <= +pi + + integer, intent(out), optional :: rc ! return code for _ASSERT() + ! ---------------------------------------------- real :: UTC_hour, time_diff @@ -9244,9 +9246,9 @@ real function LT2lon( UTC_HH, UTC_MI, UTC_S, local_hour ) longitude = time_diff/24.*360. ! [degree] -180. < longitude <= +180. - LT2lon = longitude*pi/180. ! [radians] + longitude = longitude*pi/180. ! [radians] - end function LT2lon + end subroutine LT2lon ! ***************************************************************** From 6e2ce4560f9f7e914f46719969bcb2dc02bd8a04 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 12 Sep 2025 16:52:19 -0400 Subject: [PATCH 23/29] avoid duplicate definition of index translation for C & N pools between CatchCN and CTSM (clm_varpar.F90, CNCLM_DriverMod.F90, SoilBiogeochemCarbonStateType.F90, SoilBiogeochemNitrogenStateType.F90) --- .../CLM51/CNCLM_DriverMod.F90 | 14 ++++++-------- .../CLM51/SoilBiogeochemCarbonStateType.F90 | 6 +++--- .../CLM51/SoilBiogeochemNitrogenStateType.F90 | 6 +++--- .../GEOScatchCNCLM51_GridComp/CLM51/clm_varpar.F90 | 5 +++++ 4 files changed, 17 insertions(+), 14 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 index ea1bc9975..0a7506647 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 @@ -4,8 +4,9 @@ module CNCLM_DriverMod use nanMod , only : nan use CNVegetationFacade - use clm_varpar , only : nlevsno, nlevmaxurbgrnd, num_veg, num_zon, CN_zone_weight, & - var_col, var_pft, nlevgrnd, numpft, ndecomp_pools, FVEG_MIN + use clm_varpar , only : nlevsno, nlevmaxurbgrnd, num_veg, num_zon, CN_zone_weight + use clm_varpar , only : var_col, var_pft, nlevgrnd, numpft, ndecomp_pools, FVEG_MIN + use clm_varpar , only : decomp_cpool_cncol_index, decomp_npool_cncol_index use clm_varcon , only : grav, denh2o use clm_time_manager , only : is_first_step, get_nstep use decompMod @@ -586,11 +587,6 @@ subroutine CN_exit(nch,ityp,fveg,cncol,cnpft) integer :: n, p, nv, nc, nz, np, nd - ! map between order of C & N pools in CNCOL and CTSM - - integer, dimension(8) :: decomp_cpool_cncol_index = (/ 3, 4, 5, 2, 10, 11, 12, 13 /) - integer, dimension(8) :: decomp_npool_cncol_index = (/ 18, 19, 20, 17, 25, 26, 27, 28 /) - !---------------- !rrXbo10Sep2025 ! Clean slate: zero arrays so any unassigned slots are safe in the restart @@ -608,6 +604,8 @@ subroutine CN_exit(nch,ityp,fveg,cncol,cnpft) cncol(nc,nz, 1) = soilbiogeochem_carbonstate_inst%ctrunc_vr_col(n,1) + ! C and N decomposition pools are stored in cncol entries 2:5, 10:13, 17:20, and 25:28 + do nd = 1,ndecomp_pools ! jkolassa: accounting for fact that pool order in CNCOL is different from CTSM cncol(nc,nz,decomp_cpool_cncol_index(nd)) = soilbiogeochem_carbonstate_inst%decomp_cpools_vr_col( n,1,nd) @@ -626,12 +624,12 @@ subroutine CN_exit(nch,ityp,fveg,cncol,cnpft) cncol(nc,nz,15) = soilbiogeochem_carbonstate_inst%totlitc_col( n ) cncol(nc,nz,16) = soilbiogeochem_nitrogenstate_inst%ntrunc_vr_col( n,1) - cncol(nc,nz,21) = bgc_vegetation_inst%n_products_inst%prod100_grc( nc) * CN_zone_weight(nz) cncol(nc,nz,22) = bgc_vegetation_inst%n_products_inst%prod10_grc( nc) * CN_zone_weight(nz) cncol(nc,nz,23) = bgc_vegetation_inst%cnveg_nitrogenstate_inst%seedn_grc( nc) * CN_zone_weight(nz) cncol(nc,nz,24) = soilbiogeochem_nitrogenstate_inst%sminn_vr_col( n,1) + cncol(nc,nz,29) = bgc_vegetation_inst%cnveg_nitrogenstate_inst%totn_col( n ) cncol(nc,nz,30) = soilbiogeochem_state_inst%fpg_col( n ) cncol(nc,nz,31) = bgc_vegetation_inst%cnveg_state_inst%annsum_counter_col( n ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemCarbonStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemCarbonStateType.F90 index 73475539b..bc5f03917 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemCarbonStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemCarbonStateType.F90 @@ -5,8 +5,9 @@ module SoilBiogeochemCarbonStateType use abortutils , only : endrun use nanMod , only : nan use clm_varpar , only : ndecomp_cascade_transitions, ndecomp_pools, nlevcan - use clm_varpar , only : nlevdecomp_full, nlevdecomp, nlevsoi, & - NUM_ZON, VAR_COL + use clm_varpar , only : nlevdecomp_full, nlevdecomp, nlevsoi + use clm_varpar , only : NUM_ZON, VAR_COL + use clm_varpar , only : decomp_cpool_cncol_index use clm_varcon , only : spval, ispval, dzsoi_decomp, zisoi, zsoi, c3_r2 use clm_varctl , only : iulog, use_vertsoilc, use_fates, use_soil_matrixcn, use_century_decomp use decompMod , only : bounds_type @@ -84,7 +85,6 @@ subroutine Init(this, bounds, nch, cncol) ! !LOCAL VARIABLES: integer :: begc,endc integer :: n, nc, nz, np - integer, dimension(8) :: decomp_cpool_cncol_index = (/ 3, 4, 5, 2, 10, 11, 12, 13 /) !----------------------------------- begc = bounds%begc ; endc = bounds%endc diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 index 3f7b48c53..6341de071 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/SoilBiogeochemNitrogenStateType.F90 @@ -5,8 +5,9 @@ module SoilBiogeochemNitrogenStateType use abortutils , only : endrun use nanMod , only : nan use clm_varpar , only : ndecomp_cascade_transitions, ndecomp_pools, nlevcan - use clm_varpar , only : nlevdecomp_full, nlevdecomp, nlevsoi, & - NUM_ZON, VAR_COL + use clm_varpar , only : nlevdecomp_full, nlevdecomp, nlevsoi + use clm_varpar , only : NUM_ZON, VAR_COL + use clm_varpar , only : decomp_npool_cncol_index use clm_varcon , only : spval, dzsoi_decomp, zisoi use clm_varctl , only : use_nitrif_denitrif, use_vertsoilc, use_century_decomp, use_soil_matrixcn !, iulog !rrXbo 10Sep2025 use decompMod , only : bounds_type @@ -129,7 +130,6 @@ subroutine Init(this, bounds, nch, cncol) ! !LOCAL VARIABLES: integer :: begc,endc integer :: n, nc, nz, np, l, c - integer, dimension(8) :: decomp_npool_cncol_index = (/ 18, 19, 20, 17,25, 26, 27, 28 /) logical :: no_cn51_rst = .false. !----------------------------------- diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/clm_varpar.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/clm_varpar.F90 index b566db6e2..ce3126631 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/clm_varpar.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/clm_varpar.F90 @@ -40,6 +40,11 @@ module clm_varpar integer, public :: ndecomp_cascade_transitions integer, public :: ndecomp_cascade_outtransitions + ! map between order of C & N pools in CNCOL and CTSM; *MUST* have ndecomp_pools=8!!! + + integer, public, dimension(8) :: decomp_cpool_cncol_index = (/ 3, 4, 5, 2, 10, 11, 12, 13 /) + integer, public, dimension(8) :: decomp_npool_cncol_index = (/ 18, 19, 20, 17, 25, 26, 27, 28 /) + ! for soil matrix integer, public :: ndecomp_pools_vr ! total number of pools ndecomp_pools*vertical levels From 999429858a1e97ebd90366e156c43c97ce54e6e3 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 12 Sep 2025 17:19:38 -0400 Subject: [PATCH 24/29] added comment about different decomp pool orders in CatchCN and CTSM (CNCLM_DriverMod.F90) --- .../GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 index 0a7506647..e7e754fce 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_DriverMod.F90 @@ -565,7 +565,7 @@ subroutine CN_exit(nch,ityp,fveg,cncol,cnpft) ! OUTPUT real, dimension(nch,num_zon,var_col), intent(out) :: cncol ! column-level restart variables - real, dimension(nch,num_zon,num_veg,var_pft), intent(out) :: cnpft ! PFT-level restart variables + real, dimension(nch,num_zon,num_veg,var_pft), intent(out) :: cnpft ! PFT-level restart variables ! NOTE: order of dims DIFFERS from that of ityp & fveg !!! ! LOCAL From 4442cd482791d20cd8eea7f3d4ff0cd971c6ac61 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 12 Sep 2025 17:38:22 -0400 Subject: [PATCH 25/29] fixed bug in 5-day moving average of snow depth (GEOS_CatchCNCLM51GridComp.F90) --- .../GEOS_CatchCNCLM51GridComp.F90 | 25 ++++++++++--------- 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 index f8d45d174..2ef245283 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 @@ -1872,9 +1872,9 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddInternalSpec(GC ,& - LONG_NAME = '5-day_exp_moving_avg_of_CN_sum_for_snow_depth',& + LONG_NAME = '5-day_exp_moving_avg_of_snow_depth',& UNITS = 'm' ,& - SHORT_NAME = 'SNDZM5D' ,& + SHORT_NAME = 'SNDZM5D' ,& !rrXbo, 10Sep2025 -- SHOULD BE RENAMED IN RESTART TO "SNDZ5D" DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RESTART = MAPL_RestartOptional ,& @@ -5993,12 +5993,13 @@ subroutine Driver ( RC ) ! set number of time steps within a XX-day/hour period for 2m temperature XX-day/hour "running mean" ! -------------------------------------------------------------------------------------------------- - n1d = 86400/dt - n5d = 5*86400/dt - n10d = 10*86400/dt - n30d = 30*86400/dt - n60d = 60*86400/dt + n1d = 86400/dt + n5d = 5*86400/dt + n10d = 10*86400/dt + n30d = 30*86400/dt + n60d = 60*86400/dt n365d = 365*86400/dt + ! fzeng: this is done in such way to exclude istep in the restart file if(init_accum) then istep = 0 ! set model time step index to 0 when begin to accumulate the cumulative variables, fzeng, 21 Apr 2017 @@ -7183,17 +7184,17 @@ subroutine Driver ( RC ) ! (1) 5-day exponential moving average of snow depth accper = min(istep,n5d) - SNDZM5D = ((accper-1)*SNDZM5D + SNDZM) / accper + SNDZM5D = ((accper-1)*SNDZM5D + sum(SNDZN,1) ) / accper ! (1) 10-day exponential moving average of 2-m temperature (K) and total precipitation (mm H2O/s) accper = min(istep,n10d) - T2M10D = ((accper-1)*T2M10D + TA) / accper + T2M10D = ((accper-1)*T2M10D + TA ) / accper TPREC10D = ((accper-1)*TPREC10D + PCU + PLS + SNO) / accper - TG10D = ((accper-1)*TG10D + TG(:,1)) / accper + TG10D = ((accper-1)*TG10D + TG(:,1) ) / accper ! (2) 30-day exponential moving average of relative humidity [%] accper = min(istep,n30d) - RH30D = ((accper-1)*RH30D + Qair_relative) / accper + RH30D = ((accper-1)*RH30D + Qair_relative ) / accper ! (2) 60-day exponential moving average of total precipitation (mm H2O/s) accper = min(istep,n60d) @@ -7201,7 +7202,7 @@ subroutine Driver ( RC ) else - SNDZM5D = (( n5d-1)*SNDZM5D + SNDZM ) / n5d + SNDZM5D = (( n5d-1)*SNDZM5D + sum(SNDZN,1) ) / n5d T2M10D = ((n10d-1)*T2M10D + TA ) / n10d TG10D = ((n10d-1)*TG10D + TG(:,1) ) / n10d TPREC10D = ((n10d-1)*TPREC10D + PCU + PLS + SNO) / n10d From b9d6ad4613bf99140c1cc5b7a62f7f5a4bf6e567 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 12 Sep 2025 17:44:11 -0400 Subject: [PATCH 26/29] change comment so it does not trip up GNU (GEOS_CatchCNCLM51GridComp.F90) --- .../GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 index 2ef245283..bb8c1435a 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 @@ -9204,7 +9204,7 @@ subroutine LT2lon( UTC_HH, UTC_MI, UTC_S, local_hour, longitude, rc ) real, intent(out) :: longitude ! [degree] -pi < longitude <= +pi - integer, intent(out), optional :: rc ! return code for _ASSERT() + integer, intent(out), optional :: rc ! return code for MAPL ASSERT() ! ---------------------------------------------- From acfe37bef2ec497d9045fa838af2da14a073b27a Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sat, 13 Sep 2025 07:58:58 -0400 Subject: [PATCH 27/29] use WHERE instead of indexing with logical array to make GNU happy (GEOS_CatchCNCLM51GridComp.F90) --- .../GEOS_CatchCNCLM51GridComp.F90 | 64 ++++++++++--------- 1 file changed, 33 insertions(+), 31 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 index bb8c1435a..dcc0b3549 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/GEOS_CatchCNCLM51GridComp.F90 @@ -7158,9 +7158,9 @@ subroutine Driver ( RC ) ! get (inverse) weight for exponential moving average if (init_accum) then - + accper_days = min( (istep-1)/n1d + 1, 5 ) - + else accper_days = 5 @@ -7172,44 +7172,46 @@ subroutine Driver ( RC ) if ( init_accum .and. (accper_days>1) ) then _ASSERT( all( (100. < T2MMIN5D) .and. (T2MMIN5D < 400.) ), 'error in T2MMIN5D calculation' ) end if - + ! compute 5-day exponential moving average (executed at daily time step from the perspective of a given tile) - T2MMIN5D(mask_5amLT) = ( (accper_days-1)*T2MMIN5D(mask_5amLT) + TA(mask_5amLT) )/accper_days + WHERE (mask_5amLT) + T2MMIN5D = ((accper_days-1)*T2MMIN5D + TA)/accper_days + END WHERE end if - - if(init_accum) then - - ! (1) 5-day exponential moving average of snow depth - accper = min(istep,n5d) - SNDZM5D = ((accper-1)*SNDZM5D + sum(SNDZN,1) ) / accper - ! (1) 10-day exponential moving average of 2-m temperature (K) and total precipitation (mm H2O/s) - accper = min(istep,n10d) - T2M10D = ((accper-1)*T2M10D + TA ) / accper - TPREC10D = ((accper-1)*TPREC10D + PCU + PLS + SNO) / accper - TG10D = ((accper-1)*TG10D + TG(:,1) ) / accper + if(init_accum) then - ! (2) 30-day exponential moving average of relative humidity [%] - accper = min(istep,n30d) - RH30D = ((accper-1)*RH30D + Qair_relative ) / accper + ! (1) 5-day exponential moving average of snow depth + accper = min(istep,n5d) + SNDZM5D = ((accper-1)*SNDZM5D + sum(SNDZN,1) ) / accper - ! (2) 60-day exponential moving average of total precipitation (mm H2O/s) - accper = min(istep,n60d) - TPREC60D = ((accper-1)*TPREC60D + PCU + PLS + SNO) / accper + ! (1) 10-day exponential moving average of 2-m temperature (K) and total precipitation (mm H2O/s) + accper = min(istep,n10d) + T2M10D = ((accper-1)*T2M10D + TA ) / accper + TPREC10D = ((accper-1)*TPREC10D + PCU + PLS + SNO) / accper + TG10D = ((accper-1)*TG10D + TG(:,1) ) / accper - else - - SNDZM5D = (( n5d-1)*SNDZM5D + sum(SNDZN,1) ) / n5d - T2M10D = ((n10d-1)*T2M10D + TA ) / n10d - TG10D = ((n10d-1)*TG10D + TG(:,1) ) / n10d - TPREC10D = ((n10d-1)*TPREC10D + PCU + PLS + SNO) / n10d - RH30D = ((n30d-1)*RH30D + Qair_relative ) / n30d - TPREC60D = ((n60d-1)*TPREC60D + PCU + PLS + SNO) / n60d - - endif + ! (2) 30-day exponential moving average of relative humidity [%] + accper = min(istep,n30d) + RH30D = ((accper-1)*RH30D + Qair_relative ) / accper + + ! (2) 60-day exponential moving average of total precipitation (mm H2O/s) + accper = min(istep,n60d) + TPREC60D = ((accper-1)*TPREC60D + PCU + PLS + SNO) / accper + + else + + SNDZM5D = (( n5d-1)*SNDZM5D + sum(SNDZN,1) ) / n5d + T2M10D = ((n10d-1)*T2M10D + TA ) / n10d + TG10D = ((n10d-1)*TG10D + TG(:,1) ) / n10d + TPREC10D = ((n10d-1)*TPREC10D + PCU + PLS + SNO) / n10d + RH30D = ((n30d-1)*RH30D + Qair_relative ) / n30d + TPREC60D = ((n60d-1)*TPREC60D + PCU + PLS + SNO) / n60d + + endif ! get CO2 From a7a0f1b31037d936fc7f05c55905ede75b97de10 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 15 Sep 2025 17:31:54 -0400 Subject: [PATCH 28/29] removed re-setting of carbon restart variables to zero (CatchmentCNRst.F90) --- .../Utils/mk_restarts/CatchmentCNRst.F90 | 232 ++++++++++-------- 1 file changed, 128 insertions(+), 104 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 index f9aabeac1..abe13160a 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 @@ -22,8 +22,8 @@ module CatchmentCNRstMod implicit none - integer :: iclass_40(npft_40) = (/1,1,2,3,3,4,5,5,6,7,8,9,10,11,12,11,12,11,12/) - integer :: iclass_51(npft_51) = (/1,1,2,3,3,4,5,5,6,7, 9,10,11, 11, 11 /) + integer, parameter :: iclass_40(npft_40) = (/1,1,2,3,3,4,5,5,6,7,8,9,10,11,12,11,12,11,12/) + integer, parameter :: iclass_51(npft_51) = (/1,1,2,3,3,4,5,5,6,7, 9,10,11, 11, 11 /) integer, dimension(:), allocatable :: iclass @@ -111,28 +111,33 @@ module CatchmentCNRstMod contains function CatchmentCNRst_create(filename, cnclm, time, rc) result (catch) - type(CatchmentCNRst) :: catch - character(*), intent(in) :: filename - character(*), intent(in) :: cnclm - character(*), intent(in) :: time + + type(CatchmentCNRst) :: catch + + character(*), intent(in) :: filename + character(*), intent(in) :: cnclm + character(*), intent(in) :: time + integer, optional, intent(out) :: rc - integer :: status + + integer :: status type(Netcdf4_fileformatter) :: formatter - integer :: filetype, ntiles, unit - integer :: j, dim1,dim2, myid, mpierr - type(Variable), pointer :: myVariable - character(len=:), pointer :: dname - type(FileMetadata) :: meta - character(len=256) :: Iam = "CatchmentCNRst_create" - integer :: nt_col, blk_col, off - real :: N_MAX - integer :: joff - integer :: n_fixed - + integer :: filetype, ntiles, unit + integer :: j, dim1,dim2, myid, mpierr + type(Variable), pointer :: myVariable + character(len=:), pointer :: dname + type(FileMetadata) :: meta + character(len=256) :: Iam = "CatchmentCNRst_create" + + !rrXbo, 10Sep2025 + !integer :: nt_col, blk_col, off + !real :: N_MAX + !integer :: joff + !integer :: n_fixed call MAPL_NCIOGetFileType(filename, filetype, __RC__) if (filetype /= 0) then - _ASSERT( .false., "CatchmentCN only support nc4 file restart") + _ASSERT( .false., "CatchmentCN only supports nc4 file restart") endif call MPI_COMM_RANK( MPI_COMM_WORLD, myid, mpierr ) @@ -230,59 +235,60 @@ function CatchmentCNRst_create(filename, cnclm, time, rc) result (catch) call MAPL_VarRead( formatter, "RUNSURF", catch%RUNSURF , __RC__) ! rreichle, bug fix 9 Sep 2025: added for CLM51 - endif !! <— close CLM40/CLM51 branch - - ! strip NaNs - where (.not.(catch%CNCOL == catch%CNCOL)) catch%CNCOL = 0.0 - where (.not.(catch%CNPFT == catch%CNPFT)) catch%CNPFT = 0.0 + endif - ! Targeted sanitize for N in CNCOL: 24=sminn, 36=NO3, 37=NH4 - ! Without this section model can't run values are junk! - - N_MAX = 1.0e5 ! conservative guardrail - - ! Slot 24: sminn - joff = (24-1)*nzone - where (catch%CNCOL(:, joff+1:joff+nzone) /= catch%CNCOL(:, joff+1:joff+nzone)) - catch%CNCOL(:, joff+1:joff+nzone) = 0.0 - end where - where (catch%CNCOL(:, joff+1:joff+nzone) < 0.0) - catch%CNCOL(:, joff+1:joff+nzone) = 0.0 - end where - where (catch%CNCOL(:, joff+1:joff+nzone) > N_MAX) - catch%CNCOL(:, joff+1:joff+nzone) = 0.0 - end where - - ! Slot 36: smin_no3 - joff = (36-1)*nzone - where (catch%CNCOL(:, joff+1:joff+nzone) /= catch%CNCOL(:, joff+1:joff+nzone)) - catch%CNCOL(:, joff+1:joff+nzone) = 0.0 - end where - where (catch%CNCOL(:, joff+1:joff+nzone) < 0.0) - catch%CNCOL(:, joff+1:joff+nzone) = 0.0 - end where - where (catch%CNCOL(:, joff+1:joff+nzone) > N_MAX) - catch%CNCOL(:, joff+1:joff+nzone) = 0.0 - end where - - ! Slot 37: smin_nh4 - joff = (37-1)*nzone - where (catch%CNCOL(:, joff+1:joff+nzone) /= catch%CNCOL(:, joff+1:joff+nzone)) - catch%CNCOL(:, joff+1:joff+nzone) = 0.0 - end where - where (catch%CNCOL(:, joff+1:joff+nzone) < 0.0) - catch%CNCOL(:, joff+1:joff+nzone) = 0.0 - end where - where (catch%CNCOL(:, joff+1:joff+nzone) > N_MAX) - catch%CNCOL(:, joff+1:joff+nzone) = 0.0 - end where - - n_fixed = count(.not.(catch%CNCOL(:, (24-1)*nzone+1:(24)*nzone) == catch%CNCOL(:, (24-1)*nzone+1:(24)*nzone))) & - + count(catch%CNCOL(:, (24-1)*nzone+1:(24)*nzone) < 0.0) & - + count(catch%CNCOL(:, (24-1)*nzone+1:(24)*nzone) > N_MAX) - if (n_fixed > 0) write(*,*) 'RSTCHK sanitize CNCOL(24): fixed=', n_fixed - - endif !! <— closes myid + !rrXbo, 10Sep2025 + ! ! strip NaNs + ! where (.not.(catch%CNCOL == catch%CNCOL)) catch%CNCOL = 0.0 + ! where (.not.(catch%CNPFT == catch%CNPFT)) catch%CNPFT = 0.0 + ! + ! ! Targeted sanitize for N in CNCOL: 24=sminn, 36=NO3, 37=NH4 + ! ! Without this section model can't run values are junk! + ! + ! N_MAX = 1.0e5 ! conservative guardrail + ! + ! ! Slot 24: sminn + ! joff = (24-1)*nzone + ! where (catch%CNCOL(:, joff+1:joff+nzone) /= catch%CNCOL(:, joff+1:joff+nzone)) + ! catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + ! end where + ! where (catch%CNCOL(:, joff+1:joff+nzone) < 0.0) + ! catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + ! end where + ! where (catch%CNCOL(:, joff+1:joff+nzone) > N_MAX) + ! catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + ! end where + ! + ! ! Slot 36: smin_no3 + ! joff = (36-1)*nzone + ! where (catch%CNCOL(:, joff+1:joff+nzone) /= catch%CNCOL(:, joff+1:joff+nzone)) + ! catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + ! end where + ! where (catch%CNCOL(:, joff+1:joff+nzone) < 0.0) + ! catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + ! end where + ! where (catch%CNCOL(:, joff+1:joff+nzone) > N_MAX) + ! catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + ! end where + ! + ! ! Slot 37: smin_nh4 + ! joff = (37-1)*nzone + ! where (catch%CNCOL(:, joff+1:joff+nzone) /= catch%CNCOL(:, joff+1:joff+nzone)) + ! catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + ! end where + ! where (catch%CNCOL(:, joff+1:joff+nzone) < 0.0) + ! catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + ! end where + ! where (catch%CNCOL(:, joff+1:joff+nzone) > N_MAX) + ! catch%CNCOL(:, joff+1:joff+nzone) = 0.0 + ! end where + ! + ! n_fixed = count(.not.(catch%CNCOL(:, (24-1)*nzone+1:(24)*nzone) == catch%CNCOL(:, (24-1)*nzone+1:(24)*nzone))) & + ! + count(catch%CNCOL(:, (24-1)*nzone+1:(24)*nzone) < 0.0) & + ! + count(catch%CNCOL(:, (24-1)*nzone+1:(24)*nzone) > N_MAX) + ! if (n_fixed > 0) write(*,*) 'RSTCHK sanitize CNCOL(24): fixed=', n_fixed + + endif ! myid==0 call formatter%close() @@ -1168,7 +1174,9 @@ subroutine re_tile(this, InTileFile, OutBcsDir, OutTileFile, surflay, rc) end do end do - + !rrXbo, 10Sep2025 -- the following lines were in original file, but not sure why NaNs should be turned into zeros; omit for now + !where(isnan(var_off_pft)) var_off_pft = 0. + !where(var_off_pft /= var_off_pft) var_off_pft = 0. ! should do same as previous line, but keep for now just in case print *, 'calculating regridded carbn' @@ -1222,9 +1230,11 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & integer :: N, STATUS, nv, nx, offl_cell, ityp_new, i, j, nz, iv, jj real :: fveg_new character(256) :: Iam = "write_regridded_carbon" - ! Column slots (CNCOL) that must be >= 0 + + + ! Column slots (CNCOL) that must be >= 0 integer, parameter :: nonneg_col(*) = (/ 1,2,3,4,5, 6, 10,11,12,13, 14,15,16, 21,22,23,24,25,26,27,28,29, 30,31,32,33,34,35,36,37 /) - ! PFT slots (CNPFT) that must be >= 0 (DO NOT include 76–79: veg water potentials) + ! PFT slots (CNPFT) that must be >= 0 (DO NOT include 76–79: veg water potentials) integer, parameter :: nonneg_pft(*) = (/ & 1,2,3,4,5, 6,7,8,9,10, 11,12,13,14,15, 16,17,18,19,20, 21,22,23, & 24,25,26,27, & @@ -1233,7 +1243,7 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & 69,70,71,72,73,74, & ! canopy geometry 75, & ! plant_ndemand 80,81 & ! litfall sums - /) + /) allocate (CLMC_pf1(NTILES)) @@ -1269,9 +1279,9 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & allocate (var_col_out (1: NTILES, 1 : nzone,1 : var_col)) allocate (var_pft_out (1: NTILES, 1 : nzone,1 : nveg, 1 : var_pft)) - var_col_out = 0. - !var_pft_out = NaN - var_pft_out = 0.0 + var_col_out = 0. !rrXbo, 10Sep2025 -- why not initialized to NaN (matching var_pft_out)? + var_pft_out = NaN + !var_pft_out = 0.0 !rrXbo, 10Sep2025 OUT_TILE : DO N = 1, NTILES @@ -1347,18 +1357,21 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & end do NVLOOP2 - ! reset carbon if negative < 10g + ! reset if too little (or negative) total carbon: totc_col < 10 (gC/m2) ! ------------------------ NZLOOP : do nz = 1, nzone - if(var_col_out (n, nz,14) < 10.) then + if(var_col_out (n, nz,14) < 10.) then ! totc_col<10 var_col_out(n, nz, 1) = max(var_col_out(n, nz, 1), 0.) var_col_out(n, nz, 2) = max(var_col_out(n, nz, 2), 0.) var_col_out(n, nz, 3) = max(var_col_out(n, nz, 3), 0.) var_col_out(n, nz, 4) = max(var_col_out(n, nz, 4), 0.) var_col_out(n, nz, 5) = max(var_col_out(n, nz, 5), 0.) + + ! skip 6-9 [why? - reichle, 15 Sep 2025] + var_col_out(n, nz,10) = max(var_col_out(n, nz,10), 0.) var_col_out(n, nz,11) = max(var_col_out(n, nz,11), 0.) var_col_out(n, nz,12) = max(var_col_out(n, nz,12), 0.) @@ -1377,6 +1390,9 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & var_col_out(n, nz,28) = max(var_col_out(n, nz,28), 1.) var_col_out(n, nz,29) = max(var_col_out(n, nz,29), 0.) + ! skip 30-40 (CNCLM40) [why? - reichle, 15 Sep 2025] + ! skip 30-37 (CNCLM51) [why? - reichle, 15 Sep 2025] + NVLOOP3 : do nv = 1,nveg if (this%isCLM40) then if (nv == 1) ityp_new = CLMC_pt1(n) @@ -1443,15 +1459,15 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & var_pft_out(n, nz,nv,48) = max(var_pft_out(n, nz,nv,48),0.) var_pft_out(n, nz,nv,49) = max(var_pft_out(n, nz,nv,49),0.) var_pft_out(n, nz,nv,50) = max(var_pft_out(n, nz,nv,50),0.) - var_pft_out(n, nz,nv,51) = max(var_pft_out(n, nz,nv, 5)/500.,0.) - var_pft_out(n, nz,nv,52) = max(var_pft_out(n, nz,nv,52),0.) + var_pft_out(n, nz,nv,51) = max(var_pft_out(n, nz,nv, 5)/500.,0.) ! limit: deadstemn_patch <= deadstemc_patch/500. + var_pft_out(n, nz,nv,52) = max(var_pft_out(n, nz,nv,52),0.) ! no limit? deadstemc_storage_patch var_pft_out(n, nz,nv,53) = max(var_pft_out(n, nz,nv,53),0.) var_pft_out(n, nz,nv,54) = max(var_pft_out(n, nz,nv,54),0.) var_pft_out(n, nz,nv,55) = max(var_pft_out(n, nz,nv,55),0.) - var_pft_out(n, nz,nv,56) = max(var_pft_out(n, nz,nv,56),0.) - var_pft_out(n, nz,nv,57) = max(var_pft_out(n, nz,nv,13)/25.,0.) - var_pft_out(n, nz,nv,58) = max(var_pft_out(n, nz,nv,14)/25.,0.) - var_pft_out(n, nz,nv,59) = max(var_pft_out(n, nz,nv,59),0.) + var_pft_out(n, nz,nv,56) = max(var_pft_out(n, nz,nv,56),0.) + var_pft_out(n, nz,nv,57) = max(var_pft_out(n, nz,nv,13)/25.,0.) ! limit: leafn_patch <= leafc_patch/25. + var_pft_out(n, nz,nv,58) = max(var_pft_out(n, nz,nv,14)/25.,0.) ! limit: leafn_storage_patch <= leafc_storage_patch/25. + var_pft_out(n, nz,nv,59) = max(var_pft_out(n, nz,nv,59),0.) ! var_pft_out(n, nz,nv,60) = max(var_pft_out(n, nz,nv,60),0.) var_pft_out(n, nz,nv,61) = max(var_pft_out(n, nz,nv,61),0.) var_pft_out(n, nz,nv,62) = max(var_pft_out(n, nz,nv,62),0.) @@ -1467,6 +1483,7 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & var_pft_out(n, nz,nv,72) = max(var_pft_out(n, nz,nv,72),0.) ! HTOP are heights var_pft_out(n, nz,nv,73) = max(var_pft_out(n, nz,nv,73),0.) var_pft_out(n, nz,nv,74) = max(var_pft_out(n, nz,nv,74),0.) + !if(this%isCLM45) var_pft_out(n, nz,nv,75) = max(var_pft_out(n, nz,nv,75),0.) if(this%isCLM51) then var_pft_out(n, nz,nv,75) = max(var_pft_out(n, nz,nv,75),0.) @@ -1488,22 +1505,28 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & endif ! end carbon check end do NZLOOP ! end zone loop - ! --- Enforce non-negativity on known non-negative fields (all tiles) --- - do jj = 1, size(nonneg_col) - do nz = 1, nzone - var_col_out(:,nz,nonneg_col(jj)) = max( var_col_out(:,nz,nonneg_col(jj)), 0.0 ) - end do - end do + !rrXbo, 10Sep2025 + !list of "non-neg" variables applied here differs from that in preceding block: + ! var_col_out: 6 is reset + ! var_col_out: 17-20 is *not* reset + ! var_col_out: 30-37 is reset + ! + ! ! --- Enforce non-negativity on known non-negative fields (all tiles) --- + ! do jj = 1, size(nonneg_col) + ! do nz = 1, nzone + ! var_col_out(:,nz,nonneg_col(jj)) = max( var_col_out(:,nz,nonneg_col(jj)), 0.0 ) + ! end do + ! end do + ! + ! do j = 1, size(nonneg_pft) + ! do nz = 1, nzone + ! do nv = 1, nveg + ! var_pft_out(:,nz,nv,nonneg_pft(j)) = max( var_pft_out(:,nz,nv,nonneg_pft(j)), 0.0 ) + ! end do + ! end do + ! end do - do j = 1, size(nonneg_pft) - do nz = 1, nzone - do nv = 1, nveg - var_pft_out(:,nz,nv,nonneg_pft(j)) = max( var_pft_out(:,nz,nv,nonneg_pft(j)), 0.0 ) - end do - end do - end do - - ! Update dayx variable var_pft_out (:,:,28) + ! Update dayx variable var_pft_out (:,:,28) [daylight duration in seconds] do j = 28, 28 ! 1,VAR_PFT var_pft_out (:,:,:,28) do nv = 1,nveg do nz = 1,nzone @@ -1643,9 +1666,10 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & end do OUT_TILE - ! --- Final scrub: only finite, plausible magnitudes make it to the restart - where (var_col_out /= var_col_out .or. abs(var_col_out) > 1.0e20) var_col_out = 0.0 - where (var_pft_out /= var_pft_out .or. abs(var_pft_out) > 1.0e20) var_pft_out = 0.0 + !rrXbo, 10Sep2025 + ! ! --- Final scrub: only finite, plausible magnitudes make it to the restart + ! where (var_col_out /= var_col_out .or. abs(var_col_out) > 1.0e20) var_col_out = 0.0 + ! where (var_pft_out /= var_pft_out .or. abs(var_pft_out) > 1.0e20) var_pft_out = 0.0 i = 1 From 3717bd2306316935bf7582e7cace8454d935f125 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Tue, 16 Sep 2025 09:17:36 -0400 Subject: [PATCH 29/29] fixed indexing example (README_CN51.md) --- .../GEOScatchCNCLM51_GridComp/README_CN51.md | 64 ++++++++++--------- 1 file changed, 33 insertions(+), 31 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/README_CN51.md b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/README_CN51.md index 17c2405a3..e86572d8f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/README_CN51.md +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/README_CN51.md @@ -269,38 +269,40 @@ Files that are specific to Catchment-CN5.1 are in the ./GEOScatchCNCLM51\_GridCo | WaterStateType.F90 | C1, C4 | | WaterType.F90 | C1, C4 | -**Table B2: Illustration of the layout of a CTSM variable array** In this example, tile 1 contains PFTs 7 and 10 and has data values for the corresponding entries. Only the gray shaded data is present in the CTSM variable array, the other columns of the table serve to illustrate the layout. +**Table B2: Illustration of the layout of a CTSM variable array** In this example, tile 1 contains PFTs 7 and 10 and has data values for the corresponding entries. Only the "Data" column is present in the CTSM variable array, the other columns of the table serve to illustrate the layout. | Array Index | Tile | Zone | PFT | Data | |-------------|------|------|-----|------| -| 1 | 1 | 1 | 1 | NaN | -| 2 | 1 | 1 | 2 | NaN | -| 3 | 1 | 1 | 3 | NaN | -| 4 | 1 | 1 | 4 | NaN | -| 5 | 1 | 1 | 5 | NaN | -| 6 | 1 | 1 | 6 | NaN | -| 7 | 1 | 1 | 7 | Some data | -| 8 | 1 | 1 | 8 | NaN | -| 9 | 1 | 1 | 9 | NaN | -| 10 | 1 | 1 | 10 | Some data | -| 11 | 1 | 1 | 11 | NaN | -| 12 | 1 | 1 | 12 | NaN | -| 13 | 1 | 1 | 13 | NaN | -| 14 | 1 | 1 | 14 | NaN | -| 15 | 1 | 1 | 15 | NaN | -| 16 | 1 | 2 | 1 | NaN | -| 17 | 1 | 2 | 2 | NaN | -| 18 | 1 | 2 | 3 | NaN | -| 19 | 1 | 2 | 4 | NaN | -| 20 | 1 | 2 | 5 | NaN | -| 21 | 1 | 2 | 6 | NaN | -| 22 | 1 | 2 | 7 | Some data | -| 23 | 1 | 2 | 8 | NaN | -| 24 | 1 | 2 | 9 | NaN | -| 25 | 1 | 2 | 10 | Some data | -| 26 | 1 | 2 | 11 | NaN | -| 27 | 1 | 2 | 12 | NaN | -| 28 | 1 | 2 | 13 | NaN | -| 29 | 1 | 2 | 14 | NaN | -| 30 | 1 | 2 | 15 | NaN | +| 1 | 1 | 1 | 0 | NaN | +| 2 | 1 | 1 | 1 | NaN | +| 3 | 1 | 1 | 2 | NaN | +| 4 | 1 | 1 | 3 | NaN | +| 5 | 1 | 1 | 4 | NaN | +| 6 | 1 | 1 | 5 | NaN | +| 7 | 1 | 1 | 6 | NaN | +| 8 | 1 | 1 | 7 | Some data | +| 9 | 1 | 1 | 8 | NaN | +| 10 | 1 | 1 | 9 | NaN | +| 11 | 1 | 1 | 10 | Some data | +| 12 | 1 | 1 | 11 | NaN | +| 13 | 1 | 1 | 12 | NaN | +| 14 | 1 | 1 | 13 | NaN | +| 15 | 1 | 1 | 14 | NaN | +| 16 | 1 | 1 | 15 | NaN | +| 17 | 1 | 2 | 0 | NaN | +| 18 | 1 | 2 | 1 | NaN | +| 19 | 1 | 2 | 2 | NaN | +| 20 | 1 | 2 | 3 | NaN | +| 21 | 1 | 2 | 4 | NaN | +| 22 | 1 | 2 | 5 | NaN | +| 23 | 1 | 2 | 6 | NaN | +| 24 | 1 | 2 | 7 | Some data | +| 25 | 1 | 2 | 8 | NaN | +| 26 | 1 | 2 | 9 | NaN | +| 27 | 1 | 2 | 10 | Some data | +| 28 | 1 | 2 | 11 | NaN | +| 29 | 1 | 2 | 12 | NaN | +| 30 | 1 | 2 | 13 | NaN | +| 31 | 1 | 2 | 14 | NaN | +| 32 | 1 | 2 | 15 | NaN | | ... | ... | ... | ... | ... |