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 53aff7c07..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 @@ -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 + 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 @@ -62,8 +63,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,& @@ -306,19 +308,19 @@ 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 - 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) @@ -506,7 +508,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) @@ -563,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 @@ -585,13 +587,12 @@ 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 + !rrXbo10Sep2025 cncol(:,:,:) = 0.0 + !rrXbo10Sep2025 cnpft(:,:,:,:) = 0.0 + n = 0 np = 0 @@ -603,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) @@ -614,19 +617,19 @@ 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 ) 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%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 ) @@ -644,7 +647,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 ) @@ -786,7 +789,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/CNCLM_init_mod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/CNCLM_init_mod.F90 index 0ae6280e0..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, 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 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..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 @@ -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, .true.) ! valid(v)) !rrXbo 10Sep2025 + 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..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 @@ -500,7 +499,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, .true.) ! valid(v)) !rrXbo 10Sep2025 + end function safe !--------------------------------------- subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft, carbon_type, cn5_cold_start, rc) @@ -1126,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 @@ -1143,8 +1156,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..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 @@ -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, .true.) ! valid(v)) !rrXbo 10Sep2025 + 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. !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 @@ -469,13 +484,14 @@ 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 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) @@ -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 + !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) 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 6354b3163..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 @@ -4,8 +4,9 @@ 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 !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 @@ -13,7 +14,7 @@ module CNVegNitrogenStateType use decompMod , only : bounds_type use pftconMod , only : npcropmin use PatchType , only : patch - + ! !PUBLIC TYPES: implicit none save @@ -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, .true.) ! valid(v)) !rrXbo 10Sep2025 + end function safe !------------------------------------------------------------- subroutine Init(this, bounds, nch, ityp, fveg, cncol, cnpft) @@ -445,35 +461,42 @@ 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) ) + + 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 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) + 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) ) + 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 +614,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..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,13 +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 @@ -223,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 b9b4c616c..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 @@ -4,12 +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 @@ -101,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 @@ -154,7 +158,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) @@ -167,7 +171,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/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/PhotosynthesisMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM51_GridComp/CLM51/PhotosynthesisMod.F90 index b1cc0f0c8..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 @@ -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,27 +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 - 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 - 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 @@ -756,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 @@ -1967,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 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 4d1f3057e..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,10 +5,11 @@ 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 + 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 @@ -23,6 +24,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 +75,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 +90,32 @@ 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, .true.) ! valid(v)) !rrXbo 10Sep2025 + end function safe + +!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) @@ -99,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. !----------------------------------- @@ -162,31 +192,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) = 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 - 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) = 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) = cncol(nc,nz,36); - this%smin_nh4_col(n) = cncol(nc,nz,37); + 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) 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) = 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) = 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 +242,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/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 94531140d..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 @@ -3,17 +3,17 @@ 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 ! !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,32 +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, 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: + 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: 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 - 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 @@ -145,31 +150,35 @@ subroutine Init(this, bounds, nch, 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 - 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 - 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 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..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 @@ -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 @@ -39,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 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 18a8a58eb..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 @@ -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 @@ -133,6 +134,7 @@ module GEOS_CatchCNCLM51GridCompMod ! 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 +147,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: ! @@ -1874,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 ,& @@ -5532,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 @@ -5994,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 @@ -6233,9 +6233,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. @@ -7124,73 +7124,94 @@ 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) - ! --------------------------------------------------------------------------------- - if(init_accum) then - - ! (1) 5-day exponential moving average of snow depth - accper = min(istep,n5d) - SNDZM5D = ((accper-1)*SNDZM5D + SNDZM) / accper + ! --------------------------------------------------------------------------------- + + ! 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 - ! (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 + 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 - ! (2) 30-day exponential moving average of relative humidity [%] - accper = min(istep,n30d) - RH30D = ((accper-1)*RH30D + Qair_relative) / accper + ! create mask for tiles within longitude band (4:30-5:30am local time); + ! assumes -pi <= LONS <= pi + if (tmplon11) ) 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) + + 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 + + ! (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 @@ -7356,7 +7377,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) @@ -8682,8 +8703,8 @@ subroutine Driver ( RC ) deallocate( ht ) deallocate( tp ) deallocate( soilice ) - deallocate(TA_MIN) - deallocate(ta_count) + deallocate(mask_5amLT) + call MAPL_TimerOff ( MAPL, "-CATCHCNCLM51" ) RETURN_(ESMF_SUCCESS) @@ -9168,6 +9189,72 @@ subroutine RUN0(gc, import, export, clock, rc) end subroutine RUN0 + ! ***************************************************************** + + 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. + ! + ! - reichle, 12 Sep 2025 (simplified version of subroutine localtime2longitude() in GEOSldas_GridComp + + 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, intent(out) :: longitude ! [degree] -pi < longitude <= +pi + + integer, intent(out), optional :: rc ! return code for MAPL ASSERT() + + ! ---------------------------------------------- + + real :: UTC_hour, time_diff + + character(len=*), parameter :: Iam = 'LT2lon' + character(len=400) :: err_msg + + ! --------------------------------------------------------------------------- + + ! make sure local_hour is within range: 0 <= local_hour < 24 + + if ( (local_hour < 0.) .or. (local_hour >= 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. + + longitude = longitude*pi/180. ! [radians] + + end subroutine LT2lon + +! ***************************************************************** + end module GEOS_CatchCNCLM51GridCompMod subroutine SetServices(gc, rc) 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 | | ... | ... | ... | ... | ... | 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 =========================================================== 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 2fbc0dd7d..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 @@ -15,16 +15,15 @@ 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 /) + 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 @@ -98,7 +97,6 @@ module CatchmentCNRstMod real, allocatable :: runsurf(:) contains - procedure :: write_nc4 procedure :: allocate_cn procedure :: add_bcs_to_cnrst @@ -113,23 +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 :: 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 ) @@ -224,11 +232,63 @@ 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%RUNSURF , __RC__) ! rreichle, bug fix 9 Sep 2025: added for CLM51 - endif + endif - endif + !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() @@ -1114,8 +1174,9 @@ 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 + !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' @@ -1169,7 +1230,21 @@ 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)) @@ -1204,8 +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_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 @@ -1245,17 +1321,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 @@ -1281,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.) @@ -1311,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) @@ -1328,7 +1410,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.) @@ -1377,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.) @@ -1397,27 +1479,54 @@ 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) - + !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 + + ! 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 @@ -1557,6 +1666,12 @@ SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & end do OUT_TILE + !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 deallocate(this%cncol) allocate(this%cncol(NTILES, nzone*VAR_COL))