From 7c46935524876b2c7b00f270ff6fdd503c459950 Mon Sep 17 00:00:00 2001 From: Sean Swenson Date: Thu, 10 Jul 2025 09:25:50 -0600 Subject: [PATCH 1/6] initial implementation of spatially distributed parameters --- bld/CLMBuildNamelist.pm | 1 + bld/namelist_files/namelist_defaults_ctsm.xml | 3 + .../namelist_definition_ctsm.xml | 5 + src/biogeophys/CanopyHydrologyMod.F90 | 5 +- .../InfiltrationExcessRunoffMod.F90 | 7 +- src/biogeophys/SaturatedExcessRunoffMod.F90 | 9 +- ...nowCoverFractionSwensonLawrence2012Mod.F90 | 7 +- src/biogeophys/SnowHydrologyMod.F90 | 7 +- .../SoilHydrologyInitTimeConstMod.F90 | 5 +- src/biogeophys/SoilHydrologyMod.F90 | 42 +- src/biogeophys/SoilStateInitTimeConstMod.F90 | 29 +- src/biogeophys/SoilWaterMovementMod.F90 | 16 +- src/biogeophys/SurfaceResistanceMod.F90 | 19 +- src/biogeophys/SurfaceWaterMod.F90 | 5 +- src/biogeophys/WaterDiagnosticBulkType.F90 | 3 +- src/main/DistParamType.F90 | 947 ++++++++++++++++++ src/main/clm_initializeMod.F90 | 6 +- src/main/clm_varctl.F90 | 1 + src/main/controlMod.F90 | 4 +- src/main/initVerticalMod.F90 | 24 +- src/main/readParamsMod.F90 | 4 + 21 files changed, 1062 insertions(+), 87 deletions(-) create mode 100644 src/main/DistParamType.F90 diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 06ea82d99b..90fb47089f 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -2250,6 +2250,7 @@ sub setup_logic_params_file { add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'paramfile', 'phys'=>$nl_flags->{'phys'}, 'use_flexibleCN'=>$nl_flags->{'use_flexibleCN'} ); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'distributed_paramfile'); } #------------------------------------------------------------------------------- diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index 50e3cf68ad..91ec903dff 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -544,6 +544,9 @@ attributes from the config_cache.xml file (with keys converted to upper-case). lnd/clm2/paramdata/ctsm5.3.041.Nfix_params.v13.c250221_upplim250.nc lnd/clm2/paramdata/clm50_params.c250311.nc lnd/clm2/paramdata/clm45_params.c250311.nc + +lnd/clm2/paramdata/ctsm5.3_distributed_params_placeholder.nc +lnd/clm2/paramdata/ctsm5.3_distributed_params_placeholder.nc diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 820975655d..0b5d604697 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -998,6 +998,11 @@ Full pathname datafile with plant function type (PFT) constants combined with constants for biogeochem modules + +Full pathname spatially distributed parameter file + + Full pathname datafile with fates parameters diff --git a/src/biogeophys/CanopyHydrologyMod.F90 b/src/biogeophys/CanopyHydrologyMod.F90 index 166aa6d53d..a64c2e99fa 100644 --- a/src/biogeophys/CanopyHydrologyMod.F90 +++ b/src/biogeophys/CanopyHydrologyMod.F90 @@ -31,8 +31,9 @@ module CanopyHydrologyMod use WaterStateBulkType , only : waterstatebulk_type use WaterDiagnosticBulkType , only : waterdiagnosticbulk_type use WaterTracerUtils , only : CalcTracerFromBulk - use ColumnType , only : col, column_type - use PatchType , only : patch, patch_type + use ColumnType , only : col, column_type + use PatchType , only : patch, patch_type + ! ! !PUBLIC TYPES: implicit none diff --git a/src/biogeophys/InfiltrationExcessRunoffMod.F90 b/src/biogeophys/InfiltrationExcessRunoffMod.F90 index fff6a799ab..5fbe8141b0 100644 --- a/src/biogeophys/InfiltrationExcessRunoffMod.F90 +++ b/src/biogeophys/InfiltrationExcessRunoffMod.F90 @@ -14,9 +14,10 @@ module InfiltrationExcessRunoffMod use clm_varcon , only : spval use SoilHydrologyType, only : soilhydrology_type use SoilStateType , only : soilstate_type + use DistParamType , only : distparams use SaturatedExcessRunoffMod, only : saturated_excess_runoff_type - use WaterFluxBulkType , only : waterfluxbulk_type - use WaterDiagnosticBulkType, only : waterdiagnosticbulk_type + use WaterFluxBulkType , only : waterfluxbulk_type + use WaterDiagnosticBulkType , only : waterdiagnosticbulk_type implicit none save @@ -296,7 +297,7 @@ subroutine ComputeQinmaxHksat(bounds, num_hydrologyc, filter_hydrologyc, & do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) - qinmax_on_unsaturated_area(c) = minval(10._r8**(-params_inst%e_ice*(icefrac(c,1:3)))*hksat(c,1:3)) + qinmax_on_unsaturated_area(c) = minval(10._r8**(-distparams%e_ice(c)*(icefrac(c,1:3)))*hksat(c,1:3)) end do end associate diff --git a/src/biogeophys/SaturatedExcessRunoffMod.F90 b/src/biogeophys/SaturatedExcessRunoffMod.F90 index 9956a7dfb8..193fa6a628 100644 --- a/src/biogeophys/SaturatedExcessRunoffMod.F90 +++ b/src/biogeophys/SaturatedExcessRunoffMod.F90 @@ -16,9 +16,10 @@ module SaturatedExcessRunoffMod use clm_varcon , only : spval,ispval use LandunitType , only : landunit_type use landunit_varcon , only : istcrop - use ColumnType , only : column_type + use ColumnType , only : column_type + use DistParamType , only : distparams use SoilHydrologyType, only : soilhydrology_type - use SoilStateType, only : soilstate_type + use SoilStateType , only : soilstate_type use WaterFluxBulkType, only : waterfluxbulk_type implicit none @@ -349,9 +350,9 @@ subroutine ComputeFsatTopmodel(bounds, num_hydrologyc, filter_hydrologyc, & c = filter_hydrologyc(fc) if (frost_table(c) > zwt_perched(c) .and. frost_table(c) <= zwt(c)) then ! use perched water table to determine fsat (if present) - fsat(c) = wtfact(c) * exp(-0.5_r8*params_inst%fff*zwt_perched(c)) + fsat(c) = wtfact(c) * exp(-0.5_r8*distparams%fff(c)*zwt_perched(c)) else - fsat(c) = wtfact(c) * exp(-0.5_r8*params_inst%fff*zwt(c)) + fsat(c) = wtfact(c) * exp(-0.5_r8*distparams%fff(c)*zwt(c)) end if end do diff --git a/src/biogeophys/SnowCoverFractionSwensonLawrence2012Mod.F90 b/src/biogeophys/SnowCoverFractionSwensonLawrence2012Mod.F90 index f6a5ff41ee..35a8418b63 100644 --- a/src/biogeophys/SnowCoverFractionSwensonLawrence2012Mod.F90 +++ b/src/biogeophys/SnowCoverFractionSwensonLawrence2012Mod.F90 @@ -20,6 +20,7 @@ module SnowCoverFractionSwensonLawrence2012Mod use fileutils , only : getavu, relavu, opnfil use clm_varcon , only : rpi use ColumnType , only : column_type + use DistParamType , only : distparams use glcBehaviorMod , only : glc_behavior_type use landunit_varcon, only : istice use paramUtilMod , only : readNcdioScalar @@ -121,7 +122,7 @@ subroutine UpdateSnowDepthAndFrac(this, bounds, num_c, filter_c, & ! fsca parameterization based on *changes* in swe if (h2osno_total(c) == 0._r8) then if (newsnow(c) > 0._r8) then - frac_sno(c) = tanh(this%accum_factor * newsnow(c)) + frac_sno(c) = tanh(distparams%accum_factor(c) * newsnow(c)) else ! NOTE(wjs, 2019-08-07) This resetting of frac_sno to 0 when h2osno_total is 0 ! may already be done elsewhere; if it isn't, it possibly *should* be done @@ -146,7 +147,7 @@ subroutine UpdateSnowDepthAndFrac(this, bounds, num_c, filter_c, & ! ! This form is algebraically equivalent, but simpler and less prone to ! roundoff errors (see https://github.com/ESCOMP/ctsm/issues/784) - frac_sno(c) = frac_sno(c) + tanh(this%accum_factor * newsnow(c)) * (1._r8 - frac_sno(c)) + frac_sno(c) = frac_sno(c) + tanh(distparams%accum_factor(c) * newsnow(c)) * (1._r8 - frac_sno(c)) end if end if @@ -461,7 +462,7 @@ subroutine SetDerivedParameters(this, bounds, col, glc_behavior, n_melt_coef, n_ ! value of n_melt. this%n_melt(c) = n_melt_glcmec else - this%n_melt(c) = n_melt_coef / max(10._r8, col%topo_std(c)) + this%n_melt(c) = distparams%n_melt_coef(c) / max(10._r8, col%topo_std(c)) end if end do diff --git a/src/biogeophys/SnowHydrologyMod.F90 b/src/biogeophys/SnowHydrologyMod.F90 index 578769d9ea..19cc1ccce5 100644 --- a/src/biogeophys/SnowHydrologyMod.F90 +++ b/src/biogeophys/SnowHydrologyMod.F90 @@ -36,6 +36,7 @@ module SnowHydrologyMod use LandunitType , only : landunit_type, lun use TopoMod, only : topo_type use ColumnType , only : column_type, col + use DistParamType , only : distparams use landunit_varcon , only : istsoil, istdlak, istsoil, istwet, istice, istcrop use clm_time_manager, only : get_step_size_real, get_nstep use filterColMod , only : filter_col_type, col_filter_from_filter_and_logical_array @@ -70,7 +71,7 @@ module SnowHydrologyMod public :: SnowHydrologySetControlForTesting ! Set some of the control settings type, private :: params_type - real(r8) :: wimp ! Water impremeable if porosity less than wimp (unitless) + real(r8) :: wimp ! Water impermeable if porosity less than wimp (unitless) real(r8) :: ssi ! Irreducible water saturation of snow (unitless) real(r8) :: drift_gs ! Wind drift compaction / grain size (fixed value for now) (unitless) real(r8) :: eta0_anderson ! Viscosity coefficent from Anderson1976 (kg*s/m2) @@ -1972,8 +1973,8 @@ subroutine SnowCompaction(bounds, num_snowc, filter_snowc, & ! Settling as a result of destructive metamorphism ddz1 = -c3*dexpf - if (bi > params_inst%upplim_destruct_metamorph) ddz1 = & - ddz1*exp(-46.0e-3_r8*(bi-params_inst%upplim_destruct_metamorph)) + if (bi > distparams%upplim_destruct_metamorph(c)) ddz1 = & + ddz1*exp(-46.0e-3_r8*(bi-distparams%upplim_destruct_metamorph(c))) ! Liquid water term diff --git a/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 b/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 index 5f38830891..d95d886d01 100644 --- a/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 +++ b/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 @@ -13,6 +13,7 @@ module SoilHydrologyInitTimeConstMod use WaterStateBulkType, only : waterstatebulk_type use LandunitType , only : lun use ColumnType , only : col + use DistParamType , only : distparams use SoilStateInitTimeConstMod, only: organic_max ! implicit none @@ -187,7 +188,7 @@ subroutine SoilHydrologyInitTimeConst(bounds, soilhydrology_inst, soilstate_inst if ( lev <= nlevsoi )then claycol(c,lev) = soilstate_inst%cellclay_col(c,lev) sandcol(c,lev) = soilstate_inst%cellsand_col(c,lev) - om_fraccol(c,lev) = min(params_inst%om_frac_sf*soilstate_inst%cellorg_col(c,lev) / organic_max, 1._r8) + om_fraccol(c,lev) = min(distparams%om_frac_sf(c)*soilstate_inst%cellorg_col(c,lev) / organic_max, 1._r8) else claycol(c,lev) = soilstate_inst%cellclay_col(c,nlevsoi) sandcol(c,lev) = soilstate_inst%cellsand_col(c,nlevsoi) @@ -220,7 +221,7 @@ subroutine SoilHydrologyInitTimeConst(bounds, soilhydrology_inst, soilstate_inst if (micro_sigma(c) > 1.e-6_r8 .and. (soilhydrology_inst%h2osfcflag /= 0)) then d = 0.0_r8 do p = 1,4 - fd = 0.5_r8*(1.0_r8+shr_spfn_erf(d/(micro_sigma(c)*sqrt(2.0_r8)))) - params_inst%pc + fd = 0.5_r8*(1.0_r8+shr_spfn_erf(d/(micro_sigma(c)*sqrt(2.0_r8)))) - distparams%pc(c) dfdd = exp(-d**2/(2.0_r8*micro_sigma(c)**2))/(micro_sigma(c)*sqrt(2.0_r8*shr_const_pi)) d = d - fd/dfdd enddo diff --git a/src/biogeophys/SoilHydrologyMod.F90 b/src/biogeophys/SoilHydrologyMod.F90 index 5a4aa50f6e..1e6c90a2f4 100644 --- a/src/biogeophys/SoilHydrologyMod.F90 +++ b/src/biogeophys/SoilHydrologyMod.F90 @@ -10,7 +10,7 @@ module SoilHydrologyMod use abortutils , only : endrun use decompMod , only : bounds_type, subgrid_level_column use clm_varctl , only : iulog, use_vichydro - use clm_varcon , only : ispval + use clm_varcon , only : ispval, grlnd use clm_varcon , only : denh2o, denice, rpi use clm_varcon , only : pondmx_urban use clm_varpar , only : nlevsoi, nlevgrnd, nlayer, nlayert @@ -33,6 +33,7 @@ module SoilHydrologyMod use LandunitType , only : lun use ColumnType , only : column_type, col use PatchType , only : patch + use DistParamType , only : distparams ! ! !PUBLIC TYPES: @@ -791,7 +792,7 @@ subroutine WaterTable(bounds, num_hydrologyc, filter_hydrologyc, & ! use analytical expression for aquifer specific yield rous = watsat(c,nlevsoi) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,nlevsoi))**(-1./bsw(c,nlevsoi))) - rous=max(rous, params_inst%aq_sp_yield_min) + rous=max(rous, distparams%aq_sp_yield_min(c)) !-- water table is below the soil column -------------------------------------- if(jwt(c) == nlevsoi) then @@ -806,7 +807,7 @@ subroutine WaterTable(bounds, num_hydrologyc, filter_hydrologyc, & ! use analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) - s_y=max(s_y, params_inst%aq_sp_yield_min) + s_y=max(s_y, distparams%aq_sp_yield_min(c)) qcharge_layer=min(qcharge_tot,(s_y*(zwt(c) - zi(c,j-1))*1.e3)) qcharge_layer=max(qcharge_layer,0._r8) @@ -821,7 +822,7 @@ subroutine WaterTable(bounds, num_hydrologyc, filter_hydrologyc, & ! use analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) - s_y=max(s_y, params_inst%aq_sp_yield_min) + s_y=max(s_y, distparams%aq_sp_yield_min(c)) qcharge_layer=max(qcharge_tot,-(s_y*(zi(c,j) - zwt(c))*1.e3)) qcharge_layer=min(qcharge_layer,0._r8) @@ -1088,7 +1089,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte c = filter_hydrologyc(fc) ! specify maximum drainage rate - q_perch_max = params_inst%perched_baseflow_scalar * sin(col%topo_slope(c) * (rpi/180._r8)) + q_perch_max = distparams%perched_baseflow_scalar(c) * sin(col%topo_slope(c) * (rpi/180._r8)) ! if layer containing water table is frozen, compute the following: ! frost table, perched water table, and drainage from perched saturated layer @@ -1121,7 +1122,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte wtsub = 0._r8 q_perch = 0._r8 do k = jwt(c)+1, k_frz - imped=10._r8**(-params_inst%e_ice*(0.5_r8*(icefrac(c,k)+icefrac(c,min(nlevsoi, k+1))))) + imped=10._r8**(-distparams%e_ice(c)*(0.5_r8*(icefrac(c,k)+icefrac(c,min(nlevsoi, k+1))))) q_perch = q_perch + imped*hksat(c,k)*dzmm(c,k) wtsub = wtsub + dzmm(c,k) end do @@ -1197,7 +1198,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte wtsub = 0._r8 q_perch = 0._r8 do k = k_perch, k_frz - imped=10._r8**(-params_inst%e_ice*(0.5_r8*(icefrac(c,k)+icefrac(c,min(nlevsoi, k+1))))) + imped=10._r8**(-distparams%e_ice(c)*(0.5_r8*(icefrac(c,k)+icefrac(c,min(nlevsoi, k+1))))) q_perch = q_perch + imped*hksat(c,k)*dzmm(c,k) wtsub = wtsub + dzmm(c,k) end do @@ -1242,11 +1243,11 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte end do ! add ice impedance factor to baseflow if (use_vichydro) then - imped=10._r8**(-params_inst%e_ice*min(1.0_r8,ice(c,nlayer)/max_moist(c,nlayer))) + imped=10._r8**(-distparams%e_ice(c)*min(1.0_r8,ice(c,nlayer)/max_moist(c,nlayer))) dsmax_tmp(c) = Dsmax(c) * dtime/ secspday !mm/day->mm/dtime rsub_top_max = dsmax_tmp(c) else - imped=10._r8**(-params_inst%e_ice*(icefracsum/dzsum)) + imped=10._r8**(-distparams%e_ice(c)*(icefracsum/dzsum)) rsub_top_max = 10._r8 * sin((rpi/180.) * col%topo_slope(c)) end if @@ -1271,7 +1272,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte ! use analytical expression for aquifer specific yield rous = watsat(c,nlevsoi) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,nlevsoi))**(-1./bsw(c,nlevsoi))) - rous=max(rous, params_inst%aq_sp_yield_min) + rous=max(rous, distparams%aq_sp_yield_min(c)) !-- water table is below the soil column -------------------------------------- if(jwt(c) == nlevsoi) then @@ -1308,7 +1309,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte ! use analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) - s_y=max(s_y, params_inst%aq_sp_yield_min) + s_y=max(s_y, distparams%aq_sp_yield_min(c)) rsub_top_layer=max(rsub_top_tot,-(s_y*(zi(c,j) - zwt(c))*1.e3)) rsub_top_layer=min(rsub_top_layer,0._r8) @@ -1916,7 +1917,7 @@ subroutine PerchedLateralFlow(bounds, num_hydrologyc, & else ! Non-hillslope columns ! specify maximum drainage rate - q_perch_max = params_inst%perched_baseflow_scalar & + q_perch_max = distparams%perched_baseflow_scalar(c) & * sin(col%topo_slope(c) * (rpi/180._r8)) wtsub = 0._r8 @@ -1963,7 +1964,7 @@ subroutine PerchedLateralFlow(bounds, num_hydrologyc, & s_y = watsat(c,k) & * ( 1. - (1.+1.e3*zwt_perched(c)/sucsat(c,k))**(-1./bsw(c,k))) - s_y=max(s_y,params_inst%aq_sp_yield_min) + s_y=max(s_y,distparams%aq_sp_yield_min(c)) if (k==k_perch(c)) then drainage_layer=min(drainage_tot,(s_y*(zi(c,k) - zwt_perched(c))*1.e3)) else @@ -2177,6 +2178,7 @@ subroutine SubsurfaceLateralFlow(bounds, & icefrac => soilhydrology_inst%icefrac_col , & ! Output: [real(r8) (:,:) ] fraction of ice in layer frost_table => soilhydrology_inst%frost_table_col , & ! Input: [real(r8) (:) ] frost table depth (m) zwt => soilhydrology_inst%zwt_col , & ! Input: [real(r8) (:) ] water table depth (m) + stream_water_volume => waterstatebulk_inst%stream_water_volume_lun, & ! Input: [real(r8) (:) ] stream water volume (m3) qflx_snwcp_liq => waterfluxbulk_inst%qflx_snwcp_liq_col , & ! Output: [real(r8) (:) ] excess rainfall due to snow capping (mm H2O /s) [+] @@ -2201,7 +2203,7 @@ subroutine SubsurfaceLateralFlow(bounds, & vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice)) icefrac(c,j) = min(1._r8,vol_ice/watsat(c,j)) - ice_imped(c,j)=10._r8**(-params_inst%e_ice*icefrac(c,j)) + ice_imped(c,j)=10._r8**(-distparams%e_ice(c)*icefrac(c,j)) end do end do @@ -2243,7 +2245,7 @@ subroutine SubsurfaceLateralFlow(bounds, & dzsum = dzsum + dzmm(c,j) icefracsum = icefracsum + icefrac(c,j) * dzmm(c,j) end do - ice_imped_col(c)=10._r8**(-params_inst%e_ice*(icefracsum/dzsum)) + ice_imped_col(c)=10._r8**(-distparams%e_ice(c)*(icefracsum/dzsum)) enddo do fc = 1, num_hydrologyc @@ -2391,9 +2393,9 @@ subroutine SubsurfaceLateralFlow(bounds, & ! Non-hillslope columns ! baseflow is power law expression relative to bedrock layer if(zwt(c) <= zi(c,nbedrock(c))) then - qflx_latflow_out(c) = ice_imped_col(c) * baseflow_scalar & + qflx_latflow_out(c) = ice_imped_col(c) * distparams%baseflow_scalar(c) & * tan(rpi/180._r8*col%topo_slope(c))* & - (zi(c,nbedrock(c)) - zwt(c))**(params_inst%n_baseflow) + (zi(c,nbedrock(c)) - zwt(c))**(distparams%n_baseflow(c)) endif ! convert flux to volumetric flow qflx_latflow_out_vol(c) = 1.e-3_r8*qflx_latflow_out(c)*(grc%area(g)*1.e6_r8*col%wtgcell(c)) @@ -2456,7 +2458,7 @@ subroutine SubsurfaceLateralFlow(bounds, & ! analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) - s_y=max(s_y,params_inst%aq_sp_yield_min) + s_y=max(s_y,distparams%aq_sp_yield_min(c)) drainage_layer=min(drainage_tot,(s_y*dz(c,j)*1.e3)) @@ -2483,7 +2485,7 @@ subroutine SubsurfaceLateralFlow(bounds, & ! analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) - s_y=max(s_y,params_inst%aq_sp_yield_min) + s_y=max(s_y,distparams%aq_sp_yield_min(c)) drainage_layer=max(drainage_tot,-(s_y*(zi(c,j) - zwt(c))*1.e3)) drainage_layer=min(drainage_layer,0._r8) @@ -2822,7 +2824,7 @@ subroutine CalcIrrigWithdrawals(bounds, & ! use analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) - s_y=max(s_y, params_inst%aq_sp_yield_min) + s_y=max(s_y, distparams%aq_sp_yield_min(c)) if (j==jwt(c)+1) then available_water_layer=max(0._r8,(s_y*(zi(c,j) - zwt(c))*1.e3)) diff --git a/src/biogeophys/SoilStateInitTimeConstMod.F90 b/src/biogeophys/SoilStateInitTimeConstMod.F90 index a730417315..b75d3a5c27 100644 --- a/src/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/biogeophys/SoilStateInitTimeConstMod.F90 @@ -10,6 +10,7 @@ module SoilStateInitTimeConstMod use LandunitType , only : lun use ColumnType , only : col use PatchType , only : patch + use DistParamType , only : distparams use abortUtils , only : endrun use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) ! @@ -454,20 +455,20 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) if (lev .eq. 1) then clay = clay3d(g,1) sand = sand3d(g,1) - om_frac = min(params_inst%om_frac_sf*organic3d(g,1)/organic_max, 1._r8) + om_frac = min(distparams%om_frac_sf(c)*organic3d(g,1)/organic_max, 1._r8) else if (lev <= nlevsoi) then found = 0 ! reset value if (zsoi(lev) <= zisoifl(1)) then ! Search above the dataset's range of zisoifl depths clay = clay3d(g,1) sand = sand3d(g,1) - om_frac = min(params_inst%om_frac_sf*organic3d(g,1)/organic_max, 1._r8) + om_frac = min(distparams%om_frac_sf(c)*organic3d(g,1)/organic_max, 1._r8) found = 1 else if (zsoi(lev) > zisoifl(nlevsoifl)) then ! Search below the dataset's range of zisoifl depths clay = clay3d(g,nlevsoifl) sand = sand3d(g,nlevsoifl) - om_frac = min(params_inst%om_frac_sf*organic3d(g,nlevsoifl)/organic_max, 1._r8) + om_frac = min(distparams%om_frac_sf(c)*organic3d(g,nlevsoifl)/organic_max, 1._r8) found = 1 else ! For remaining model soil levels, search within dataset's @@ -477,7 +478,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) if (zsoi(lev) > zisoifl(j) .AND. zsoi(lev) <= zisoifl(j+1)) then clay = clay3d(g,j+1) sand = sand3d(g,j+1) - om_frac = min(params_inst%om_frac_sf*organic3d(g,j+1)/organic_max, 1._r8) + om_frac = min(distparams%om_frac_sf(c)*organic3d(g,j+1)/organic_max, 1._r8) found = 1 endif if (found == 1) exit ! no need to stay in the loop @@ -547,13 +548,13 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%bd_col(c,lev) = (1._r8 - soilstate_inst%watsat_col(c,lev))*params_inst%pd ! do not allow watsat_sf to push watsat above 0.93 - soilstate_inst%watsat_col(c,lev) = min(params_inst%watsat_sf * ( (1._r8 - om_frac) * & + soilstate_inst%watsat_col(c,lev) = min(distparams%watsat_sf(c) * ( (1._r8 - om_frac) * & soilstate_inst%watsat_col(c,lev) + om_watsat*om_frac ), 0.93_r8) tkm = (1._r8-om_frac) * (params_inst%tkd_sand*sand+params_inst%tkd_clay*clay)/ & (sand+clay)+params_inst%tkm_om*om_frac ! W/(m K) - soilstate_inst%bsw_col(c,lev) = params_inst%bsw_sf * ( (1._r8-om_frac) * & + soilstate_inst%bsw_col(c,lev) = distparams%bsw_sf(c) * ( (1._r8-om_frac) * & (2.91_r8 + 0.159_r8*clay) + om_frac*om_b ) - soilstate_inst%sucsat_col(c,lev) = params_inst%sucsat_sf * ( (1._r8-om_frac) * & + soilstate_inst%sucsat_col(c,lev) = distparams%sucsat_sf(c) * ( (1._r8-om_frac) * & soilstate_inst%sucsat_col(c,lev) + om_sucsat*om_frac ) soilstate_inst%hksat_min_col(c,lev) = xksat @@ -575,7 +576,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) else uncon_hksat = 0._r8 end if - soilstate_inst%hksat_col(c,lev) = params_inst%hksat_sf * ( uncon_frac*uncon_hksat + & + soilstate_inst%hksat_col(c,lev) = distparams%hksat_sf(c) * ( uncon_frac*uncon_hksat + & (perc_frac*om_frac)*om_hksat ) soilstate_inst%tkmg_col(c,lev) = tkm ** (1._r8- soilstate_inst%watsat_col(c,lev)) @@ -630,9 +631,9 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) clay = soilstate_inst%cellclay_col(c,lev) sand = soilstate_inst%cellsand_col(c,lev) if ( organic_frac_squared )then - om_frac = min( params_inst%om_frac_sf*(soilstate_inst%cellorg_col(c,lev)/organic_max)**2._r8, 1._r8) + om_frac = min( distparams%om_frac_sf(c)*(soilstate_inst%cellorg_col(c,lev)/organic_max)**2._r8, 1._r8) else - om_frac = min(params_inst%om_frac_sf*soilstate_inst%cellorg_col(c,lev)/organic_max, 1._r8) + om_frac = min(distparams%om_frac_sf(c)*soilstate_inst%cellorg_col(c,lev)/organic_max, 1._r8) end if else clay = soilstate_inst%cellclay_col(c,nlevsoi) @@ -649,16 +650,16 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) bd = (1._r8-soilstate_inst%watsat_col(c,lev))*params_inst%pd ! do not allow watsat_sf to push watsat above 0.93 - soilstate_inst%watsat_col(c,lev) = min(params_inst%watsat_sf * ( (1._r8 - om_frac) * & + soilstate_inst%watsat_col(c,lev) = min(distparams%watsat_sf(c) * ( (1._r8 - om_frac) * & soilstate_inst%watsat_col(c,lev) + om_watsat_lake * om_frac), 0.93_r8) tkm = (1._r8-om_frac)*(params_inst%tkd_sand*sand+params_inst%tkd_clay*clay)/(sand+clay) + & params_inst%tkm_om * om_frac ! W/(m K) - soilstate_inst%bsw_col(c,lev) = params_inst%bsw_sf * ( (1._r8-om_frac) * & + soilstate_inst%bsw_col(c,lev) = distparams%bsw_sf(c) * ( (1._r8-om_frac) * & (2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake ) - soilstate_inst%sucsat_col(c,lev) = params_inst%sucsat_sf * ( (1._r8-om_frac) * & + soilstate_inst%sucsat_col(c,lev) = distparams%sucsat_sf(c) * ( (1._r8-om_frac) * & soilstate_inst%sucsat_col(c,lev) + om_sucsat_lake * om_frac ) xksat = 0.0070556_r8 *( 10._r8**(-0.884_r8+0.0153_r8*sand) ) ! mm/s @@ -682,7 +683,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) uncon_hksat = 0._r8 end if - soilstate_inst%hksat_col(c,lev) = params_inst%hksat_sf * ( uncon_frac*uncon_hksat + & + soilstate_inst%hksat_col(c,lev) = distparams%hksat_sf(c) * ( uncon_frac*uncon_hksat + & (perc_frac*om_frac)*om_hksat_lake ) soilstate_inst%tkmg_col(c,lev) = tkm ** (1._r8- soilstate_inst%watsat_col(c,lev)) soilstate_inst%tksatu_col(c,lev) = soilstate_inst%tkmg_col(c,lev)*0.57_r8**soilstate_inst%watsat_col(c,lev) diff --git a/src/biogeophys/SoilWaterMovementMod.F90 b/src/biogeophys/SoilWaterMovementMod.F90 index 85bcf42c5e..3164de6b15 100644 --- a/src/biogeophys/SoilWaterMovementMod.F90 +++ b/src/biogeophys/SoilWaterMovementMod.F90 @@ -9,6 +9,7 @@ module SoilWaterMovementMod ! created by Jinyun Tang, Mar 12, 2014 use shr_kind_mod , only : r8 => shr_kind_r8 use shr_sys_mod , only : shr_sys_flush + use DistParamType , only : distparams ! implicit none @@ -723,7 +724,7 @@ subroutine soilwater_zengdecker2009(bounds, num_hydrologyc, filter_hydrologyc, & s1 = min(1._r8, s1) s2 = hksat(c,j)*s1**(2._r8*bsw(c,j)+2._r8) - imped(c,j)=10._r8**(-params_inst%e_ice*(0.5_r8*(icefrac(c,j)+icefrac(c,min(nlevsoi, j+1))))) + imped(c,j)=10._r8**(-distparams%e_ice(c)*(0.5_r8*(icefrac(c,j)+icefrac(c,min(nlevsoi, j+1))))) hk(c,j) = imped(c,j)*s1*s2 dhkdw(c,j) = imped(c,j)*(2._r8*bsw(c,j)+3._r8)*s2* & @@ -1504,10 +1505,10 @@ subroutine compute_hydraulic_properties(c, nlayers, & ! s1 is interface value, s2 is node value if(j==nlayers)then s1 = s2(j) - call IceImpedance(icefrac(c,j), imped(j) ) + call IceImpedance(c,icefrac(c,j), imped(j) ) else s1 = 0.5_r8 * (s2(j) + s2(j+1)) - call IceImpedance(0.5_r8*(icefrac(c,j) + icefrac(c,j+1)), imped(j) ) + call IceImpedance(c,0.5_r8*(icefrac(c,j) + icefrac(c,j+1)), imped(j) ) endif ! impose constraints on relative saturation at the layer interface @@ -2122,7 +2123,7 @@ end subroutine compute_qcharge !#13 !----------------------------------------------------------------------- - subroutine IceImpedance(icefrac, imped) + subroutine IceImpedance(c, icefrac, imped) ! !DESCRIPTION ! compute soil suction potential @@ -2133,15 +2134,16 @@ subroutine IceImpedance(icefrac, imped) ! ! !ARGUMENTS: implicit none - real(r8), intent(in) :: icefrac !fraction of pore space filled with ice + integer , intent(in) :: c ! column index + real(r8), intent(in) :: icefrac ! fraction of pore space filled with ice - real(r8), intent(out) :: imped !hydraulic conductivity reduction due to the presence of ice in pore space + real(r8), intent(out) :: imped ! hydraulic conductivity reduction due to the presence of ice in pore space ! ! !LOCAL VARIABLES: character(len=32) :: subname = 'IceImpedance' ! subroutine name !------------------------------------------------------------------------------ - imped = 10._r8**(-params_inst%e_ice*icefrac) + imped = 10._r8**(-distparams%e_ice(c)*icefrac) end subroutine IceImpedance diff --git a/src/biogeophys/SurfaceResistanceMod.F90 b/src/biogeophys/SurfaceResistanceMod.F90 index 9b04327252..d5fdbb5b9b 100644 --- a/src/biogeophys/SurfaceResistanceMod.F90 +++ b/src/biogeophys/SurfaceResistanceMod.F90 @@ -8,13 +8,14 @@ module SurfaceResistanceMod ! transported with BeTR. The surface here refers to water and soil, not including canopy ! ! !USES: - use shr_kind_mod , only: r8 => shr_kind_r8 - use shr_const_mod , only: SHR_CONST_TKFRZ - use clm_varctl , only: iulog - use SoilStateType , only: soilstate_type - use WaterStateBulkType, only: waterstatebulk_type - use WaterDiagnosticBulkType, only: waterdiagnosticbulk_type - use TemperatureType , only : temperature_type + use shr_kind_mod , only : r8 => shr_kind_r8 + use shr_const_mod , only : SHR_CONST_TKFRZ + use clm_varctl , only : iulog + use SoilStateType , only : soilstate_type + use DistParamType , only : distparams + use WaterStateBulkType, only : waterstatebulk_type + use WaterDiagnosticBulkType, only : waterdiagnosticbulk_type + use TemperatureType , only : temperature_type implicit none save private @@ -401,9 +402,9 @@ subroutine calc_soil_resistance_sl14(bounds, num_nolakec, filter_nolakec, & ! dsl(c) = dzmm(c,1)*max(0.001_r8,(0.8*eff_porosity(c,1) - vwc_liq)) & ! try arbitrary scaling (not top layer thickness) ! dsl(c) = 15._r8*max(0.001_r8,(0.8*eff_porosity(c,1) - vwc_liq)) & - dsl(c) = params_inst%d_max * max(0.001_r8, (params_inst%frac_sat_soil_dsl_init * eff_por_top - vwc_liq)) & + dsl(c) = distparams%d_max(c) * max(0.001_r8, (distparams%frac_sat_soil_dsl_init(c) * eff_por_top - vwc_liq)) & ! /max(0.001_r8,(watsat(c,1)- aird)) - / max(0.001_r8, (params_inst%frac_sat_soil_dsl_init * watsat(c,1) - aird)) + / max(0.001_r8, (distparams%frac_sat_soil_dsl_init(c) * watsat(c,1) - aird)) dsl(c)=max(dsl(c),0._r8) dsl(c)=min(dsl(c),200._r8) diff --git a/src/biogeophys/SurfaceWaterMod.F90 b/src/biogeophys/SurfaceWaterMod.F90 index 562c64cc18..8cd0d1c96c 100644 --- a/src/biogeophys/SurfaceWaterMod.F90 +++ b/src/biogeophys/SurfaceWaterMod.F90 @@ -15,6 +15,7 @@ module SurfaceWaterMod use column_varcon , only : icol_roof, icol_road_imperv, icol_sunwall, icol_shadewall, icol_road_perv use decompMod , only : bounds_type, subgrid_level_column use ColumnType , only : col + use DistParamType , only : distparams use NumericsMod , only : truncate_small_values use InfiltrationExcessRunoffMod , only : infiltration_excess_runoff_type use EnergyFluxType , only : energyflux_type @@ -473,10 +474,10 @@ subroutine QflxH2osfcSurf(bounds, num_hydrologyc, filter_hydrologyc, & c = filter_hydrologyc(fc) if (h2osfcflag==1) then - if (frac_h2osfc_nosnow(c) <= params_inst%pc) then + if (frac_h2osfc_nosnow(c) <= distparams%pc(c)) then frac_infclust=0.0_r8 else - frac_infclust=(frac_h2osfc_nosnow(c)-params_inst%pc)**params_inst%mu + frac_infclust=(frac_h2osfc_nosnow(c)-distparams%pc(c))**distparams%mu(c) endif endif diff --git a/src/biogeophys/WaterDiagnosticBulkType.F90 b/src/biogeophys/WaterDiagnosticBulkType.F90 index f91aaca761..bab5ab87fc 100644 --- a/src/biogeophys/WaterDiagnosticBulkType.F90 +++ b/src/biogeophys/WaterDiagnosticBulkType.F90 @@ -21,6 +21,7 @@ module WaterDiagnosticBulkType use clm_varcon , only : spval use LandunitType , only : lun use ColumnType , only : col + use DistParamType , only : distparams use filterColMod , only : filter_col_type, col_filter_from_ltypes use WaterDiagnosticType, only : waterdiagnostic_type use WaterInfoBaseType, only : water_info_base_type @@ -757,7 +758,7 @@ subroutine InitBulkCold(this, bounds, & fmelt = (snowbd/100.)**1. ! 100 is the assumed fresh snow density; 1 is a melting factor that could be ! reconsidered, optimal value of 1.5 in Niu et al., 2007 - this%frac_sno_col(c) = tanh( this%snow_depth_col(c) / (2.5 * params_inst%zlnd * fmelt) ) + this%frac_sno_col(c) = tanh( this%snow_depth_col(c) / (2.5 * distparams%zlnd(c) * fmelt) ) endif end if end do diff --git a/src/main/DistParamType.F90 b/src/main/DistParamType.F90 new file mode 100644 index 0000000000..c80a51959d --- /dev/null +++ b/src/main/DistParamType.F90 @@ -0,0 +1,947 @@ +module DistParamType + + !----------------------------------------------------------------------- + ! !DESCRIPTION: + ! Spatially distributed parameter data type allocation and initialization + ! -------------------------------------------------------- + ! + use shr_kind_mod , only : r8 => shr_kind_r8 + use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) + use shr_sys_mod , only : shr_sys_abort + use clm_varcon , only : spval, ispval, grlnd + use clm_varctl , only : iulog + use clm_varctl , only : paramfile, distributed_paramfile + use ColumnType , only : col + use spmdMod , only : masterproc, mpicom + use fileutils , only : getfil, opnfil, getavu, relavu + use shr_nl_mod , only : shr_nl_find_group_name + use shr_mpi_mod , only : shr_mpi_bcast + use decompMod , only : bounds_type + use ncdio_pio , only : ncd_pio_closefile, ncd_pio_openfile + use ncdio_pio , only : file_desc_t, ncd_inqdid, ncd_inqdlen + ! + ! !PUBLIC TYPES: + implicit none + save + private + character(len=*), parameter, private :: sourcefile = __FILE__ + ! + type, public :: distparam_type + ! CanopyHydrologyMod + real(r8), pointer :: liq_canopy_storage_scalar (:) ! Canopy-storage-of-liquid-water parameter (kg/m2) + real(r8), pointer :: snow_canopy_storage_scalar (:) ! Canopy-storage-of-snow parameter (kg/m2) + real(r8), pointer :: snowcan_unload_temp_fact (:) ! Temperature canopy snow unload scaling (C2 in Eq. 14, Roesch et al. 2001) (K*s) + real(r8), pointer :: snowcan_unload_wind_fact (:) ! Wind canopy snow unload scaling (modifies 1.56e5, where 1.56e5 is C3 in Eq. 15, Roesch et al. 2001) (-) + real(r8), pointer :: interception_fraction (:) ! Fraction of intercepted precipitation (-) + real(r8), pointer :: maximum_leaf_wetted_fraction (:) ! Maximum fraction of leaf that may be wet (-) + + ! SoilHydrologyMod + real(r8), pointer :: aq_sp_yield_min (:) ! Minimum aquifer specific yield (unitless) + real(r8), pointer :: n_baseflow (:) ! Drainage power law exponent (unitless) + real(r8), pointer :: perched_baseflow_scalar (:) ! Scalar multiplier for perched base flow rate (kg/m2/s) + real(r8), pointer :: e_ice (:) ! Soil ice impedance factor (unitless) + real(r8), pointer :: baseflow_scalar (:) ! Scalar multiplier for base flow rate () + + ! SaturatedExcessRunoff + real(r8), pointer :: fff (:) ! Decay factor for fractional saturated area (1/m) + + ! initVerticalMod + real(r8), pointer :: slopebeta (:) ! exponent for microtopography pdf sigma (unitless) + real(r8), pointer :: slopemax (:) ! max topographic slope for microtopography pdf sigma (unitless) + real(r8), pointer :: zbedrock_sf (:) ! parameter to scale zbedrock (m) + + ! SnowCoverFractionSwensonLawrence2012Mod + real(r8), pointer :: n_melt_coef (:) ! SCA shape parameter + real(r8), pointer :: accum_factor (:) ! Accumulation constant for fractional snow covered area (unitless) + + ! SnowHydrologyMod + real(r8), pointer :: upplim_destruct_metamorph (:) ! Upper limit on destructive metamorphism compaction (kg/m3) + + ! SoilHydrologyInitTimeConstMod + real(r8), pointer :: pc (:) ! Threshold probability for surface water (unitless) + real(r8), pointer :: om_frac_sf (:) ! Scale factor for organic matter fraction (unitless) + + ! SoilStateInitTimeConstMod + real(r8), pointer :: bsw_sf (:) ! Scale factor for bsw (unitless) + real(r8), pointer :: hksat_sf (:) ! Scale factor for hksat (unitless) + real(r8), pointer :: sucsat_sf (:) ! Scale factor for sucsat (unitless) + real(r8), pointer :: watsat_sf (:) ! Scale factor for watsat (unitless) + + ! SurfaceResistanceMod + real(r8), pointer :: d_max (:) ! Dry surface layer parameter (mm) + real(r8), pointer :: frac_sat_soil_dsl_init (:) ! Fraction of saturated soil for moisture value at which DSL initiates (unitless) + + ! SurfaceWaterMod + real(r8), pointer :: mu (:) ! Connectivity exponent for surface water (unitless) + + ! WaterDiagnosticBulkType + real(r8), pointer :: zlnd (:) ! Momentum roughness length for soil, glacier, wetland (m) + real(r8), pointer :: snw_rds_min (:) ! minimum allowed snow effective radius (also cold "fresh snow" value) [microns] + + ! atm2lndType + real(r8), pointer :: precip_repartition_nonglc_all_rain_t (:) ! Rain temperature threshold for non-glacier landunits (C) + real(r8), pointer :: precip_repartition_nonglc_all_snow_t (:) ! Snow temperature threshold for non-glacier landunits (C) + + contains + + procedure, public :: Init + procedure, public :: readDistributedParams + procedure, public :: Clean + + end type distparam_type + + type(distparam_type), public, target :: distparams ! + + !------------------------------------------------------------------------ + +contains + + !------------------------------------------------------------------------ + subroutine Init(this, bounds) + ! + ! !ARGUMENTS: + class(distparam_type) :: this + type(bounds_type), intent(in) :: bounds + ! + ! !LOCAL VARIABLES: + integer :: begc, endc + integer :: begg, endg + !------------------------------------------------------------------------ + + begc = bounds%begc; endc = bounds%endc + begg = bounds%begg; endg = bounds%endg + + ! these should be patch; todo + ! CanopyHydrology + allocate(this%liq_canopy_storage_scalar (begc:endc)) ; this%liq_canopy_storage_scalar (:) = nan + allocate(this%snow_canopy_storage_scalar (begc:endc)) ; this%snow_canopy_storage_scalar (:) = nan + allocate(this%snowcan_unload_temp_fact (begc:endc)) ; this%snowcan_unload_temp_fact (:) = nan + allocate(this%snowcan_unload_wind_fact (begc:endc)) ; this%snowcan_unload_wind_fact (:) = nan + allocate(this%interception_fraction (begc:endc)) ; this%interception_fraction (:) = nan + allocate(this%maximum_leaf_wetted_fraction (begc:endc)) ; this%maximum_leaf_wetted_fraction (:) = nan + + ! SoilHydrology + allocate(this%aq_sp_yield_min (begc:endc)) ; this%aq_sp_yield_min (:) = nan + allocate(this%n_baseflow (begc:endc)) ; this%n_baseflow (:) = nan + allocate(this%perched_baseflow_scalar (begc:endc)) ; this%perched_baseflow_scalar (:) = nan + allocate(this%e_ice (begc:endc)) ; this%e_ice (:) = nan + allocate(this%baseflow_scalar (begc:endc)) ; this%baseflow_scalar (:) = nan + + ! SaturatedExcessRunoff + allocate(this%fff (begc:endc)) ; this%fff (:) = nan + + ! initVertical + allocate(this%slopebeta (begg:endc)) ; this%slopebeta (:) = nan + allocate(this%slopemax (begg:endc)) ; this%slopemax (:) = nan + allocate(this%zbedrock_sf (begg:endg)) ; this%zbedrock_sf (:) = nan + + ! SnowCoverFractionSwensonLawrence2012Mod + allocate(this%n_melt_coef (begg:endc)) ; this%n_melt_coef (:) = nan + allocate(this%accum_factor (begg:endc)) ; this%accum_factor (:) = nan + + ! SnowHydrologyMod + allocate(this%upplim_destruct_metamorph (begg:endc)) ; this%upplim_destruct_metamorph (:) = nan + + ! SoilHydrologyInitTimeConstMod + allocate(this%pc (begg:endc)) ; this%pc (:) = nan + allocate(this%om_frac_sf (begg:endc)) ; this% om_frac_sf (:) = nan + + ! SoilStateInitTimeConstMod + allocate(this%bsw_sf (begg:endc)) ; this%bsw_sf (:) = nan + allocate(this%hksat_sf (begg:endc)) ; this%hksat_sf (:) = nan + allocate(this%sucsat_sf (begg:endc)) ; this%sucsat_sf (:) = nan + allocate(this%watsat_sf (begg:endc)) ; this%watsat_sf (:) = nan + + ! SurfaceResistanceMod + allocate(this%d_max (begg:endc)) ; this%d_max (:) = nan + allocate(this%frac_sat_soil_dsl_init (begg:endc)) ; this%frac_sat_soil_dsl_init (:) = nan + + ! SurfaceWaterMod + allocate(this%mu (begg:endc)) ; this%mu (:) = nan + + ! WaterDiagnosticBulkType + allocate(this%zlnd (begg:endc)) ; this%zlnd (:) = nan + allocate(this%snw_rds_min (begg:endc)) ; this%snw_rds_min (:) = nan + + ! atm2lndType + allocate(this%precip_repartition_nonglc_all_rain_t (begg:endc)) ; this%precip_repartition_nonglc_all_rain_t (:) = nan + allocate(this%precip_repartition_nonglc_all_snow_t (begg:endc)) ; this%precip_repartition_nonglc_all_snow_t (:) = nan + + ! allocate(this% (begg:endc)) ; this% (:) = nan + + end subroutine Init + + !----------------------------------------------------------------------- + subroutine readDistributedParams(this,bounds) + ! + ! !USES: + use ncdio_pio , only : check_var, ncd_io + use paramUtilMod , only : readNcdioScalar + use abortutils , only : endrun + use shr_log_mod , only : errMsg => shr_log_errMsg + use domainMod , only : ldomain + use clm_varctl , only : NLFilename => NLFilename_in + + ! + ! !ARGUMENTS: + implicit none + class(distparam_type) :: this + type(bounds_type), intent(in) :: bounds + ! + ! !LOCAL VARIABLES: + logical :: readvar ! whether the variable was found + integer :: c, g + integer :: dimid ! netCDF dimension id + integer :: dim1, dim2 ! dimension to compare + ! namelist variable names must match file + real(r8) :: baseflow_scalar ! read in namelist + real(r8) :: precip_repartition_nonglc_all_rain_t ! read in namelist + real(r8) :: precip_repartition_nonglc_all_snow_t ! read in namelist + real(r8) :: fscalar_in ! read in - scalar - float + real(r8), pointer :: fparam_in(:) ! read in - 1D - float + type(file_desc_t) :: ncids, ncidd ! pio netCDF file ids + character(len=256) :: locfns, locfnd ! local filenames + integer :: ierr ! error code + integer :: unitn ! unit for namelist file + character(len=256) :: nmlname ! namelist + character(len=*), parameter :: subname = 'readDistributedParams' + !-------------------------------------------------------------------- + + if (masterproc) then + write(iulog,*) trim(subname)//' :: reading CLM spatially distributed parameters' + end if + + ! global (scalar) parameters + call getfil (paramfile, locfns, 0) + call ncd_pio_openfile (ncids, trim(locfns), 0) + ! spatially distributed parameters + call getfil (distributed_paramfile, locfnd, 0) + call ncd_pio_openfile (ncidd, trim(locfnd), 0) + + ! should a check be added to make sure dimensions match surface data / mesh files? + call ncd_inqdid(ncidd,'lon',dimid) + call ncd_inqdlen(ncidd,dimid,dim1) + call ncd_inqdid(ncidd,'lat',dimid) + call ncd_inqdlen(ncidd,dimid,dim2) + if (dim1 /= ldomain%ni .or. dim2 /= ldomain%nj) then + call endrun(msg="ERROR dimensions of distributed parameter file incompatible with domain dimensions"//errmsg(sourcefile, __LINE__)) + endif + + !----------------------------------------------------------- + ! CanopyHydrology ! + !----------------------------------------------------------- + + ! these should be patch/pft not column... leaving in for now. + + ! Canopy-storage-of-liquid-water parameter (kg/m2) + call check_var(ncidd, 'liq_canopy_storage_scalar', readvar) + if (.false.) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%liq_canopy_storage_scalar(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'liq_canopy_storage_scalar', subname, fscalar_in) + this%liq_canopy_storage_scalar(:) = fscalar_in + endif + + ! Canopy-storage-of-snow parameter (kg/m2) + call check_var(ncidd, 'snow_canopy_storage_scalar', readvar) + if (.false.) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%snow_canopy_storage_scalar(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'snow_canopy_storage_scalar', subname, fscalar_in) + this%snow_canopy_storage_scalar(:) = fscalar_in + endif + + ! Temperature canopy snow unload scaling (C2 in Eq. 14, Roesch et al. 2001) (K*s) + call check_var(ncidd, 'snowcan_unload_temp_fact', readvar) + if (.false.) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='snowcan_unload_temp_fact', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%snowcan_unload_temp_fact(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'snowcan_unload_temp_fact', subname, fscalar_in) + this%snowcan_unload_temp_fact(:) = fscalar_in + endif + + ! Wind canopy snow unload scaling (modifies 1.56e5, where 1.56e5 is C3 in Eq. 15, Roesch et al. 2001) (-) + call check_var(ncidd, 'snowcan_unload_wind_fact', readvar) + if (.false.) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='snowcan_unload_wind_fact', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%snowcan_unload_wind_fact(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'snowcan_unload_wind_fact', subname, fscalar_in) + this%snowcan_unload_wind_fact(:) = fscalar_in + endif + + ! Fraction of intercepted precipitation (-) + call check_var(ncidd, 'interception_fraction', readvar) + if (.false.) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='interception_fraction', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%interception_fraction(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'interception_fraction', subname, fscalar_in) + this%interception_fraction(:) = fscalar_in + endif + + ! Maximum fraction of leaf that may be wet (-) + call check_var(ncidd, 'maximum_leaf_wetted_fraction', readvar) + if (.false.) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='maximum_leaf_wetted_fraction', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%maximum_leaf_wetted_fraction(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'maximum_leaf_wetted_fraction', subname, fscalar_in) + this%maximum_leaf_wetted_fraction(:) = fscalar_in + endif + + !----------------------------------------------------------- + ! SoilHydrology ! + !----------------------------------------------------------- + + ! Minimum aquifer specific yield (unitless) + call check_var(ncidd, 'aq_sp_yield_min', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='aq_sp_yield_min', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%aq_sp_yield_min(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'aq_sp_yield_min', subname, fscalar_in) + this%aq_sp_yield_min(:) = fscalar_in + endif + + ! Drainage power law exponent (unitless) + call check_var(ncidd, 'n_baseflow', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='n_baseflow', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%n_baseflow(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'n_baseflow', subname, fscalar_in) + this%n_baseflow(:) = fscalar_in + endif + + ! Scalar multiplier for perched base flow rate (kg/m2/s) + call check_var(ncidd, 'perched_baseflow_scalar', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='perched_baseflow_scalar', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%perched_baseflow_scalar(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'perched_baseflow_scalar', subname, fscalar_in) + this%perched_baseflow_scalar(:) = fscalar_in + endif + + ! Soil ice impedance factor (unitless) + call check_var(ncidd, 'e_ice', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='e_ice', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%e_ice(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'e_ice', subname, fscalar_in) + this%e_ice(:) = fscalar_in + endif + + ! Scalar multiplier for base flow rate () + call check_var(ncidd, 'baseflow_scalar', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='baseflow_scalar', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%baseflow_scalar(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! global parameter value from namelist + nmlname = 'soilhydrology_inparm' + namelist /soilhydrology_inparm/ baseflow_scalar + + if (masterproc) then + unitn = getavu() + write(iulog,*) 'Read in '//nmlname//' namelist' + call opnfil (NLFilename, unitn, 'F') + call shr_nl_find_group_name(unitn, nmlname, status=ierr) + if (ierr == 0) then + read(unitn, nml=soilhydrology_inparm, iostat=ierr) + if (ierr /= 0) then + call endrun(msg="ERROR reading "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) + end if + else + call endrun(msg="ERROR could NOT find "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) + end if + call relavu( unitn ) + end if + + call shr_mpi_bcast (baseflow_scalar, mpicom) + this%baseflow_scalar(:) = baseflow_scalar + endif + + !----------------------------------------------------------- + ! SaturatedExcessRunoff ! + !----------------------------------------------------------- + ! Decay factor for fractional saturated area (1/m) + call check_var(ncidd, 'fff', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='fff', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%fff(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'fff', subname, fscalar_in) + this%fff(:) = fscalar_in + endif + + + !----------------------------------------------------------- + ! initVertical ! + !----------------------------------------------------------- + + ! exponent for microtopography pdf sigma (unitless) + call check_var(ncidd, 'slopebeta', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='slopebeta', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%slopebeta(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'slopebeta', subname, fscalar_in) + this%slopebeta(:) = fscalar_in + endif + + ! max topographic slope for microtopography pdf sigma (unitless) + call check_var(ncidd, 'slopemax', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='slopemax', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%slopemax(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'slopemax', subname, fscalar_in) + this%slopemax(:) = fscalar_in + endif + + ! parameter to scale zbedrock (m) + call check_var(ncidd, 'zbedrock_sf', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='zbedrock_sf', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do g = bounds%begg,bounds%endg + this%zbedrock_sf(g) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'zbedrock_sf', subname, fscalar_in) + this%zbedrock_sf(:) = fscalar_in + endif + + !----------------------------------------------------------- + ! SnowCoverFractionSwensonLawrence2012 ! + !----------------------------------------------------------- + + ! SCA shape parameter + call check_var(ncidd, 'n_melt_coef', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='n_melt_coef', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%n_melt_coef(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'n_melt_coef', subname, fscalar_in) + this%n_melt_coef(:) = fscalar_in + endif + + ! Accumulation constant for fractional snow covered area (unitless) + call check_var(ncidd, 'accum_factor', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='accum_factor', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%accum_factor(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'accum_factor', subname, fscalar_in) + this%accum_factor(:) = fscalar_in + endif + + !----------------------------------------------------------- + ! SnowHydrologyMod ! + !----------------------------------------------------------- + + ! Upper limit on destructive metamorphism compaction (kg/m3) + call check_var(ncidd, 'upplim_destruct_metamorph', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='upplim_destruct_metamorph', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%upplim_destruct_metamorph(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'upplim_destruct_metamorph', subname, fscalar_in) + this%upplim_destruct_metamorph(:) = fscalar_in + endif + + !----------------------------------------------------------- + ! SoilHydrologyInitTimeConstMod ! + !----------------------------------------------------------- + + ! Scale factor for organic matter fraction (unitless) + call check_var(ncidd, 'om_frac_sf', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='om_frac_sf', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%om_frac_sf(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'om_frac_sf', subname, fscalar_in) + this%om_frac_sf(:) = fscalar_in + endif + + !----------------------------------------------------------- + ! SoilStateInitTimeConstMod ! + !----------------------------------------------------------- + + call check_var(ncidd, 'bsw_sf', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='bsw_sf', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%bsw_sf(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'bsw_sf', subname, fscalar_in) + this%bsw_sf(:) = fscalar_in + endif + + call check_var(ncidd, 'hksat_sf', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='hksat_sf', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%hksat_sf(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'hksat_sf', subname, fscalar_in) + this%hksat_sf(:) = fscalar_in + endif + + call check_var(ncidd, 'sucsat_sf', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='sucsat_sf', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%sucsat_sf(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'sucsat_sf', subname, fscalar_in) + this%sucsat_sf(:) = fscalar_in + endif + + call check_var(ncidd, 'watsat_sf', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='watsat_sf', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%watsat_sf(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'watsat_sf', subname, fscalar_in) + this%watsat_sf(:) = fscalar_in + endif + + !----------------------------------------------------------- + ! SurfaceResistanceMod ! + !----------------------------------------------------------- + + ! Dry surface layer parameter (mm) + call check_var(ncidd, 'd_max', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='d_max', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%d_max(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'd_max', subname, fscalar_in) + this%d_max(:) = fscalar_in + endif + + ! Fraction of saturated soil for moisture value at which DSL initiates (unitless) + call check_var(ncidd, 'frac_sat_soil_dsl_init', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='frac_sat_soil_dsl_init', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%frac_sat_soil_dsl_init(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'frac_sat_soil_dsl_init', subname, fscalar_in) + this%frac_sat_soil_dsl_init(:) = fscalar_in + endif + + !----------------------------------------------------------- + ! SurfaceWaterMod ! + !----------------------------------------------------------- + + + ! Threshold probability for surface water (unitless) + call check_var(ncidd, 'pc', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='pc', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%pc(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'pc', subname, fscalar_in) + this%pc(:) = fscalar_in + endif + + ! Connectivity exponent for surface water (unitless) + call check_var(ncidd, 'mu', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='mu', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%mu(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'mu', subname, fscalar_in) + this%mu(:) = fscalar_in + endif + + !----------------------------------------------------------- + ! WaterDiagnosticBulkType ! + !----------------------------------------------------------- + + ! Momentum roughness length for soil, glacier, wetland (m) + call check_var(ncidd, 'zlnd', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='zlnd', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%zlnd(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'zlnd', subname, fscalar_in) + this%zlnd(:) = fscalar_in + endif + + ! minimum allowed snow effective radius (also cold "fresh snow" value) [microns] + call check_var(ncidd, 'snw_rds_min', readvar) + if (.false.) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='snw_rds_min', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%snw_rds_min(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! use global parameter value + call readNcdioScalar(ncids, 'snw_rds_min', subname, fscalar_in) + this%snw_rds_min(:) = fscalar_in + endif + + !----------------------------------------------------------- + ! atm2lndType ! + !----------------------------------------------------------- + call check_var(ncidd, 'precip_repartition_nonglc_all_snow_t', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='precip_repartition_nonglc_all_snow_t', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%precip_repartition_nonglc_all_snow_t(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! global parameter value from namelist + nmlname = 'atm2lnd_inparm' + namelist /atm2lnd_inparm/ precip_repartition_nonglc_all_snow_t + + if (masterproc) then + unitn = getavu() + write(iulog,*) 'Read in '//nmlname//' namelist' + call opnfil (NLFilename, unitn, 'F') + call shr_nl_find_group_name(unitn, nmlname, status=ierr) + if (ierr == 0) then + read(unitn, nml=atm2lnd_inparm, iostat=ierr) + if (ierr /= 0) then + call endrun(msg="ERROR reading "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) + end if + else + call endrun(msg="ERROR could NOT find "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) + end if + call relavu( unitn ) + end if + + call shr_mpi_bcast (precip_repartition_nonglc_all_snow_t, mpicom) + this%precip_repartition_nonglc_all_snow_t(:) = precip_repartition_nonglc_all_snow_t + endif + + call check_var(ncidd, 'precip_repartition_nonglc_all_rain_t', readvar) + if (readvar) then + ! if distributed parameter found, overwrite global value + allocate(fparam_in(bounds%begg:bounds%endg)) + call ncd_io(ncid=ncidd, varname='precip_repartition_nonglc_all_rain_t', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) + + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + this%precip_repartition_nonglc_all_rain_t(c) = fparam_in(g) + enddo + + deallocate(fparam_in) + else ! global parameter value from namelist + nmlname = 'atm2lnd_inparm' + namelist /atm2lnd_inparm/ precip_repartition_nonglc_all_rain_t + + if (masterproc) then + unitn = getavu() + write(iulog,*) 'Read in '//nmlname//' namelist' + call opnfil (NLFilename, unitn, 'F') + call shr_nl_find_group_name(unitn, nmlname, status=ierr) + if (ierr == 0) then + read(unitn, nml=atm2lnd_inparm, iostat=ierr) + if (ierr /= 0) then + call endrun(msg="ERROR reading "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) + end if + else + call endrun(msg="ERROR could NOT find "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) + end if + call relavu( unitn ) + end if + + call shr_mpi_bcast (precip_repartition_nonglc_all_rain_t, mpicom) + this%precip_repartition_nonglc_all_rain_t(:) = precip_repartition_nonglc_all_rain_t + endif + + ! close parameter files + call ncd_pio_closefile(ncids) + call ncd_pio_closefile(ncidd) + + end subroutine readDistributedParams + + !------------------------------------------------------------------------ + subroutine Clean(this) + ! + ! !ARGUMENTS: + class(distparam_type) :: this + !------------------------------------------------------------------------ + + deallocate(this%liq_canopy_storage_scalar) + deallocate(this%snow_canopy_storage_scalar) + deallocate(this%snowcan_unload_temp_fact) + deallocate(this%snowcan_unload_wind_fact) + deallocate(this%interception_fraction) + deallocate(this%maximum_leaf_wetted_fraction) + + deallocate(this%aq_sp_yield_min) + deallocate(this%n_baseflow) + deallocate(this%perched_baseflow_scalar) + deallocate(this%e_ice) + + deallocate(this%fff) + + deallocate(this%slopebeta) + deallocate(this%slopemax) + deallocate(this%zbedrock_sf) + + deallocate(this%n_melt_coef) + deallocate(this%accum_factor) + + deallocate(this%upplim_destruct_metamorph) + + deallocate(this%pc) + deallocate(this%mu) + + deallocate(this%bsw_sf) + deallocate(this%hksat_sf) + deallocate(this%sucsat_sf) + deallocate(this%watsat_sf) + deallocate(this%om_frac_sf) + + deallocate(this%d_max) + deallocate(this%frac_sat_soil_dsl_init) + + deallocate(this%zlnd) + deallocate(this%snw_rds_min) + + deallocate(this%precip_repartition_nonglc_all_rain_t) + deallocate(this%precip_repartition_nonglc_all_snow_t) + + ! deallocate(this%) + + + end subroutine Clean + +end module Distparamtype diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90 index 8c0b50230b..a703bfeb00 100644 --- a/src/main/clm_initializeMod.F90 +++ b/src/main/clm_initializeMod.F90 @@ -25,6 +25,7 @@ module clm_initializeMod use LandunitType , only : lun ! instance use ColumnType , only : col ! instance use PatchType , only : patch ! instance + use DistParamType , only : distparams use reweightMod , only : reweight_wrapup use filterMod , only : allocFilters, filter, filter_inactive_and_active use CLMFatesInterfaceMod , only : CLMFatesGlobals1,CLMFatesGlobals2 @@ -342,12 +343,15 @@ subroutine initialize2(ni,nj, currtime) call get_proc_bounds(bounds_proc) nclumps = get_proc_clumps() - ! Read in parameters files + ! Read in parameters call clm_instReadNML( NLFilename ) allocate(nutrient_competition_method, & source=create_nutrient_competition_method(bounds_proc)) call readParameters(photosyns_inst) + ! Read in spatially distributed parameters + call distparams%Init(bounds_proc) + call distparams%readDistributedParams(bounds_proc) ! Initialize time manager if (nsrest == nsrStartup) then diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90 index 41978ae695..a4d7b59f87 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -111,6 +111,7 @@ module clm_varctl character(len=fname_len), public :: fsurdat = ' ' ! surface data file name character(len=fname_len), public :: hillslope_file = ' ' ! hillslope data file name character(len=fname_len), public :: paramfile = ' ' ! ASCII data file with PFT physiological constants + character(len=fname_len), public :: distributed_paramfile = ' ' ! spatially distributed parameter file name character(len=fname_len), public :: nrevsn = ' ' ! restart data file name for branch run character(len=fname_len), public :: fsnowoptics = ' ' ! snow optical properties file name character(len=fname_len), public :: fsnowaging = ' ' ! snow aging parameters file name diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90 index 6d363a9a6e..6e0be9055d 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -149,7 +149,7 @@ subroutine control_init(dtime) namelist /clm_inparm/ & fsurdat, hillslope_file, & - paramfile, fsnowoptics, fsnowaging + paramfile, distributed_paramfile, fsnowoptics, fsnowaging ! History, restart options @@ -755,6 +755,7 @@ subroutine control_spmd() call mpi_bcast (hillslope_file, len(hillslope_file), MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (fatmlndfrc,len(fatmlndfrc),MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (paramfile, len(paramfile) , MPI_CHARACTER, 0, mpicom, ier) + call mpi_bcast (distributed_paramfile, len(distributed_paramfile) , MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (fsnowoptics, len(fsnowoptics), MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (fsnowaging, len(fsnowaging), MPI_CHARACTER, 0, mpicom, ier) @@ -1045,6 +1046,7 @@ subroutine control_print () write(iulog,*) ' use_SSRE = ', use_SSRE write(iulog,*) 'input data files:' write(iulog,*) ' PFT physiology and parameters file = ',trim(paramfile) + write(iulog,*) ' distributed parameter file = ',trim(distributed_paramfile) if (fsurdat == ' ') then write(iulog,*) ' fsurdat, surface dataset not set' else diff --git a/src/main/initVerticalMod.F90 b/src/main/initVerticalMod.F90 index 64383e7a7c..666e421d06 100644 --- a/src/main/initVerticalMod.F90 +++ b/src/main/initVerticalMod.F90 @@ -27,6 +27,7 @@ module initVerticalMod use LandunitType , only : lun use GridcellType , only : grc use ColumnType , only : col + use DistParamType , only : distparams use glcBehaviorMod , only : glc_behavior_type use abortUtils , only : endrun use ncdio_pio @@ -45,10 +46,7 @@ module initVerticalMod ! !PRIVATE MEMBER FUNCTIONS: private :: hasBedrock ! true if the given column type includes bedrock layers type, private :: params_type - real(r8) :: slopebeta ! exponent for microtopography pdf sigma (unitless) - real(r8) :: slopemax ! max topographic slope for microtopography pdf sigma (unitless) real(r8) :: zbedrock ! parameter to substitute for zbedrock (m) - real(r8) :: zbedrock_sf ! parameter to scale zbedrock (m) end type params_type type(params_type), private :: params_inst ! @@ -76,13 +74,7 @@ subroutine readParams( ncid ) character(len=*), parameter :: subname = 'readParams_initVertical' !-------------------------------------------------------------------- - ! Exponent for microtopography pdf sigma (unitless) - call readNcdioScalar(ncid, 'slopebeta', subname, params_inst%slopebeta) - ! Max topographic slope for microtopography pdf sigma (unitless) - call readNcdioScalar(ncid, 'slopemax', subname, params_inst%slopemax) - call readNcdioScalar(ncid, 'zbedrock', subname, params_inst%zbedrock) - call readNcdioScalar(ncid, 'zbedrock_sf', subname, params_inst%zbedrock_sf) end subroutine readParams @@ -455,9 +447,11 @@ subroutine initVertical(bounds, glc_behavior, thick_wall, thick_roof) if (params_inst%zbedrock>=0._r8) then zbedrock_in(:) = params_inst%zbedrock end if - if (params_inst%zbedrock_sf/=1._r8) then - zbedrock_in(:) = params_inst%zbedrock_sf*zbedrock_in(:) - end if + do g = bounds%begg,bounds%endg + if (distparams%zbedrock_sf(g)/=1._r8) then + zbedrock_in(g) = distparams%zbedrock_sf(g)*zbedrock_in(g) + end if + enddo ! if use_bedrock = false, set zbedrock to lowest layer bottom interface else @@ -657,12 +651,12 @@ subroutine initVertical(bounds, glc_behavior, thick_wall, thick_roof) do c = begc,endc ! microtopographic parameter, units are meters (try smooth function of slope) - slope0 = params_inst%slopemax**(1._r8/params_inst%slopebeta) + slope0 = distparams%slopemax(c)**(1._r8/distparams%slopebeta(c)) if (col%is_hillslope_column(c)) then - col%micro_sigma(c) = (atan(col%hill_slope(c)) + slope0)**(params_inst%slopebeta) + col%micro_sigma(c) = (atan(col%hill_slope(c)) + slope0)**(distparams%slopebeta(c)) else - col%micro_sigma(c) = (col%topo_slope(c) + slope0)**(params_inst%slopebeta) + col%micro_sigma(c) = (col%topo_slope(c) + slope0)**(distparams%slopebeta(c)) endif end do diff --git a/src/main/readParamsMod.F90 b/src/main/readParamsMod.F90 index c11a741198..2467185822 100644 --- a/src/main/readParamsMod.F90 +++ b/src/main/readParamsMod.F90 @@ -8,14 +8,17 @@ module readParamsMod ! ! ! USES: use clm_varctl , only : paramfile, iulog, use_fates, use_cn + use clm_varctl , only : distributed_paramfile use SoilBiogeochemDecompCascadeConType, only : mimics_decomp, century_decomp, decomp_method use spmdMod , only : masterproc use fileutils , only : getfil use ncdio_pio , only : ncd_pio_closefile, ncd_pio_openfile use ncdio_pio , only : file_desc_t , ncd_inqdid, ncd_inqdlen + use decompMod , only : bounds_type implicit none private + character(len=*), parameter, private :: sourcefile = __FILE__ ! public :: readParameters @@ -63,6 +66,7 @@ subroutine readParameters (photosyns_inst) use SoilHydrologyInitTimeConstMod , only : readParams_SoilHydrologyInitTimeConst => readParams use clm_varctl, only : NLFilename_in use PhotosynthesisMod , only : photosyns_type + ! ! !ARGUMENTS: type(photosyns_type) , intent(in) :: photosyns_inst From 1fe77d5913f5bbda427a99668bd55137f378c7e9 Mon Sep 17 00:00:00 2001 From: Sean Swenson Date: Thu, 10 Jul 2025 09:50:08 -0600 Subject: [PATCH 2/6] revert zbedrock_sf and om_frac_sf --- .../SoilHydrologyInitTimeConstMod.F90 | 2 +- src/biogeophys/SoilStateInitTimeConstMod.F90 | 12 ++--- src/main/DistParamType.F90 | 53 +------------------ src/main/initVerticalMod.F90 | 10 ++-- 4 files changed, 14 insertions(+), 63 deletions(-) diff --git a/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 b/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 index d95d886d01..59b9fd4439 100644 --- a/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 +++ b/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 @@ -188,7 +188,7 @@ subroutine SoilHydrologyInitTimeConst(bounds, soilhydrology_inst, soilstate_inst if ( lev <= nlevsoi )then claycol(c,lev) = soilstate_inst%cellclay_col(c,lev) sandcol(c,lev) = soilstate_inst%cellsand_col(c,lev) - om_fraccol(c,lev) = min(distparams%om_frac_sf(c)*soilstate_inst%cellorg_col(c,lev) / organic_max, 1._r8) + om_fraccol(c,lev) = min(params_inst%om_frac_sf*soilstate_inst%cellorg_col(c,lev) / organic_max, 1._r8) else claycol(c,lev) = soilstate_inst%cellclay_col(c,nlevsoi) sandcol(c,lev) = soilstate_inst%cellsand_col(c,nlevsoi) diff --git a/src/biogeophys/SoilStateInitTimeConstMod.F90 b/src/biogeophys/SoilStateInitTimeConstMod.F90 index b75d3a5c27..96547d2c57 100644 --- a/src/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/biogeophys/SoilStateInitTimeConstMod.F90 @@ -455,20 +455,20 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) if (lev .eq. 1) then clay = clay3d(g,1) sand = sand3d(g,1) - om_frac = min(distparams%om_frac_sf(c)*organic3d(g,1)/organic_max, 1._r8) + om_frac = min(params_inst%om_frac_sf*organic3d(g,1)/organic_max, 1._r8) else if (lev <= nlevsoi) then found = 0 ! reset value if (zsoi(lev) <= zisoifl(1)) then ! Search above the dataset's range of zisoifl depths clay = clay3d(g,1) sand = sand3d(g,1) - om_frac = min(distparams%om_frac_sf(c)*organic3d(g,1)/organic_max, 1._r8) + om_frac = min(params_inst%om_frac_sf*organic3d(g,1)/organic_max, 1._r8) found = 1 else if (zsoi(lev) > zisoifl(nlevsoifl)) then ! Search below the dataset's range of zisoifl depths clay = clay3d(g,nlevsoifl) sand = sand3d(g,nlevsoifl) - om_frac = min(distparams%om_frac_sf(c)*organic3d(g,nlevsoifl)/organic_max, 1._r8) + om_frac = min(params_inst%om_frac_sf*organic3d(g,nlevsoifl)/organic_max, 1._r8) found = 1 else ! For remaining model soil levels, search within dataset's @@ -478,7 +478,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) if (zsoi(lev) > zisoifl(j) .AND. zsoi(lev) <= zisoifl(j+1)) then clay = clay3d(g,j+1) sand = sand3d(g,j+1) - om_frac = min(distparams%om_frac_sf(c)*organic3d(g,j+1)/organic_max, 1._r8) + om_frac = min(params_inst%om_frac_sf*organic3d(g,j+1)/organic_max, 1._r8) found = 1 endif if (found == 1) exit ! no need to stay in the loop @@ -631,9 +631,9 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) clay = soilstate_inst%cellclay_col(c,lev) sand = soilstate_inst%cellsand_col(c,lev) if ( organic_frac_squared )then - om_frac = min( distparams%om_frac_sf(c)*(soilstate_inst%cellorg_col(c,lev)/organic_max)**2._r8, 1._r8) + om_frac = min( params_inst%om_frac_sf*(soilstate_inst%cellorg_col(c,lev)/organic_max)**2._r8, 1._r8) else - om_frac = min(distparams%om_frac_sf(c)*soilstate_inst%cellorg_col(c,lev)/organic_max, 1._r8) + om_frac = min(params_inst%om_frac_sf*soilstate_inst%cellorg_col(c,lev)/organic_max, 1._r8) end if else clay = soilstate_inst%cellclay_col(c,nlevsoi) diff --git a/src/main/DistParamType.F90 b/src/main/DistParamType.F90 index c80a51959d..8a9cec4a56 100644 --- a/src/main/DistParamType.F90 +++ b/src/main/DistParamType.F90 @@ -48,7 +48,6 @@ module DistParamType ! initVerticalMod real(r8), pointer :: slopebeta (:) ! exponent for microtopography pdf sigma (unitless) real(r8), pointer :: slopemax (:) ! max topographic slope for microtopography pdf sigma (unitless) - real(r8), pointer :: zbedrock_sf (:) ! parameter to scale zbedrock (m) ! SnowCoverFractionSwensonLawrence2012Mod real(r8), pointer :: n_melt_coef (:) ! SCA shape parameter @@ -57,10 +56,6 @@ module DistParamType ! SnowHydrologyMod real(r8), pointer :: upplim_destruct_metamorph (:) ! Upper limit on destructive metamorphism compaction (kg/m3) - ! SoilHydrologyInitTimeConstMod - real(r8), pointer :: pc (:) ! Threshold probability for surface water (unitless) - real(r8), pointer :: om_frac_sf (:) ! Scale factor for organic matter fraction (unitless) - ! SoilStateInitTimeConstMod real(r8), pointer :: bsw_sf (:) ! Scale factor for bsw (unitless) real(r8), pointer :: hksat_sf (:) ! Scale factor for hksat (unitless) @@ -72,6 +67,7 @@ module DistParamType real(r8), pointer :: frac_sat_soil_dsl_init (:) ! Fraction of saturated soil for moisture value at which DSL initiates (unitless) ! SurfaceWaterMod + real(r8), pointer :: pc (:) ! Threshold probability for surface water (unitless) real(r8), pointer :: mu (:) ! Connectivity exponent for surface water (unitless) ! WaterDiagnosticBulkType @@ -133,7 +129,6 @@ subroutine Init(this, bounds) ! initVertical allocate(this%slopebeta (begg:endc)) ; this%slopebeta (:) = nan allocate(this%slopemax (begg:endc)) ; this%slopemax (:) = nan - allocate(this%zbedrock_sf (begg:endg)) ; this%zbedrock_sf (:) = nan ! SnowCoverFractionSwensonLawrence2012Mod allocate(this%n_melt_coef (begg:endc)) ; this%n_melt_coef (:) = nan @@ -142,10 +137,6 @@ subroutine Init(this, bounds) ! SnowHydrologyMod allocate(this%upplim_destruct_metamorph (begg:endc)) ; this%upplim_destruct_metamorph (:) = nan - ! SoilHydrologyInitTimeConstMod - allocate(this%pc (begg:endc)) ; this%pc (:) = nan - allocate(this%om_frac_sf (begg:endc)) ; this% om_frac_sf (:) = nan - ! SoilStateInitTimeConstMod allocate(this%bsw_sf (begg:endc)) ; this%bsw_sf (:) = nan allocate(this%hksat_sf (begg:endc)) ; this%hksat_sf (:) = nan @@ -157,6 +148,7 @@ subroutine Init(this, bounds) allocate(this%frac_sat_soil_dsl_init (begg:endc)) ; this%frac_sat_soil_dsl_init (:) = nan ! SurfaceWaterMod + allocate(this%pc (begg:endc)) ; this%pc (:) = nan allocate(this%mu (begg:endc)) ; this%mu (:) = nan ! WaterDiagnosticBulkType @@ -515,23 +507,6 @@ subroutine readDistributedParams(this,bounds) call readNcdioScalar(ncids, 'slopemax', subname, fscalar_in) this%slopemax(:) = fscalar_in endif - - ! parameter to scale zbedrock (m) - call check_var(ncidd, 'zbedrock_sf', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='zbedrock_sf', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do g = bounds%begg,bounds%endg - this%zbedrock_sf(g) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'zbedrock_sf', subname, fscalar_in) - this%zbedrock_sf(:) = fscalar_in - endif !----------------------------------------------------------- ! SnowCoverFractionSwensonLawrence2012 ! @@ -595,28 +570,6 @@ subroutine readDistributedParams(this,bounds) this%upplim_destruct_metamorph(:) = fscalar_in endif - !----------------------------------------------------------- - ! SoilHydrologyInitTimeConstMod ! - !----------------------------------------------------------- - - ! Scale factor for organic matter fraction (unitless) - call check_var(ncidd, 'om_frac_sf', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='om_frac_sf', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%om_frac_sf(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'om_frac_sf', subname, fscalar_in) - this%om_frac_sf(:) = fscalar_in - endif - !----------------------------------------------------------- ! SoilStateInitTimeConstMod ! !----------------------------------------------------------- @@ -914,7 +867,6 @@ subroutine Clean(this) deallocate(this%slopebeta) deallocate(this%slopemax) - deallocate(this%zbedrock_sf) deallocate(this%n_melt_coef) deallocate(this%accum_factor) @@ -928,7 +880,6 @@ subroutine Clean(this) deallocate(this%hksat_sf) deallocate(this%sucsat_sf) deallocate(this%watsat_sf) - deallocate(this%om_frac_sf) deallocate(this%d_max) deallocate(this%frac_sat_soil_dsl_init) diff --git a/src/main/initVerticalMod.F90 b/src/main/initVerticalMod.F90 index 666e421d06..dd67ff8169 100644 --- a/src/main/initVerticalMod.F90 +++ b/src/main/initVerticalMod.F90 @@ -47,6 +47,7 @@ module initVerticalMod private :: hasBedrock ! true if the given column type includes bedrock layers type, private :: params_type real(r8) :: zbedrock ! parameter to substitute for zbedrock (m) + real(r8) :: zbedrock_sf ! parameter to scale zbedrock (m) end type params_type type(params_type), private :: params_inst ! @@ -75,6 +76,7 @@ subroutine readParams( ncid ) !-------------------------------------------------------------------- call readNcdioScalar(ncid, 'zbedrock', subname, params_inst%zbedrock) + call readNcdioScalar(ncid, 'zbedrock_sf', subname, params_inst%zbedrock_sf) end subroutine readParams @@ -447,11 +449,9 @@ subroutine initVertical(bounds, glc_behavior, thick_wall, thick_roof) if (params_inst%zbedrock>=0._r8) then zbedrock_in(:) = params_inst%zbedrock end if - do g = bounds%begg,bounds%endg - if (distparams%zbedrock_sf(g)/=1._r8) then - zbedrock_in(g) = distparams%zbedrock_sf(g)*zbedrock_in(g) - end if - enddo + if (params_inst%zbedrock_sf/=1._r8) then + zbedrock_in(:) = params_inst%zbedrock_sf*zbedrock_in(:) + end if ! if use_bedrock = false, set zbedrock to lowest layer bottom interface else From f661443019df122791e7f91b02e9e3db8bcea478 Mon Sep 17 00:00:00 2001 From: Sean Swenson Date: Thu, 10 Jul 2025 10:30:49 -0600 Subject: [PATCH 3/6] remove some unused readparam routines --- src/biogeophys/SaturatedExcessRunoffMod.F90 | 26 ------------- ...nowCoverFractionSwensonLawrence2012Mod.F90 | 34 +---------------- src/biogeophys/SoilHydrologyMod.F90 | 37 +------------------ src/biogeophys/SurfaceResistanceMod.F90 | 29 --------------- src/biogeophys/SurfaceWaterMod.F90 | 28 -------------- src/main/readParamsMod.F90 | 8 ---- 6 files changed, 2 insertions(+), 160 deletions(-) diff --git a/src/biogeophys/SaturatedExcessRunoffMod.F90 b/src/biogeophys/SaturatedExcessRunoffMod.F90 index 193fa6a628..7b7af15986 100644 --- a/src/biogeophys/SaturatedExcessRunoffMod.F90 +++ b/src/biogeophys/SaturatedExcessRunoffMod.F90 @@ -51,12 +51,6 @@ module SaturatedExcessRunoffMod procedure, private, nopass :: ComputeFsatTopmodel procedure, private, nopass :: ComputeFsatVic end type saturated_excess_runoff_type - public :: readParams - - type, private :: params_type - real(r8) :: fff ! Decay factor for fractional saturated area (1/m) - end type params_type - type(params_type), private :: params_inst ! !PRIVATE DATA MEMBERS: @@ -176,26 +170,6 @@ subroutine InitCold(this, bounds) end subroutine InitCold - !----------------------------------------------------------------------- - subroutine readParams( ncid ) - ! - ! !USES: - use ncdio_pio, only: file_desc_t - use paramUtilMod, only: readNcdioScalar - ! - ! !ARGUMENTS: - implicit none - type(file_desc_t),intent(inout) :: ncid ! pio netCDF file id - ! - ! !LOCAL VARIABLES: - character(len=*), parameter :: subname = 'readParams_SaturatedExcessRunoff' - !-------------------------------------------------------------------- - - ! Decay factor for fractional saturated area (1/m) - call readNcdioScalar(ncid, 'fff', subname, params_inst%fff) - - end subroutine readParams - ! ======================================================================== ! Science routines ! ======================================================================== diff --git a/src/biogeophys/SnowCoverFractionSwensonLawrence2012Mod.F90 b/src/biogeophys/SnowCoverFractionSwensonLawrence2012Mod.F90 index 35a8418b63..c007b906f2 100644 --- a/src/biogeophys/SnowCoverFractionSwensonLawrence2012Mod.F90 +++ b/src/biogeophys/SnowCoverFractionSwensonLawrence2012Mod.F90 @@ -45,7 +45,6 @@ module SnowCoverFractionSwensonLawrence2012Mod procedure, public :: Init procedure, private :: ReadNamelist - procedure, private :: ReadParams procedure, private :: CheckValidInputs procedure, private :: SetDerivedParameters end type snow_cover_fraction_swenson_lawrence_2012_type @@ -292,7 +291,6 @@ subroutine Init(this, bounds, col, glc_behavior, NLFilename, params_ncid) type(file_desc_t) , intent(inout) :: params_ncid ! pio netCDF file id for parameter file ! ! !LOCAL VARIABLES: - real(r8) :: n_melt_coef ! n_melt parameter (unitless) real(r8) :: n_melt_glcmec ! SCA shape parameter for glc_mec columns character(len=*), parameter :: subname = 'Init' @@ -302,10 +300,6 @@ subroutine Init(this, bounds, col, glc_behavior, NLFilename, params_ncid) NLFilename = NLFilename, & n_melt_glcmec = n_melt_glcmec) - call this%ReadParams( & - params_ncid = params_ncid, & - n_melt_coef = n_melt_coef) - if (masterproc) then call this%CheckValidInputs( & n_melt_glcmec = n_melt_glcmec) @@ -315,7 +309,6 @@ subroutine Init(this, bounds, col, glc_behavior, NLFilename, params_ncid) bounds = bounds, & col = col, & glc_behavior = glc_behavior, & - n_melt_coef = n_melt_coef, & n_melt_glcmec = n_melt_glcmec) end subroutine Init @@ -377,30 +370,6 @@ subroutine ReadNamelist(this, NLFilename, n_melt_glcmec) end subroutine ReadNamelist - !----------------------------------------------------------------------- - subroutine ReadParams(this, params_ncid, n_melt_coef) - ! - ! !DESCRIPTION: - ! Read netCDF parameters needed for the SwensonLawrence2012 method - ! - ! !ARGUMENTS: - class(snow_cover_fraction_swenson_lawrence_2012_type), intent(inout) :: this - type(file_desc_t) , intent(inout) :: params_ncid ! pio netCDF file id for parameter file - real(r8) , intent(out) :: n_melt_coef ! n_melt parameter (unitless) - ! - ! !LOCAL VARIABLES: - - character(len=*), parameter :: subname = 'ReadParams' - !----------------------------------------------------------------------- - - ! Accumulation constant for fractional snow covered area (unitless) - call readNcdioScalar(params_ncid, 'accum_factor', subname, this%accum_factor) - - ! n_melt parameter (unitless) - call readNcdioScalar(params_ncid, 'n_melt_coef', subname, n_melt_coef) - - end subroutine ReadParams - !----------------------------------------------------------------------- subroutine CheckValidInputs(this, n_melt_glcmec) ! @@ -431,7 +400,7 @@ subroutine CheckValidInputs(this, n_melt_glcmec) end subroutine CheckValidInputs !----------------------------------------------------------------------- - subroutine SetDerivedParameters(this, bounds, col, glc_behavior, n_melt_coef, n_melt_glcmec) + subroutine SetDerivedParameters(this, bounds, col, glc_behavior, n_melt_glcmec) ! ! !DESCRIPTION: ! Set parameters that are derived from other inputs @@ -441,7 +410,6 @@ subroutine SetDerivedParameters(this, bounds, col, glc_behavior, n_melt_coef, n_ type(bounds_type) , intent(in) :: bounds type(column_type) , intent(in) :: col type(glc_behavior_type) , intent(in) :: glc_behavior - real(r8) , intent(in) :: n_melt_coef real(r8) , intent(in) :: n_melt_glcmec ! SCA shape parameter for glc_mec columns ! ! !LOCAL VARIABLES: diff --git a/src/biogeophys/SoilHydrologyMod.F90 b/src/biogeophys/SoilHydrologyMod.F90 index 1e6c90a2f4..01e53b631c 100644 --- a/src/biogeophys/SoilHydrologyMod.F90 +++ b/src/biogeophys/SoilHydrologyMod.F90 @@ -58,16 +58,7 @@ module SoilHydrologyMod public :: RenewCondensation ! Misc. corrections public :: CalcIrrigWithdrawals ! Calculate irrigation withdrawals from groundwater by layer public :: WithdrawGroundwaterIrrigation ! Remove groundwater irrigation from unconfined and confined aquifers - public :: readParams - - type, private :: params_type - real(r8) :: aq_sp_yield_min ! Minimum aquifer specific yield (unitless) - real(r8) :: n_baseflow ! Drainage power law exponent (unitless) - real(r8) :: perched_baseflow_scalar ! Scalar multiplier for perched base flow rate (kg/m2/s) - real(r8) :: e_ice ! Soil ice impedance factor (unitless) - end type params_type - type(params_type), public :: params_inst - + !----------------------------------------------------------------------- real(r8), private :: baseflow_scalar = 1.e-2_r8 real(r8), parameter :: tolerance = 1.e-12_r8 ! tolerance for checking whether sublimation is greater than ice in top soil layer @@ -172,32 +163,6 @@ subroutine hillslope_hydrology_ReadNML(NLFilename) end subroutine hillslope_hydrology_ReadNML - !----------------------------------------------------------------------- - subroutine readParams( ncid ) - ! - ! !USES: - use ncdio_pio, only: file_desc_t - use paramUtilMod, only: readNcdioScalar - ! - ! !ARGUMENTS: - implicit none - type(file_desc_t),intent(inout) :: ncid ! pio netCDF file id - ! - ! !LOCAL VARIABLES: - character(len=*), parameter :: subname = 'readParams_SoilHydrology' - !-------------------------------------------------------------------- - - ! Minimum aquifer specific yield (unitless) - call readNcdioScalar(ncid, 'aq_sp_yield_min', subname, params_inst%aq_sp_yield_min) - ! Drainage power law exponent (unitless) - call readNcdioScalar(ncid, 'n_baseflow', subname, params_inst%n_baseflow) - ! Scalar multiplier for perched base flow rate (kg/m2/s) - call readNcdioScalar(ncid, 'perched_baseflow_scalar', subname, params_inst%perched_baseflow_scalar) - ! Soil ice impedance factor (unitless) - call readNcdioScalar(ncid, 'e_ice', subname, params_inst%e_ice) - - end subroutine readParams - !----------------------------------------------------------------------- subroutine soilHydReadNML( NLFilename ) ! diff --git a/src/biogeophys/SurfaceResistanceMod.F90 b/src/biogeophys/SurfaceResistanceMod.F90 index d5fdbb5b9b..bb54a040b9 100644 --- a/src/biogeophys/SurfaceResistanceMod.F90 +++ b/src/biogeophys/SurfaceResistanceMod.F90 @@ -30,13 +30,6 @@ module SurfaceResistanceMod public :: do_soilevap_beta, do_soil_resistance_sl14 ! public :: init_soil_resistance public :: soil_resistance_readNL - public :: readParams - - type, private :: params_type - real(r8) :: d_max ! Dry surface layer parameter (mm) - real(r8) :: frac_sat_soil_dsl_init ! Fraction of saturated soil for moisture value at which DSL initiates (unitless) - end type params_type - type(params_type), private :: params_inst character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -166,28 +159,6 @@ subroutine soil_resistance_readNL(NLFilename) endif end subroutine soil_resistance_readNL - - !------------------------------------------------------------------------------ - subroutine readParams( ncid ) - ! - ! !USES: - use ncdio_pio, only: file_desc_t - use paramUtilMod, only: readNcdioScalar - ! - ! !ARGUMENTS: - implicit none - type(file_desc_t),intent(inout) :: ncid ! pio netCDF file id - ! - ! !LOCAL VARIABLES: - character(len=*), parameter :: subname = 'readParams_SurfaceResistance' - !-------------------------------------------------------------------- - - ! Dry surface layer parameter (mm) - call readNcdioScalar(ncid, 'd_max', subname, params_inst%d_max) - ! Fraction of saturated soil for moisture value at which DSL initiates (unitless) - call readNcdioScalar(ncid, 'frac_sat_soil_dsl_init', subname, params_inst%frac_sat_soil_dsl_init) - - end subroutine readParams !------------------------------------------------------------------------------ subroutine calc_soilevap_resis(bounds, num_nolakec, filter_nolakec, & diff --git a/src/biogeophys/SurfaceWaterMod.F90 b/src/biogeophys/SurfaceWaterMod.F90 index 8cd0d1c96c..60ce1f67b6 100644 --- a/src/biogeophys/SurfaceWaterMod.F90 +++ b/src/biogeophys/SurfaceWaterMod.F90 @@ -33,45 +33,17 @@ module SurfaceWaterMod ! !PUBLIC MEMBER FUNCTIONS: public :: UpdateFracH2oSfc ! Determine fraction of land surfaces which are submerged public :: UpdateH2osfc ! Calculate fluxes out of h2osfc and update the h2osfc state - public :: readParams ! !PRIVATE MEMBER FUNCTIONS: private :: BulkDiag_FracH2oSfc ! Determine fraction of land surfaces which are submerged private :: QflxH2osfcSurf ! Compute qflx_h2osfc_surf private :: QflxH2osfcDrain ! Compute qflx_h2osfc_drain - type, private :: params_type - real(r8) :: pc ! Threshold probability for surface water (unitless) - real(r8) :: mu ! Connectivity exponent for surface water (unitless) - end type params_type - type(params_type), private :: params_inst character(len=*), parameter, private :: sourcefile = & __FILE__ contains - !----------------------------------------------------------------------- - subroutine readParams( ncid ) - ! - ! !USES: - use ncdio_pio, only: file_desc_t - use paramUtilMod, only: readNcdioScalar - ! - ! !ARGUMENTS: - implicit none - type(file_desc_t),intent(inout) :: ncid ! pio netCDF file id - ! - ! !LOCAL VARIABLES: - character(len=*), parameter :: subname = 'readParams_SurfaceWater' - !-------------------------------------------------------------------- - - ! Threshold probability for surface water (unitless) - call readNcdioScalar(ncid, 'pc', subname, params_inst%pc) - ! Connectivity exponent for surface water (unitless) - call readNcdioScalar(ncid, 'mu', subname, params_inst%mu) - - end subroutine readParams - !----------------------------------------------------------------------- subroutine UpdateFracH2oSfc(bounds, num_soilc, filter_soilc, & water_inst) diff --git a/src/main/readParamsMod.F90 b/src/main/readParamsMod.F90 index 2467185822..ffe0b6a123 100644 --- a/src/main/readParamsMod.F90 +++ b/src/main/readParamsMod.F90 @@ -52,17 +52,13 @@ subroutine readParameters (photosyns_inst) use CanopyFluxesMod , only : readParams_CanopyFluxes => readParams use UrbanFluxesMod , only : readParams_UrbanFluxes => readParams use CanopyHydrologyMod , only : readParams_CanopyHydrology => readParams - use SoilHydrologyMod , only : readParams_SoilHydrology => readParams use SoilStateInitTimeConstMod , only : readParams_SoilStateInitTimeConst => readParams use SoilWaterMovementMod , only : readParams_SoilWaterMovement => readParams - use SaturatedExcessRunoffMod , only : readParams_SaturatedExcessRunoff => readParams use InfiltrationExcessRunoffMod , only : readParams_InfiltrationExcessRunoff => readParams - use SurfaceResistanceMod , only : readParams_SurfaceResistance => readParams use WaterDiagnosticBulkType , only : readParams_WaterDiagnosticBulk => readParams use SnowHydrologyMod , only : readParams_SnowHydrology => readParams use SnowSnicarMod , only : readParams_SnowSnicar => readParams use initVerticalMod , only : readParams_initVertical => readParams - use SurfaceWaterMod , only : readParams_SurfaceWater => readParams use SoilHydrologyInitTimeConstMod , only : readParams_SoilHydrologyInitTimeConst => readParams use clm_varctl, only : NLFilename_in use PhotosynthesisMod , only : photosyns_type @@ -130,17 +126,13 @@ subroutine readParameters (photosyns_inst) call readParams_CanopyFluxes ( ncid ) call readParams_UrbanFluxes ( ncid ) call readParams_CanopyHydrology ( ncid ) - call readParams_SoilHydrology ( ncid ) call readParams_SoilStateInitTimeConst ( ncid ) - call readParams_SaturatedExcessRunoff ( ncid ) call readParams_SoilWaterMovement ( ncid ) call readParams_InfiltrationExcessRunoff ( ncid ) - call readParams_SurfaceResistance ( ncid ) call readParams_WaterDiagnosticBulk ( ncid ) call readParams_SnowHydrology ( ncid ) call readParams_SnowSnicar ( ncid ) call readParams_initVertical ( ncid ) - call readParams_SurfaceWater ( ncid ) call readParams_SoilHydrologyInitTimeConst ( ncid ) ! call ncd_pio_closefile(ncid) From c6d1fd763762596e04fc269146f638a00c2426eb Mon Sep 17 00:00:00 2001 From: Sean Swenson Date: Tue, 5 Aug 2025 11:46:05 -0600 Subject: [PATCH 4/6] use data streams --- src/cpl/share_esmf/DistParamStreamMod.F90 | 362 ++++++++++++++++++++++ src/main/DistParamMod.F90 | 66 ++++ 2 files changed, 428 insertions(+) create mode 100644 src/cpl/share_esmf/DistParamStreamMod.F90 create mode 100644 src/main/DistParamMod.F90 diff --git a/src/cpl/share_esmf/DistParamStreamMod.F90 b/src/cpl/share_esmf/DistParamStreamMod.F90 new file mode 100644 index 0000000000..83f7197c57 --- /dev/null +++ b/src/cpl/share_esmf/DistParamStreamMod.F90 @@ -0,0 +1,362 @@ +module DistParamsStreamMod + + + !----------------------------------------------------------------------- + ! !DESCRIPTION: + ! Contains methods for reading from a spatially distributed parameter file + + ! !USES + use ESMF + use dshr_strdata_mod , only : shr_strdata_type + use shr_kind_mod , only : r8 => shr_kind_r8, CL => shr_kind_cl + use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) + use shr_log_mod , only : errMsg => shr_log_errMsg + use spmdMod , only : mpicom, masterproc + use clm_varctl , only : iulog, FL => fname_len + use clm_varctl , only : use_distributed_parameters + use abortutils , only : endrun + use decompMod , only : bounds_type + use DistParamType , only : distributed_parameter_type, distparam_class + + ! !PUBLIC TYPES: + implicit none + private + + ! + type, public :: distributed_parameter_stream_type + ! + contains + + ! !PUBLIC MEMBER FUNCTIONS: + procedure, public :: Init ! Initialize and read data in + procedure, public :: Clean ! Clean and deallocate the object + + end type distributed_parameter_stream_type + + type(distributed_parameter_stream_type), public, target :: distributed_parameter_stream + + ! ! PRIVATE DATA: + type, private :: streamcontrol_type + character(len=FL) :: stream_fldFileName_distparams ! data Filename + character(len=FL) :: stream_meshfile_distparams ! mesh Filename + character(len=CL) :: mapalgo_distparams ! map algo + contains + procedure, private :: ReadNML ! Read in control namelist + end type streamcontrol_type + + type(streamcontrol_type), private :: control ! Stream control data + logical , private :: NMLRead = .false. ! If namelist has been read + logical , private :: InitDone = .false. ! If initialization of streams has been done + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + +!============================================================================== +contains +!============================================================================== + + !------------------------------------------------------------------------ + subroutine AssignDistributedParameter(this,bounds,dataptr1d,stream_varname) + ! + ! !USES: + ! + ! !ARGUMENTS: + implicit none + class(distparam_class), intent(inout) :: this ! distributed parameter + type(bounds_type) , intent(in) :: bounds + real(r8) , pointer, intent(in) :: dataptr1d(:) ! stream data + character(len=*) , intent(in) :: stream_varname ! name of parameter on stream file + ! + ! local variables + integer :: g, ig + character(len=*), parameter :: subname = 'AssignDistributedParameter' + !-------------------------------------------------------------------- + + if (trim(stream_varname) == trim(this%name)) then + allocate(this%val(bounds%begg:bounds%endg)) + this%val(:) = nan + ig = 0 + do g = bounds%begg,bounds%endg + ig = ig+1 + this%val(g) = dataptr1d(ig) + end do + this%is_distributed = .true. + end if + + end subroutine AssignDistributedParameter + + !------------------------------------------------------------------------ + subroutine Init(this, bounds, NLFilename, distparams) + ! + ! Initialize the distributed parameters stream object + ! + ! Uses: + use spmdMod , only : iam + use lnd_comp_shr , only : mesh, model_clock + use dshr_strdata_mod , only : shr_strdata_init_from_inline, shr_strdata_print + use dshr_strdata_mod , only : shr_strdata_advance + use dshr_methods_mod , only : dshr_fldbun_getfldptr + use fileutils , only : getfil + use ncdio_pio , only : ncd_pio_closefile, ncd_pio_openfile + use ncdio_pio , only : file_desc_t, var_desc_t + use ncdio_pio , only : ncd_inqvname, ncd_inqnvars + ! + ! arguments + implicit none + class(distributed_parameter_stream_type) :: this + type(bounds_type), intent(in) :: bounds + character(len=*), intent(in) :: NLFilename ! Namelist filename + type(distributed_parameter_type), intent(inout) :: distparams + ! + ! local variables + integer :: ig, g, n ! Indices + integer :: begg, endg + integer :: year ! year (0, ...) for nstep+1 + integer :: mon ! month (1, ..., 12) for nstep+1 + integer :: day ! day of month (1, ..., 31) for nstep+1 + integer :: sec ! seconds into current date for nstep+1 + integer :: mcdate ! Current model date (yyyymmdd) + integer :: nparams ! number of parameters on file + integer :: nvariables ! number of variables on file + integer :: dimid ! netCDF dimension id + type(file_desc_t) :: ncid ! pio netCDF file id + character(len=256) :: locfn ! local filename + type(shr_strdata_type) :: sdat_distparams ! input data stream + character(len=256) :: varname ! variable name + type(var_desc_t) :: var_desc ! variable descriptor + character(len=CL), allocatable :: stream_varnames_in(:) ! array of stream field names + character(len=CL), allocatable :: stream_varnames(:) ! array of stream field names + integer :: rc ! error code + real(r8), pointer :: dataptr1d(:) ! temporary pointer + character(len=*), parameter :: stream_name = 'Distributed parameters' + character(len=*), parameter :: subname = 'DistParamsStream::Init' + !----------------------------------------------------------------------- + + call control%ReadNML( bounds, NLFileName ) + + if ( use_distributed_parameters )then + + ! query parameter file + call getfil (control%stream_fldFileName_distparams, locfn, 0) + call ncd_pio_openfile (ncid, trim(locfn), 0) + + call ncd_inqnvars(ncid,nvariables=nvariables) + + allocate(stream_varnames_in(nvariables)) + nparams = 0 + do n=1,nvariables + call ncd_inqvname(ncid,n,varname,var_desc) + if (.not. any(varname == (/'time','lat','lon'/))) then + nparams = nparams + 1 + write(stream_varnames_in(nparams),'(a)') trim(varname) + endif + enddo + + ! close parameter file + call ncd_pio_closefile(ncid) + + ! keep valid variable names (not entirely satisfying, but avoids having to explicitly index array later + allocate(stream_varnames(nparams)) + stream_varnames = stream_varnames_in(1:nparams) + deallocate(stream_varnames_in) + + if (masterproc) then + write(iulog,*) ' stream_varnames = ',stream_varnames + end if + + ! Initialize the cdeps data type sdat_distparams + call shr_strdata_init_from_inline(sdat_distparams, & + my_task = iam, & + logunit = iulog, & + compname = 'LND', & + model_clock = model_clock, & + model_mesh = mesh, & + stream_meshfile = control%stream_meshfile_distparams, & + stream_lev_dimname = 'null', & + stream_mapalgo = control%mapalgo_distparams, & + stream_filenames = (/trim(control%stream_fldFileName_distparams)/), & +! stream_fldlistFile = stream_varnames(:nparams), & +! stream_fldListModel = stream_varnames(:nparams), & + stream_fldlistFile = stream_varnames, & + stream_fldListModel = stream_varnames, & + stream_yearFirst = 2000, & + stream_yearLast = 2000, & + stream_yearAlign = 1, & + stream_offset = 0, & + stream_taxmode = 'extend', & + stream_dtlimit = 1.0e30_r8, & + stream_tintalgo = 'linear', & + stream_name = 'Distributed parameters', & + rc = rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + + ! Explicitly set current date to a hardcoded constant value. Otherwise + ! using the real date can cause roundoff differences that are + ! detrected as issues with exact restart. EBK M05/20/2017 + ! call get_curr_date(year, mon, day, sec) + year = 2000 + mon = 12 + day = 31 + sec = 0 + mcdate = year*10000 + mon*100 + day + + call shr_strdata_advance(sdat_distparams, ymd=mcdate, tod=sec, logunit=iulog, istr='distparams', rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + + ! Get pointer for stream data that is time and spatially interpolate to model time and grid + begg = bounds%begg; endg = bounds%endg + + do n = 1,size(stream_varnames) + call dshr_fldbun_getFldPtr(sdat_distparams%pstrm(1)%fldbun_model, stream_varnames(n), fldptr1=dataptr1d, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + + ! Assign values to distributed parameter arrays if stream data exist + call AssignDistributedParameter(distparams%aq_sp_yield_min,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%n_baseflow,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%perched_baseflow_scalar,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%e_ice,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%baseflow_scalar,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%fff,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%slopebeta,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%slopemax,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%n_melt_coef,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%accum_factor,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%upplim_destruct_metamorph,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%pc,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%mu,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%bsw_sf,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%hksat_sf,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%sucsat_sf,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%watsat_sf,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%d_max,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%frac_sat_soil_dsl_init,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%zlnd,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%snw_rds_min,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%precip_repartition_nonglc_all_rain_t,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%precip_repartition_nonglc_all_snow_t,bounds,dataptr1d,stream_varnames(n)) + + end do + deallocate(stream_varnames) + InitDone = .true. + end if + + end subroutine Init + + !============================================================================== + subroutine Clean(this) + ! + ! Deallocate and clean the object + ! + ! Uses: + ! + ! arguments + implicit none + class(distributed_parameter_stream_type) :: this + ! + ! local variables + !----------------------------------------------------------------------- + !deallocate(this%distparams) + !this%distparams => NULL() + InitDone = .false. + + end subroutine Clean + + !============================================================================== + subroutine ReadNML(this, bounds, NLFilename) + ! + ! Read the namelist data stream information. + ! + ! Uses: + use shr_nl_mod , only : shr_nl_find_group_name + use shr_log_mod , only : errMsg => shr_log_errMsg + use shr_mpi_mod , only : shr_mpi_bcast + ! + ! arguments + implicit none + class(streamcontrol_type) :: this + type(bounds_type), intent(in) :: bounds + character(len=*), intent(in) :: NLFilename ! Namelist filename + ! + ! local variables + integer :: nu_nml ! unit for namelist file + integer :: nml_error ! namelist i/o error flag + character(len=FL) :: stream_fldFileName_distparams = ' ' + character(len=FL) :: stream_meshfile_distparams = ' ' + character(len=CL) :: mapalgo_distparams = 'bilinear' + character(len=*), parameter :: namelist_name = 'distparams_streams' ! MUST agree with group name in namelist definition to read. + character(len=*), parameter :: subName = "('distparams_streams::ReadNML')" + !----------------------------------------------------------------------- + + namelist /distparams_streams/ & ! MUST agree with namelist_name above + use_distributed_parameters, stream_fldFileName_distparams, & + stream_meshfile_distparams + + ! Default values for namelist + + ! Read distparams_streams namelist + if (masterproc) then + open( newunit=nu_nml, file=trim(NLFilename), status='old', iostat=nml_error ) + call shr_nl_find_group_name(nu_nml, namelist_name, status=nml_error) + if (nml_error == 0) then + read(nu_nml, nml=distparams_streams,iostat=nml_error) ! MUST agree with namelist_name above + if (nml_error /= 0) then + call endrun(msg=' ERROR reading '//namelist_name//' namelist'//errMsg(sourcefile, __LINE__)) + end if + else + call endrun(msg=' ERROR finding '//namelist_name//' namelist'//errMsg(sourcefile, __LINE__)) + end if + close(nu_nml) + endif + + call shr_mpi_bcast(use_distributed_parameters , mpicom) + call shr_mpi_bcast(mapalgo_distparams , mpicom) + call shr_mpi_bcast(stream_fldFileName_distparams , mpicom) + call shr_mpi_bcast(stream_meshfile_distparams , mpicom) + + ! Error checking + if ( .not. use_distributed_parameters )then + if ( len_trim(stream_fldFileName_distparams) /= 0 )then + call endrun(msg=' ERROR stream_fldFileName_distparams is set, but use_distributed_parameters is FALSE' & + //errMsg(sourcefile, __LINE__)) + end if + if ( len_trim(stream_meshfile_distparams) /= 0 )then + call endrun(msg=' ERROR stream_meshfile_distparams is set, but use_distributed_parameters is FALSE' & + //errMsg(sourcefile, __LINE__)) + end if + else + if ( len_trim(stream_fldFileName_distparams) == 0 )then + call endrun(msg=' ERROR stream_fldFileName_distparams is NOT set, but use_distributed_parameters is TRUE' & + //errMsg(sourcefile, __LINE__)) + end if + if ( len_trim(stream_meshfile_distparams) == 0 )then + call endrun(msg=' ERROR stream_meshfile_distparams is NOT set, but use_distributed_parameters is TRUE' & + //errMsg(sourcefile, __LINE__)) + end if + end if + + if (masterproc) then + write(iulog,*) ' ' + write(iulog,*) namelist_name, ' stream settings:' + write(iulog,*) ' use_distributed_parameters = ',use_distributed_parameters + if ( use_distributed_parameters )then + write(iulog,*) ' stream_fldFileName_distparams = ',trim(stream_fldFileName_distparams) + write(iulog,*) ' stream_meshfile_distparams = ',trim(stream_meshfile_distparams) + write(iulog,*) ' mapalgo_distparams = ',trim(mapalgo_distparams) + end if + endif + this%stream_fldFileName_distparams = stream_fldFileName_distparams + this%stream_meshfile_distparams = stream_meshfile_distparams + this%mapalgo_distparams = mapalgo_distparams + + ! Mark namelist read as having been done + NMLRead = .true. + + end subroutine ReadNML + +end module DistParamsStreamMod diff --git a/src/main/DistParamMod.F90 b/src/main/DistParamMod.F90 new file mode 100644 index 0000000000..f3ac886475 --- /dev/null +++ b/src/main/DistParamMod.F90 @@ -0,0 +1,66 @@ +module DistParamMod + + !----------------------------------------------------------------------- + ! !DESCRIPTION: + ! Spatially distributed parameter data type allocation and initialization + ! -------------------------------------------------------- + ! + use shr_kind_mod , only : r8 => shr_kind_r8 + use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) + use shr_sys_mod , only : shr_sys_abort + use clm_varcon , only : spval, ispval, grlnd + use clm_varctl , only : iulog + use clm_varctl , only : use_distributed_parameters + use spmdMod , only : masterproc, mpicom + use shr_mpi_mod , only : shr_mpi_bcast + use decompMod , only : bounds_type + + ! + ! !PUBLIC TYPES: + implicit none + save + private + character(len=*), parameter, private :: sourcefile = __FILE__ + ! + + public :: InitDistributedParameters + + !------------------------------------------------------------------------ + +contains + + !------------------------------------------------------------------------ + subroutine InitDistributedParameters(bounds) + ! + ! !USES: + use DistParamType , only : distributed_parameters, InitGlobalParameters + use DistParamsStreamMod, only : distributed_parameter_stream + use clm_varctl , only : NLFilename => NLFilename_in + use abortutils , only : endrun + + ! + ! !ARGUMENTS: + implicit none + type(bounds_type), intent(in) :: bounds + ! + ! !LOCAL VARIABLES: + logical :: readvar ! whether the variable was found + character(len=*), parameter :: subname = 'InitDistributedParameters' + !-------------------------------------------------------------------- + + if (masterproc) then + write(iulog,*) trim(subname)//' :: reading CLM parameters' + end if + + ! Initialize distributed_parameters object + call distributed_parameters%Init() + + ! Initialize distributed parameters based on stream data + call distributed_parameter_stream%Init(bounds, NLFilename, distributed_parameters) + + ! Initialize global parameters + call InitGlobalParameters(bounds) + + end subroutine InitDistributedParameters + +end module DistParamMod From b654f54f62a5543733ba844f8f57d36bb33a02c7 Mon Sep 17 00:00:00 2001 From: Sean Swenson Date: Mon, 11 Aug 2025 08:15:21 -0600 Subject: [PATCH 5/6] refactor precip repartitioning --- bld/CLMBuildNamelist.pm | 39 +- bld/namelist_files/namelist_defaults_ctsm.xml | 21 +- .../namelist_definition_ctsm.xml | 44 +- .../InfiltrationExcessRunoffMod.F90 | 6 +- src/biogeophys/SaturatedExcessRunoffMod.F90 | 8 +- ...nowCoverFractionSwensonLawrence2012Mod.F90 | 10 +- src/biogeophys/SnowHydrologyMod.F90 | 6 +- .../SoilHydrologyInitTimeConstMod.F90 | 4 +- src/biogeophys/SoilHydrologyMod.F90 | 103 +- src/biogeophys/SoilStateInitTimeConstMod.F90 | 18 +- src/biogeophys/SoilWaterMovementMod.F90 | 8 +- src/biogeophys/SurfaceResistanceMod.F90 | 6 +- src/biogeophys/SurfaceWaterMod.F90 | 6 +- src/biogeophys/WaterDiagnosticBulkType.F90 | 4 +- src/cpl/share_esmf/DistParamStreamMod.F90 | 6 +- src/main/DistParamType.F90 | 977 +++++------------- src/main/atm2lndMod.F90 | 13 +- src/main/atm2lndType.F90 | 189 ++-- src/main/clm_initializeMod.F90 | 8 +- src/main/clm_varctl.F90 | 6 + src/main/controlMod.F90 | 3 +- src/main/initVerticalMod.F90 | 8 +- src/main/ncdio_pio.F90.in | 23 +- 23 files changed, 522 insertions(+), 994 deletions(-) diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 90fb47089f..ca79a93b04 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -1814,6 +1814,11 @@ sub process_namelist_inline_logic { ################################## setup_logic_cropcal_streams($opts, $nl_flags, $definition, $defaults, $nl); + ################################## + # namelist group: distparams_streams # + ################################## + setup_logic_distparams_streams($opts, $nl_flags, $definition, $defaults, $nl); + ########################################## # namelist group: soil_moisture_streams # ########################################## @@ -4421,6 +4426,25 @@ sub setup_logic_cropcal_streams { #------------------------------------------------------------------------------- +sub setup_logic_distparams_streams { + my ($opts, $nl_flags, $definition, $defaults, $nl) = @_; + + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'use_distributed_parameters'); + + # Add defaults if using any crop calendar input files + my $distparams_file = $nl->get_value('stream_fldFileName_distparams') ; + my $mesh_file = $nl->get_value('stream_meshfile_distparams') ; + + # User provided an input file but set mesh file to empty + if ( !&string_is_undef_or_empty($distparams_file)) { + if ( &string_is_undef_or_empty($mesh_file) ) { + $log->fatal_error("If providing a spatially distributed parameter file, you must provide stream_meshfile_distparams" ); + } + } +} + +#------------------------------------------------------------------------------- + sub setup_logic_soilwater_movement { my ($opts, $nl_flags, $definition, $defaults, $nl) = @_; @@ -4659,19 +4683,6 @@ sub setup_logic_atm_forcing { } } } - - foreach $var ("precip_repartition_glc_all_snow_t", - "precip_repartition_glc_all_rain_t", - "precip_repartition_nonglc_all_snow_t", - "precip_repartition_nonglc_all_rain_t") { - if ( &value_is_true($nl->get_value("repartition_rain_snow")) ) { - add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var); - } else { - if (defined($nl->get_value($var))) { - $log->fatal_error("$var can only be set if repartition_rain_snow is true"); - } - } - } } #------------------------------------------------------------------------------- @@ -5207,7 +5218,7 @@ sub write_output_files { @groups = qw(clm_inparm ndepdyn_nml popd_streams urbantv_streams light_streams soil_moisture_streams lai_streams atm2lnd_inparm lnd2atm_inparm clm_canopyhydrology_inparm cnphenology - cropcal_streams + cropcal_streams distparams_streams clm_soilhydrology_inparm dynamic_subgrid cnvegcarbonstate finidat_consistency_checks dynpft_consistency_checks clm_initinterp_inparm century_soilbgcdecompcascade diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index 91ec903dff..651835dd0a 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -227,20 +227,6 @@ attributes from the config_cache.xml file (with keys converted to upper-case). 0.5 -0. -2. - --2. -0. - -0. -2. - .true. .false. @@ -2085,6 +2071,13 @@ lnd/clm2/surfdata_esmf/NEON/ctsm5.3.0/surfdata_1x1_NEON_TOOL_hist_2000_78pfts_c2 lnd/clm2/cropdata/calendars/processed/hdates_ggcmi_crop_calendar_phase3_v1.01_nninterp-hcru_hcru_mt13.2000-2000.20230728_165845.tweaked_latlons.nc lnd/clm2/cropdata/calendars/processed/360x720_120830_ESMFmesh_c20210507_cdf5.tweaked_latlons.nc + +.false. + + +no_file_yet +no_file_yet + none diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 0b5d604697..227e8e10a7 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -1699,30 +1699,6 @@ This parameter must be in the range [0,1] Only relevant if glcmec_downscale_longwave is .true. - -Temperature below which all precipitation falls as snow, for glacier columns (deg C) -Only relevant if repartition_rain_snow is .true. - - - -Temperature above which all precipitation falls as rain, for glacier columns (deg C) -Only relevant if repartition_rain_snow is .true. - - - -Temperature below which all precipitation falls as snow, for non-glacier columns (deg C) -Only relevant if repartition_rain_snow is .true. - - - -Temperature above which all precipitation falls as rain, for non-glacier columns (deg C) -Only relevant if repartition_rain_snow is .true. - - @@ -2060,6 +2036,26 @@ Filename of input stream data for date (day of year) of end of gdd20 accumulatio Filename of input stream data for crop calendar inputs + + + + + + +Flag to enable spatially distributed parameters + + + +Filename of input stream data for spatially distributed parameters + + + +Filename of mesh data for spatially distributed parameters + + diff --git a/src/biogeophys/InfiltrationExcessRunoffMod.F90 b/src/biogeophys/InfiltrationExcessRunoffMod.F90 index 5fbe8141b0..d5e4acc3b4 100644 --- a/src/biogeophys/InfiltrationExcessRunoffMod.F90 +++ b/src/biogeophys/InfiltrationExcessRunoffMod.F90 @@ -14,7 +14,7 @@ module InfiltrationExcessRunoffMod use clm_varcon , only : spval use SoilHydrologyType, only : soilhydrology_type use SoilStateType , only : soilstate_type - use DistParamType , only : distparams + use DistParamType , only : distparams => distributed_parameters use SaturatedExcessRunoffMod, only : saturated_excess_runoff_type use WaterFluxBulkType , only : waterfluxbulk_type use WaterDiagnosticBulkType , only : waterdiagnosticbulk_type @@ -268,6 +268,8 @@ subroutine ComputeQinmaxHksat(bounds, num_hydrologyc, filter_hydrologyc, & soilhydrology_inst, soilstate_inst, & qinmax_on_unsaturated_area) ! + ! USES + use ColumnType , only : col ! !DESCRIPTION: ! Compute qinmax using a parameterization based on hksat ! @@ -297,7 +299,7 @@ subroutine ComputeQinmaxHksat(bounds, num_hydrologyc, filter_hydrologyc, & do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) - qinmax_on_unsaturated_area(c) = minval(10._r8**(-distparams%e_ice(c)*(icefrac(c,1:3)))*hksat(c,1:3)) + qinmax_on_unsaturated_area(c) = minval(10._r8**(-distparams%e_ice%param_val(col%gridcell(c))*(icefrac(c,1:3)))*hksat(c,1:3)) end do end associate diff --git a/src/biogeophys/SaturatedExcessRunoffMod.F90 b/src/biogeophys/SaturatedExcessRunoffMod.F90 index 7b7af15986..cdee4c7ec2 100644 --- a/src/biogeophys/SaturatedExcessRunoffMod.F90 +++ b/src/biogeophys/SaturatedExcessRunoffMod.F90 @@ -17,7 +17,7 @@ module SaturatedExcessRunoffMod use LandunitType , only : landunit_type use landunit_varcon , only : istcrop use ColumnType , only : column_type - use DistParamType , only : distparams + use DistParamType , only : distparams => distributed_parameters use SoilHydrologyType, only : soilhydrology_type use SoilStateType , only : soilstate_type use WaterFluxBulkType, only : waterfluxbulk_type @@ -290,6 +290,8 @@ end subroutine SaturatedExcessRunoff subroutine ComputeFsatTopmodel(bounds, num_hydrologyc, filter_hydrologyc, & soilhydrology_inst, soilstate_inst, fsat) ! + ! USES + use ColumnType , only : col ! !DESCRIPTION: ! Compute fsat using the TOPModel-based parameterization ! @@ -324,9 +326,9 @@ subroutine ComputeFsatTopmodel(bounds, num_hydrologyc, filter_hydrologyc, & c = filter_hydrologyc(fc) if (frost_table(c) > zwt_perched(c) .and. frost_table(c) <= zwt(c)) then ! use perched water table to determine fsat (if present) - fsat(c) = wtfact(c) * exp(-0.5_r8*distparams%fff(c)*zwt_perched(c)) + fsat(c) = wtfact(c) * exp(-0.5_r8*distparams%fff%param_val(col%gridcell(c))*zwt_perched(c)) else - fsat(c) = wtfact(c) * exp(-0.5_r8*distparams%fff(c)*zwt(c)) + fsat(c) = wtfact(c) * exp(-0.5_r8*distparams%fff%param_val(col%gridcell(c))*zwt(c)) end if end do diff --git a/src/biogeophys/SnowCoverFractionSwensonLawrence2012Mod.F90 b/src/biogeophys/SnowCoverFractionSwensonLawrence2012Mod.F90 index c007b906f2..a791b7ffe1 100644 --- a/src/biogeophys/SnowCoverFractionSwensonLawrence2012Mod.F90 +++ b/src/biogeophys/SnowCoverFractionSwensonLawrence2012Mod.F90 @@ -20,7 +20,7 @@ module SnowCoverFractionSwensonLawrence2012Mod use fileutils , only : getavu, relavu, opnfil use clm_varcon , only : rpi use ColumnType , only : column_type - use DistParamType , only : distparams + use DistParamType , only : distparams => distributed_parameters use glcBehaviorMod , only : glc_behavior_type use landunit_varcon, only : istice use paramUtilMod , only : readNcdioScalar @@ -66,6 +66,8 @@ subroutine UpdateSnowDepthAndFrac(this, bounds, num_c, filter_c, & lun_itype_col, urbpoi, h2osno_total, snowmelt, int_snow, newsnow, bifall, & snow_depth, frac_sno, frac_sno_eff) ! + ! USES + use ColumnType , only : col ! !DESCRIPTION: ! Update snow depth and snow fraction using the SwensonLawrence2012 parameterization ! @@ -121,7 +123,7 @@ subroutine UpdateSnowDepthAndFrac(this, bounds, num_c, filter_c, & ! fsca parameterization based on *changes* in swe if (h2osno_total(c) == 0._r8) then if (newsnow(c) > 0._r8) then - frac_sno(c) = tanh(distparams%accum_factor(c) * newsnow(c)) + frac_sno(c) = tanh(distparams%accum_factor%param_val(col%gridcell(c)) * newsnow(c)) else ! NOTE(wjs, 2019-08-07) This resetting of frac_sno to 0 when h2osno_total is 0 ! may already be done elsewhere; if it isn't, it possibly *should* be done @@ -146,7 +148,7 @@ subroutine UpdateSnowDepthAndFrac(this, bounds, num_c, filter_c, & ! ! This form is algebraically equivalent, but simpler and less prone to ! roundoff errors (see https://github.com/ESCOMP/ctsm/issues/784) - frac_sno(c) = frac_sno(c) + tanh(distparams%accum_factor(c) * newsnow(c)) * (1._r8 - frac_sno(c)) + frac_sno(c) = frac_sno(c) + tanh(distparams%accum_factor%param_val(col%gridcell(c)) * newsnow(c)) * (1._r8 - frac_sno(c)) end if end if @@ -430,7 +432,7 @@ subroutine SetDerivedParameters(this, bounds, col, glc_behavior, n_melt_glcmec) ! value of n_melt. this%n_melt(c) = n_melt_glcmec else - this%n_melt(c) = distparams%n_melt_coef(c) / max(10._r8, col%topo_std(c)) + this%n_melt(c) = distparams%n_melt_coef%param_val(col%gridcell(c)) / max(10._r8, col%topo_std(c)) end if end do diff --git a/src/biogeophys/SnowHydrologyMod.F90 b/src/biogeophys/SnowHydrologyMod.F90 index 19cc1ccce5..870d846e25 100644 --- a/src/biogeophys/SnowHydrologyMod.F90 +++ b/src/biogeophys/SnowHydrologyMod.F90 @@ -36,7 +36,7 @@ module SnowHydrologyMod use LandunitType , only : landunit_type, lun use TopoMod, only : topo_type use ColumnType , only : column_type, col - use DistParamType , only : distparams + use DistParamType , only : distparams => distributed_parameters use landunit_varcon , only : istsoil, istdlak, istsoil, istwet, istice, istcrop use clm_time_manager, only : get_step_size_real, get_nstep use filterColMod , only : filter_col_type, col_filter_from_filter_and_logical_array @@ -1973,8 +1973,8 @@ subroutine SnowCompaction(bounds, num_snowc, filter_snowc, & ! Settling as a result of destructive metamorphism ddz1 = -c3*dexpf - if (bi > distparams%upplim_destruct_metamorph(c)) ddz1 = & - ddz1*exp(-46.0e-3_r8*(bi-distparams%upplim_destruct_metamorph(c))) + if (bi > distparams%upplim_destruct_metamorph%param_val(col%gridcell(c))) ddz1 = & + ddz1*exp(-46.0e-3_r8*(bi-distparams%upplim_destruct_metamorph%param_val(col%gridcell(c)))) ! Liquid water term diff --git a/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 b/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 index 59b9fd4439..7289987fff 100644 --- a/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 +++ b/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 @@ -13,7 +13,7 @@ module SoilHydrologyInitTimeConstMod use WaterStateBulkType, only : waterstatebulk_type use LandunitType , only : lun use ColumnType , only : col - use DistParamType , only : distparams + use DistParamType , only : distparams => distributed_parameters use SoilStateInitTimeConstMod, only: organic_max ! implicit none @@ -221,7 +221,7 @@ subroutine SoilHydrologyInitTimeConst(bounds, soilhydrology_inst, soilstate_inst if (micro_sigma(c) > 1.e-6_r8 .and. (soilhydrology_inst%h2osfcflag /= 0)) then d = 0.0_r8 do p = 1,4 - fd = 0.5_r8*(1.0_r8+shr_spfn_erf(d/(micro_sigma(c)*sqrt(2.0_r8)))) - distparams%pc(c) + fd = 0.5_r8*(1.0_r8+shr_spfn_erf(d/(micro_sigma(c)*sqrt(2.0_r8)))) - distparams%pc%param_val(col%gridcell(c)) dfdd = exp(-d**2/(2.0_r8*micro_sigma(c)**2))/(micro_sigma(c)*sqrt(2.0_r8*shr_const_pi)) d = d - fd/dfdd enddo diff --git a/src/biogeophys/SoilHydrologyMod.F90 b/src/biogeophys/SoilHydrologyMod.F90 index 01e53b631c..a13d971977 100644 --- a/src/biogeophys/SoilHydrologyMod.F90 +++ b/src/biogeophys/SoilHydrologyMod.F90 @@ -33,7 +33,7 @@ module SoilHydrologyMod use LandunitType , only : lun use ColumnType , only : column_type, col use PatchType , only : patch - use DistParamType , only : distparams + use DistParamType , only : distparams => distributed_parameters ! ! !PUBLIC TYPES: @@ -41,7 +41,6 @@ module SoilHydrologyMod save ! ! !PUBLIC MEMBER FUNCTIONS: - public :: SoilHydReadNML ! Read in the Soil hydrology namelist public :: SetSoilWaterFractions ! Set diagnostic variables related to the fraction of water and ice in each layer public :: SetFloodc ! Apply gridcell flood water flux to non-lake columns public :: SetQflxInputs ! Set the flux of water into the soil from the top @@ -60,7 +59,6 @@ module SoilHydrologyMod public :: WithdrawGroundwaterIrrigation ! Remove groundwater irrigation from unconfined and confined aquifers !----------------------------------------------------------------------- - real(r8), private :: baseflow_scalar = 1.e-2_r8 real(r8), parameter :: tolerance = 1.e-12_r8 ! tolerance for checking whether sublimation is greater than ice in top soil layer integer, private :: head_gradient_method ! Method for calculating hillslope saturated head gradient @@ -77,7 +75,7 @@ module SoilHydrologyMod __FILE__ contains - + !----------------------------------------------------------------------- subroutine hillslope_hydrology_ReadNML(NLFilename) ! @@ -163,65 +161,6 @@ subroutine hillslope_hydrology_ReadNML(NLFilename) end subroutine hillslope_hydrology_ReadNML - !----------------------------------------------------------------------- - subroutine soilHydReadNML( NLFilename ) - ! - ! !DESCRIPTION: - ! Read the namelist for soil hydrology - ! - ! !USES: - use fileutils , only : getavu, relavu, opnfil - use shr_nl_mod , only : shr_nl_find_group_name - use spmdMod , only : masterproc, mpicom - use shr_mpi_mod , only : shr_mpi_bcast - use clm_varctl , only : iulog - use shr_log_mod , only : errMsg => shr_log_errMsg - ! - ! !ARGUMENTS: - character(len=*), intent(in) :: NLFilename ! Namelist filename - ! - ! !LOCAL VARIABLES: - integer :: ierr ! error code - integer :: unitn ! unit for namelist file - - character(len=*), parameter :: subname = 'soilHydReadNML' - character(len=*), parameter :: nmlname = 'soilhydrology_inparm' - !----------------------------------------------------------------------- - namelist /soilhydrology_inparm/ baseflow_scalar - - ! Initialize options to default values, in case they are not specified in - ! the namelist - - - if (masterproc) then - unitn = getavu() - write(iulog,*) 'Read in '//nmlname//' namelist' - call opnfil (NLFilename, unitn, 'F') - call shr_nl_find_group_name(unitn, nmlname, status=ierr) - if (ierr == 0) then - read(unitn, nml=soilhydrology_inparm, iostat=ierr) - if (ierr /= 0) then - call endrun(msg="ERROR reading "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) - end if - else - call endrun(msg="ERROR could NOT find "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) - end if - call relavu( unitn ) - end if - - call shr_mpi_bcast (baseflow_scalar, mpicom) - - if (masterproc) then - write(iulog,*) ' ' - write(iulog,*) nmlname//' settings:' - write(iulog,nml=soilhydrology_inparm) - write(iulog,*) ' ' - end if - - end subroutine soilhydReadNML - - - !----------------------------------------------------------------------- subroutine SetSoilWaterFractions(bounds, num_hydrologyc, filter_hydrologyc, & soilhydrology_inst, soilstate_inst, waterstatebulk_inst) @@ -757,7 +696,7 @@ subroutine WaterTable(bounds, num_hydrologyc, filter_hydrologyc, & ! use analytical expression for aquifer specific yield rous = watsat(c,nlevsoi) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,nlevsoi))**(-1./bsw(c,nlevsoi))) - rous=max(rous, distparams%aq_sp_yield_min(c)) + rous=max(rous, distparams%aq_sp_yield_min%param_val(col%gridcell(c))) !-- water table is below the soil column -------------------------------------- if(jwt(c) == nlevsoi) then @@ -772,7 +711,7 @@ subroutine WaterTable(bounds, num_hydrologyc, filter_hydrologyc, & ! use analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) - s_y=max(s_y, distparams%aq_sp_yield_min(c)) + s_y=max(s_y, distparams%aq_sp_yield_min%param_val(col%gridcell(c))) qcharge_layer=min(qcharge_tot,(s_y*(zwt(c) - zi(c,j-1))*1.e3)) qcharge_layer=max(qcharge_layer,0._r8) @@ -787,7 +726,7 @@ subroutine WaterTable(bounds, num_hydrologyc, filter_hydrologyc, & ! use analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) - s_y=max(s_y, distparams%aq_sp_yield_min(c)) + s_y=max(s_y, distparams%aq_sp_yield_min%param_val(col%gridcell(c))) qcharge_layer=max(qcharge_tot,-(s_y*(zi(c,j) - zwt(c))*1.e3)) qcharge_layer=min(qcharge_layer,0._r8) @@ -1054,7 +993,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte c = filter_hydrologyc(fc) ! specify maximum drainage rate - q_perch_max = distparams%perched_baseflow_scalar(c) * sin(col%topo_slope(c) * (rpi/180._r8)) + q_perch_max = distparams%perched_baseflow_scalar%param_val(col%gridcell(c)) * sin(col%topo_slope(c) * (rpi/180._r8)) ! if layer containing water table is frozen, compute the following: ! frost table, perched water table, and drainage from perched saturated layer @@ -1087,7 +1026,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte wtsub = 0._r8 q_perch = 0._r8 do k = jwt(c)+1, k_frz - imped=10._r8**(-distparams%e_ice(c)*(0.5_r8*(icefrac(c,k)+icefrac(c,min(nlevsoi, k+1))))) + imped=10._r8**(-distparams%e_ice%param_val(col%gridcell(c))*(0.5_r8*(icefrac(c,k)+icefrac(c,min(nlevsoi, k+1))))) q_perch = q_perch + imped*hksat(c,k)*dzmm(c,k) wtsub = wtsub + dzmm(c,k) end do @@ -1163,7 +1102,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte wtsub = 0._r8 q_perch = 0._r8 do k = k_perch, k_frz - imped=10._r8**(-distparams%e_ice(c)*(0.5_r8*(icefrac(c,k)+icefrac(c,min(nlevsoi, k+1))))) + imped=10._r8**(-distparams%e_ice%param_val(col%gridcell(c))*(0.5_r8*(icefrac(c,k)+icefrac(c,min(nlevsoi, k+1))))) q_perch = q_perch + imped*hksat(c,k)*dzmm(c,k) wtsub = wtsub + dzmm(c,k) end do @@ -1208,11 +1147,11 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte end do ! add ice impedance factor to baseflow if (use_vichydro) then - imped=10._r8**(-distparams%e_ice(c)*min(1.0_r8,ice(c,nlayer)/max_moist(c,nlayer))) + imped=10._r8**(-distparams%e_ice%param_val(col%gridcell(c))*min(1.0_r8,ice(c,nlayer)/max_moist(c,nlayer))) dsmax_tmp(c) = Dsmax(c) * dtime/ secspday !mm/day->mm/dtime rsub_top_max = dsmax_tmp(c) else - imped=10._r8**(-distparams%e_ice(c)*(icefracsum/dzsum)) + imped=10._r8**(-distparams%e_ice%param_val(col%gridcell(c))*(icefracsum/dzsum)) rsub_top_max = 10._r8 * sin((rpi/180.) * col%topo_slope(c)) end if @@ -1237,7 +1176,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte ! use analytical expression for aquifer specific yield rous = watsat(c,nlevsoi) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,nlevsoi))**(-1./bsw(c,nlevsoi))) - rous=max(rous, distparams%aq_sp_yield_min(c)) + rous=max(rous, distparams%aq_sp_yield_min%param_val(col%gridcell(c))) !-- water table is below the soil column -------------------------------------- if(jwt(c) == nlevsoi) then @@ -1274,7 +1213,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte ! use analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) - s_y=max(s_y, distparams%aq_sp_yield_min(c)) + s_y=max(s_y, distparams%aq_sp_yield_min%param_val(col%gridcell(c))) rsub_top_layer=max(rsub_top_tot,-(s_y*(zi(c,j) - zwt(c))*1.e3)) rsub_top_layer=min(rsub_top_layer,0._r8) @@ -1882,7 +1821,7 @@ subroutine PerchedLateralFlow(bounds, num_hydrologyc, & else ! Non-hillslope columns ! specify maximum drainage rate - q_perch_max = distparams%perched_baseflow_scalar(c) & + q_perch_max = distparams%perched_baseflow_scalar%param_val(col%gridcell(c)) & * sin(col%topo_slope(c) * (rpi/180._r8)) wtsub = 0._r8 @@ -1929,7 +1868,7 @@ subroutine PerchedLateralFlow(bounds, num_hydrologyc, & s_y = watsat(c,k) & * ( 1. - (1.+1.e3*zwt_perched(c)/sucsat(c,k))**(-1./bsw(c,k))) - s_y=max(s_y,distparams%aq_sp_yield_min(c)) + s_y=max(s_y,distparams%aq_sp_yield_min%param_val(col%gridcell(c))) if (k==k_perch(c)) then drainage_layer=min(drainage_tot,(s_y*(zi(c,k) - zwt_perched(c))*1.e3)) else @@ -2168,7 +2107,7 @@ subroutine SubsurfaceLateralFlow(bounds, & vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice)) icefrac(c,j) = min(1._r8,vol_ice/watsat(c,j)) - ice_imped(c,j)=10._r8**(-distparams%e_ice(c)*icefrac(c,j)) + ice_imped(c,j)=10._r8**(-distparams%e_ice%param_val(col%gridcell(c))*icefrac(c,j)) end do end do @@ -2210,7 +2149,7 @@ subroutine SubsurfaceLateralFlow(bounds, & dzsum = dzsum + dzmm(c,j) icefracsum = icefracsum + icefrac(c,j) * dzmm(c,j) end do - ice_imped_col(c)=10._r8**(-distparams%e_ice(c)*(icefracsum/dzsum)) + ice_imped_col(c)=10._r8**(-distparams%e_ice%param_val(col%gridcell(c))*(icefracsum/dzsum)) enddo do fc = 1, num_hydrologyc @@ -2358,9 +2297,9 @@ subroutine SubsurfaceLateralFlow(bounds, & ! Non-hillslope columns ! baseflow is power law expression relative to bedrock layer if(zwt(c) <= zi(c,nbedrock(c))) then - qflx_latflow_out(c) = ice_imped_col(c) * distparams%baseflow_scalar(c) & + qflx_latflow_out(c) = ice_imped_col(c) * distparams%baseflow_scalar%param_val(col%gridcell(c)) & * tan(rpi/180._r8*col%topo_slope(c))* & - (zi(c,nbedrock(c)) - zwt(c))**(distparams%n_baseflow(c)) + (zi(c,nbedrock(c)) - zwt(c))**(distparams%n_baseflow%param_val(col%gridcell(c))) endif ! convert flux to volumetric flow qflx_latflow_out_vol(c) = 1.e-3_r8*qflx_latflow_out(c)*(grc%area(g)*1.e6_r8*col%wtgcell(c)) @@ -2423,7 +2362,7 @@ subroutine SubsurfaceLateralFlow(bounds, & ! analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) - s_y=max(s_y,distparams%aq_sp_yield_min(c)) + s_y=max(s_y,distparams%aq_sp_yield_min%param_val(col%gridcell(c))) drainage_layer=min(drainage_tot,(s_y*dz(c,j)*1.e3)) @@ -2450,7 +2389,7 @@ subroutine SubsurfaceLateralFlow(bounds, & ! analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) - s_y=max(s_y,distparams%aq_sp_yield_min(c)) + s_y=max(s_y,distparams%aq_sp_yield_min%param_val(col%gridcell(c))) drainage_layer=max(drainage_tot,-(s_y*(zi(c,j) - zwt(c))*1.e3)) drainage_layer=min(drainage_layer,0._r8) @@ -2789,7 +2728,7 @@ subroutine CalcIrrigWithdrawals(bounds, & ! use analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) - s_y=max(s_y, distparams%aq_sp_yield_min(c)) + s_y=max(s_y, distparams%aq_sp_yield_min%param_val(col%gridcell(c))) if (j==jwt(c)+1) then available_water_layer=max(0._r8,(s_y*(zi(c,j) - zwt(c))*1.e3)) diff --git a/src/biogeophys/SoilStateInitTimeConstMod.F90 b/src/biogeophys/SoilStateInitTimeConstMod.F90 index 96547d2c57..3f5311f8ca 100644 --- a/src/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/biogeophys/SoilStateInitTimeConstMod.F90 @@ -10,7 +10,7 @@ module SoilStateInitTimeConstMod use LandunitType , only : lun use ColumnType , only : col use PatchType , only : patch - use DistParamType , only : distparams + use DistParamType , only : distparams => distributed_parameters use abortUtils , only : endrun use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) ! @@ -548,13 +548,13 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%bd_col(c,lev) = (1._r8 - soilstate_inst%watsat_col(c,lev))*params_inst%pd ! do not allow watsat_sf to push watsat above 0.93 - soilstate_inst%watsat_col(c,lev) = min(distparams%watsat_sf(c) * ( (1._r8 - om_frac) * & + soilstate_inst%watsat_col(c,lev) = min(distparams%watsat_sf%param_val(col%gridcell(c)) * ( (1._r8 - om_frac) * & soilstate_inst%watsat_col(c,lev) + om_watsat*om_frac ), 0.93_r8) tkm = (1._r8-om_frac) * (params_inst%tkd_sand*sand+params_inst%tkd_clay*clay)/ & (sand+clay)+params_inst%tkm_om*om_frac ! W/(m K) - soilstate_inst%bsw_col(c,lev) = distparams%bsw_sf(c) * ( (1._r8-om_frac) * & + soilstate_inst%bsw_col(c,lev) = distparams%bsw_sf%param_val(col%gridcell(c)) * ( (1._r8-om_frac) * & (2.91_r8 + 0.159_r8*clay) + om_frac*om_b ) - soilstate_inst%sucsat_col(c,lev) = distparams%sucsat_sf(c) * ( (1._r8-om_frac) * & + soilstate_inst%sucsat_col(c,lev) = distparams%sucsat_sf%param_val(col%gridcell(c)) * ( (1._r8-om_frac) * & soilstate_inst%sucsat_col(c,lev) + om_sucsat*om_frac ) soilstate_inst%hksat_min_col(c,lev) = xksat @@ -576,7 +576,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) else uncon_hksat = 0._r8 end if - soilstate_inst%hksat_col(c,lev) = distparams%hksat_sf(c) * ( uncon_frac*uncon_hksat + & + soilstate_inst%hksat_col(c,lev) = distparams%hksat_sf%param_val(col%gridcell(c)) * ( uncon_frac*uncon_hksat + & (perc_frac*om_frac)*om_hksat ) soilstate_inst%tkmg_col(c,lev) = tkm ** (1._r8- soilstate_inst%watsat_col(c,lev)) @@ -650,16 +650,16 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) bd = (1._r8-soilstate_inst%watsat_col(c,lev))*params_inst%pd ! do not allow watsat_sf to push watsat above 0.93 - soilstate_inst%watsat_col(c,lev) = min(distparams%watsat_sf(c) * ( (1._r8 - om_frac) * & + soilstate_inst%watsat_col(c,lev) = min(distparams%watsat_sf%param_val(col%gridcell(c)) * ( (1._r8 - om_frac) * & soilstate_inst%watsat_col(c,lev) + om_watsat_lake * om_frac), 0.93_r8) tkm = (1._r8-om_frac)*(params_inst%tkd_sand*sand+params_inst%tkd_clay*clay)/(sand+clay) + & params_inst%tkm_om * om_frac ! W/(m K) - soilstate_inst%bsw_col(c,lev) = distparams%bsw_sf(c) * ( (1._r8-om_frac) * & + soilstate_inst%bsw_col(c,lev) = distparams%bsw_sf%param_val(col%gridcell(c)) * ( (1._r8-om_frac) * & (2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake ) - soilstate_inst%sucsat_col(c,lev) = distparams%sucsat_sf(c) * ( (1._r8-om_frac) * & + soilstate_inst%sucsat_col(c,lev) = distparams%sucsat_sf%param_val(col%gridcell(c)) * ( (1._r8-om_frac) * & soilstate_inst%sucsat_col(c,lev) + om_sucsat_lake * om_frac ) xksat = 0.0070556_r8 *( 10._r8**(-0.884_r8+0.0153_r8*sand) ) ! mm/s @@ -683,7 +683,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) uncon_hksat = 0._r8 end if - soilstate_inst%hksat_col(c,lev) = distparams%hksat_sf(c) * ( uncon_frac*uncon_hksat + & + soilstate_inst%hksat_col(c,lev) = distparams%hksat_sf%param_val(col%gridcell(c)) * ( uncon_frac*uncon_hksat + & (perc_frac*om_frac)*om_hksat_lake ) soilstate_inst%tkmg_col(c,lev) = tkm ** (1._r8- soilstate_inst%watsat_col(c,lev)) soilstate_inst%tksatu_col(c,lev) = soilstate_inst%tkmg_col(c,lev)*0.57_r8**soilstate_inst%watsat_col(c,lev) diff --git a/src/biogeophys/SoilWaterMovementMod.F90 b/src/biogeophys/SoilWaterMovementMod.F90 index 3164de6b15..3c168cd8a2 100644 --- a/src/biogeophys/SoilWaterMovementMod.F90 +++ b/src/biogeophys/SoilWaterMovementMod.F90 @@ -9,7 +9,7 @@ module SoilWaterMovementMod ! created by Jinyun Tang, Mar 12, 2014 use shr_kind_mod , only : r8 => shr_kind_r8 use shr_sys_mod , only : shr_sys_flush - use DistParamType , only : distparams + use DistParamType , only : distparams => distributed_parameters ! implicit none @@ -724,7 +724,7 @@ subroutine soilwater_zengdecker2009(bounds, num_hydrologyc, filter_hydrologyc, & s1 = min(1._r8, s1) s2 = hksat(c,j)*s1**(2._r8*bsw(c,j)+2._r8) - imped(c,j)=10._r8**(-distparams%e_ice(c)*(0.5_r8*(icefrac(c,j)+icefrac(c,min(nlevsoi, j+1))))) + imped(c,j)=10._r8**(-distparams%e_ice%param_val(col%gridcell(c))*(0.5_r8*(icefrac(c,j)+icefrac(c,min(nlevsoi, j+1))))) hk(c,j) = imped(c,j)*s1*s2 dhkdw(c,j) = imped(c,j)*(2._r8*bsw(c,j)+3._r8)*s2* & @@ -2125,6 +2125,8 @@ end subroutine compute_qcharge !----------------------------------------------------------------------- subroutine IceImpedance(c, icefrac, imped) ! + ! USES + use ColumnType , only : col !DESCRIPTION ! compute soil suction potential ! @@ -2143,7 +2145,7 @@ subroutine IceImpedance(c, icefrac, imped) character(len=32) :: subname = 'IceImpedance' ! subroutine name !------------------------------------------------------------------------------ - imped = 10._r8**(-distparams%e_ice(c)*icefrac) + imped = 10._r8**(-distparams%e_ice%param_val(col%gridcell(c))*icefrac) end subroutine IceImpedance diff --git a/src/biogeophys/SurfaceResistanceMod.F90 b/src/biogeophys/SurfaceResistanceMod.F90 index bb54a040b9..64d2bdba79 100644 --- a/src/biogeophys/SurfaceResistanceMod.F90 +++ b/src/biogeophys/SurfaceResistanceMod.F90 @@ -12,7 +12,7 @@ module SurfaceResistanceMod use shr_const_mod , only : SHR_CONST_TKFRZ use clm_varctl , only : iulog use SoilStateType , only : soilstate_type - use DistParamType , only : distparams + use DistParamType , only : distparams => distributed_parameters use WaterStateBulkType, only : waterstatebulk_type use WaterDiagnosticBulkType, only : waterdiagnosticbulk_type use TemperatureType , only : temperature_type @@ -373,9 +373,9 @@ subroutine calc_soil_resistance_sl14(bounds, num_nolakec, filter_nolakec, & ! dsl(c) = dzmm(c,1)*max(0.001_r8,(0.8*eff_porosity(c,1) - vwc_liq)) & ! try arbitrary scaling (not top layer thickness) ! dsl(c) = 15._r8*max(0.001_r8,(0.8*eff_porosity(c,1) - vwc_liq)) & - dsl(c) = distparams%d_max(c) * max(0.001_r8, (distparams%frac_sat_soil_dsl_init(c) * eff_por_top - vwc_liq)) & + dsl(c) = distparams%d_max%param_val(col%gridcell(c)) * max(0.001_r8, (distparams%frac_sat_soil_dsl_init%param_val(col%gridcell(c)) * eff_por_top - vwc_liq)) & ! /max(0.001_r8,(watsat(c,1)- aird)) - / max(0.001_r8, (distparams%frac_sat_soil_dsl_init(c) * watsat(c,1) - aird)) + / max(0.001_r8, (distparams%frac_sat_soil_dsl_init%param_val(col%gridcell(c)) * watsat(c,1) - aird)) dsl(c)=max(dsl(c),0._r8) dsl(c)=min(dsl(c),200._r8) diff --git a/src/biogeophys/SurfaceWaterMod.F90 b/src/biogeophys/SurfaceWaterMod.F90 index 60ce1f67b6..e59471d514 100644 --- a/src/biogeophys/SurfaceWaterMod.F90 +++ b/src/biogeophys/SurfaceWaterMod.F90 @@ -15,7 +15,7 @@ module SurfaceWaterMod use column_varcon , only : icol_roof, icol_road_imperv, icol_sunwall, icol_shadewall, icol_road_perv use decompMod , only : bounds_type, subgrid_level_column use ColumnType , only : col - use DistParamType , only : distparams + use DistParamType , only : distparams => distributed_parameters use NumericsMod , only : truncate_small_values use InfiltrationExcessRunoffMod , only : infiltration_excess_runoff_type use EnergyFluxType , only : energyflux_type @@ -446,10 +446,10 @@ subroutine QflxH2osfcSurf(bounds, num_hydrologyc, filter_hydrologyc, & c = filter_hydrologyc(fc) if (h2osfcflag==1) then - if (frac_h2osfc_nosnow(c) <= distparams%pc(c)) then + if (frac_h2osfc_nosnow(c) <= distparams%pc%param_val(col%gridcell(c))) then frac_infclust=0.0_r8 else - frac_infclust=(frac_h2osfc_nosnow(c)-distparams%pc(c))**distparams%mu(c) + frac_infclust=(frac_h2osfc_nosnow(c)-distparams%pc%param_val(col%gridcell(c)))**distparams%mu%param_val(col%gridcell(c)) endif endif diff --git a/src/biogeophys/WaterDiagnosticBulkType.F90 b/src/biogeophys/WaterDiagnosticBulkType.F90 index bab5ab87fc..b8e9ed4250 100644 --- a/src/biogeophys/WaterDiagnosticBulkType.F90 +++ b/src/biogeophys/WaterDiagnosticBulkType.F90 @@ -21,7 +21,7 @@ module WaterDiagnosticBulkType use clm_varcon , only : spval use LandunitType , only : lun use ColumnType , only : col - use DistParamType , only : distparams + use DistParamType , only : distparams => distributed_parameters use filterColMod , only : filter_col_type, col_filter_from_ltypes use WaterDiagnosticType, only : waterdiagnostic_type use WaterInfoBaseType, only : water_info_base_type @@ -758,7 +758,7 @@ subroutine InitBulkCold(this, bounds, & fmelt = (snowbd/100.)**1. ! 100 is the assumed fresh snow density; 1 is a melting factor that could be ! reconsidered, optimal value of 1.5 in Niu et al., 2007 - this%frac_sno_col(c) = tanh( this%snow_depth_col(c) / (2.5 * distparams%zlnd(c) * fmelt) ) + this%frac_sno_col(c) = tanh( this%snow_depth_col(c) / (2.5 * distparams%zlnd%param_val(col%gridcell(c)) * fmelt) ) endif end if end do diff --git a/src/cpl/share_esmf/DistParamStreamMod.F90 b/src/cpl/share_esmf/DistParamStreamMod.F90 index 83f7197c57..ea57cb956f 100644 --- a/src/cpl/share_esmf/DistParamStreamMod.F90 +++ b/src/cpl/share_esmf/DistParamStreamMod.F90 @@ -238,8 +238,10 @@ subroutine Init(this, bounds, NLFilename, distparams) call AssignDistributedParameter(distparams%frac_sat_soil_dsl_init,bounds,dataptr1d,stream_varnames(n)) call AssignDistributedParameter(distparams%zlnd,bounds,dataptr1d,stream_varnames(n)) call AssignDistributedParameter(distparams%snw_rds_min,bounds,dataptr1d,stream_varnames(n)) - call AssignDistributedParameter(distparams%precip_repartition_nonglc_all_rain_t,bounds,dataptr1d,stream_varnames(n)) - call AssignDistributedParameter(distparams%precip_repartition_nonglc_all_snow_t,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%precip_repartition_nonglc_all_rain_t_celsius,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%precip_repartition_nonglc_all_snow_t_celsius,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%precip_repartition_glc_all_rain_t_celsius,bounds,dataptr1d,stream_varnames(n)) + call AssignDistributedParameter(distparams%precip_repartition_glc_all_snow_t_celsius,bounds,dataptr1d,stream_varnames(n)) end do deallocate(stream_varnames) diff --git a/src/main/DistParamType.F90 b/src/main/DistParamType.F90 index 8a9cec4a56..e17a89cd04 100644 --- a/src/main/DistParamType.F90 +++ b/src/main/DistParamType.F90 @@ -10,7 +10,7 @@ module DistParamType use shr_sys_mod , only : shr_sys_abort use clm_varcon , only : spval, ispval, grlnd use clm_varctl , only : iulog - use clm_varctl , only : paramfile, distributed_paramfile + use clm_varctl , only : paramfile use ColumnType , only : col use spmdMod , only : masterproc, mpicom use fileutils , only : getfil, opnfil, getavu, relavu @@ -19,6 +19,7 @@ module DistParamType use decompMod , only : bounds_type use ncdio_pio , only : ncd_pio_closefile, ncd_pio_openfile use ncdio_pio , only : file_desc_t, ncd_inqdid, ncd_inqdlen + ! ! !PUBLIC TYPES: implicit none @@ -26,145 +27,122 @@ module DistParamType private character(len=*), parameter, private :: sourcefile = __FILE__ ! - type, public :: distparam_type - ! CanopyHydrologyMod - real(r8), pointer :: liq_canopy_storage_scalar (:) ! Canopy-storage-of-liquid-water parameter (kg/m2) - real(r8), pointer :: snow_canopy_storage_scalar (:) ! Canopy-storage-of-snow parameter (kg/m2) - real(r8), pointer :: snowcan_unload_temp_fact (:) ! Temperature canopy snow unload scaling (C2 in Eq. 14, Roesch et al. 2001) (K*s) - real(r8), pointer :: snowcan_unload_wind_fact (:) ! Wind canopy snow unload scaling (modifies 1.56e5, where 1.56e5 is C3 in Eq. 15, Roesch et al. 2001) (-) - real(r8), pointer :: interception_fraction (:) ! Fraction of intercepted precipitation (-) - real(r8), pointer :: maximum_leaf_wetted_fraction (:) ! Maximum fraction of leaf that may be wet (-) + + type, public :: distparam_class + logical :: is_distributed = .false. ! is parameter spatially distributed? + character(len=256) :: name ! name on parameter files + real(r8), allocatable :: val(:) ! parameter value array (either distributed or length 1) + contains + procedure, public :: param_val ! return parameter value + end type distparam_class + + type, public :: distributed_parameter_type + ! all parameters that could *potentially* be spatially distributed ! SoilHydrologyMod - real(r8), pointer :: aq_sp_yield_min (:) ! Minimum aquifer specific yield (unitless) - real(r8), pointer :: n_baseflow (:) ! Drainage power law exponent (unitless) - real(r8), pointer :: perched_baseflow_scalar (:) ! Scalar multiplier for perched base flow rate (kg/m2/s) - real(r8), pointer :: e_ice (:) ! Soil ice impedance factor (unitless) - real(r8), pointer :: baseflow_scalar (:) ! Scalar multiplier for base flow rate () + class(distparam_class), pointer :: aq_sp_yield_min => NULL() ! Minimum aquifer specific yield (unitless) + class(distparam_class), pointer :: n_baseflow => NULL() ! Drainage power law exponent (unitless) + class(distparam_class), pointer :: perched_baseflow_scalar => NULL() ! Scalar multiplier for perched base flow rate (kg/m2/s) + class(distparam_class), pointer :: e_ice => NULL() ! Soil ice impedance factor (unitless) + class(distparam_class), pointer :: baseflow_scalar => NULL() ! Scalar multiplier for base flow rate () ! SaturatedExcessRunoff - real(r8), pointer :: fff (:) ! Decay factor for fractional saturated area (1/m) + class(distparam_class), pointer :: fff => NULL() ! Decay factor for fractional saturated area (1/m) ! initVerticalMod - real(r8), pointer :: slopebeta (:) ! exponent for microtopography pdf sigma (unitless) - real(r8), pointer :: slopemax (:) ! max topographic slope for microtopography pdf sigma (unitless) + class(distparam_class), pointer :: slopebeta => NULL() ! exponent for microtopography pdf sigma (unitless) + class(distparam_class), pointer :: slopemax => NULL() ! max topographic slope for microtopography pdf sigma (unitless) ! SnowCoverFractionSwensonLawrence2012Mod - real(r8), pointer :: n_melt_coef (:) ! SCA shape parameter - real(r8), pointer :: accum_factor (:) ! Accumulation constant for fractional snow covered area (unitless) + class(distparam_class), pointer :: n_melt_coef => NULL() ! SCA shape parameter + class(distparam_class), pointer :: accum_factor => NULL() ! Accumulation constant for fractional snow covered area (unitless) ! SnowHydrologyMod - real(r8), pointer :: upplim_destruct_metamorph (:) ! Upper limit on destructive metamorphism compaction (kg/m3) + class(distparam_class), pointer :: upplim_destruct_metamorph => NULL() ! Upper limit on destructive metamorphism compaction (kg/m3) ! SoilStateInitTimeConstMod - real(r8), pointer :: bsw_sf (:) ! Scale factor for bsw (unitless) - real(r8), pointer :: hksat_sf (:) ! Scale factor for hksat (unitless) - real(r8), pointer :: sucsat_sf (:) ! Scale factor for sucsat (unitless) - real(r8), pointer :: watsat_sf (:) ! Scale factor for watsat (unitless) + class(distparam_class), pointer :: bsw_sf => NULL() ! Scale factor for bsw (unitless) + class(distparam_class), pointer :: hksat_sf => NULL() ! Scale factor for hksat (unitless) + class(distparam_class), pointer :: sucsat_sf => NULL() ! Scale factor for sucsat (unitless) + class(distparam_class), pointer :: watsat_sf => NULL() ! Scale factor for watsat (unitless) ! SurfaceResistanceMod - real(r8), pointer :: d_max (:) ! Dry surface layer parameter (mm) - real(r8), pointer :: frac_sat_soil_dsl_init (:) ! Fraction of saturated soil for moisture value at which DSL initiates (unitless) + class(distparam_class), pointer :: d_max => NULL() ! Dry surface layer parameter (mm) + class(distparam_class), pointer :: frac_sat_soil_dsl_init => NULL() ! Fraction of saturated soil for moisture value at which DSL initiates (unitless) ! SurfaceWaterMod - real(r8), pointer :: pc (:) ! Threshold probability for surface water (unitless) - real(r8), pointer :: mu (:) ! Connectivity exponent for surface water (unitless) + class(distparam_class), pointer :: pc => NULL() ! Threshold probability for surface water (unitless) + class(distparam_class), pointer :: mu => NULL() ! Connectivity exponent for surface water (unitless) ! WaterDiagnosticBulkType - real(r8), pointer :: zlnd (:) ! Momentum roughness length for soil, glacier, wetland (m) - real(r8), pointer :: snw_rds_min (:) ! minimum allowed snow effective radius (also cold "fresh snow" value) [microns] + class(distparam_class), pointer :: zlnd => NULL() ! Momentum roughness length for soil, glacier, wetland (m) + class(distparam_class), pointer :: snw_rds_min => NULL() ! minimum allowed snow effective radius (also cold "fresh snow" value) [microns] ! atm2lndType - real(r8), pointer :: precip_repartition_nonglc_all_rain_t (:) ! Rain temperature threshold for non-glacier landunits (C) - real(r8), pointer :: precip_repartition_nonglc_all_snow_t (:) ! Snow temperature threshold for non-glacier landunits (C) - + class(distparam_class), pointer :: precip_repartition_nonglc_all_rain_t_celsius => NULL() ! Rain temperature threshold for non-glacier landunits (C) + class(distparam_class), pointer :: precip_repartition_nonglc_all_snow_t_celsius => NULL() ! Snow temperature threshold for non-glacier landunits (C) + class(distparam_class), pointer :: precip_repartition_glc_all_rain_t_celsius => NULL() ! Rain temperature threshold for glacier landunits (C) + class(distparam_class), pointer :: precip_repartition_glc_all_snow_t_celsius => NULL() ! Snow temperature threshold for glacier landunits (C) + contains procedure, public :: Init - procedure, public :: readDistributedParams procedure, public :: Clean - end type distparam_type + end type distributed_parameter_type + + type(distributed_parameter_type), public, target :: distributed_parameters + + public :: InitGlobalParameters - type(distparam_type), public, target :: distparams ! - !------------------------------------------------------------------------ contains !------------------------------------------------------------------------ - subroutine Init(this, bounds) + function param_val(this,g) + ! + ! !DESCRIPTION: + ! Return distributed parameter value if allocated, otherwise return scalar value ! ! !ARGUMENTS: - class(distparam_type) :: this - type(bounds_type), intent(in) :: bounds + implicit none + class(distparam_class) :: this + real(r8) :: param_val + integer :: g ! - ! !LOCAL VARIABLES: - integer :: begc, endc - integer :: begg, endg - !------------------------------------------------------------------------ - - begc = bounds%begc; endc = bounds%endc - begg = bounds%begg; endg = bounds%endg - - ! these should be patch; todo - ! CanopyHydrology - allocate(this%liq_canopy_storage_scalar (begc:endc)) ; this%liq_canopy_storage_scalar (:) = nan - allocate(this%snow_canopy_storage_scalar (begc:endc)) ; this%snow_canopy_storage_scalar (:) = nan - allocate(this%snowcan_unload_temp_fact (begc:endc)) ; this%snowcan_unload_temp_fact (:) = nan - allocate(this%snowcan_unload_wind_fact (begc:endc)) ; this%snowcan_unload_wind_fact (:) = nan - allocate(this%interception_fraction (begc:endc)) ; this%interception_fraction (:) = nan - allocate(this%maximum_leaf_wetted_fraction (begc:endc)) ; this%maximum_leaf_wetted_fraction (:) = nan - - ! SoilHydrology - allocate(this%aq_sp_yield_min (begc:endc)) ; this%aq_sp_yield_min (:) = nan - allocate(this%n_baseflow (begc:endc)) ; this%n_baseflow (:) = nan - allocate(this%perched_baseflow_scalar (begc:endc)) ; this%perched_baseflow_scalar (:) = nan - allocate(this%e_ice (begc:endc)) ; this%e_ice (:) = nan - allocate(this%baseflow_scalar (begc:endc)) ; this%baseflow_scalar (:) = nan - - ! SaturatedExcessRunoff - allocate(this%fff (begc:endc)) ; this%fff (:) = nan - - ! initVertical - allocate(this%slopebeta (begg:endc)) ; this%slopebeta (:) = nan - allocate(this%slopemax (begg:endc)) ; this%slopemax (:) = nan - - ! SnowCoverFractionSwensonLawrence2012Mod - allocate(this%n_melt_coef (begg:endc)) ; this%n_melt_coef (:) = nan - allocate(this%accum_factor (begg:endc)) ; this%accum_factor (:) = nan - - ! SnowHydrologyMod - allocate(this%upplim_destruct_metamorph (begg:endc)) ; this%upplim_destruct_metamorph (:) = nan - - ! SoilStateInitTimeConstMod - allocate(this%bsw_sf (begg:endc)) ; this%bsw_sf (:) = nan - allocate(this%hksat_sf (begg:endc)) ; this%hksat_sf (:) = nan - allocate(this%sucsat_sf (begg:endc)) ; this%sucsat_sf (:) = nan - allocate(this%watsat_sf (begg:endc)) ; this%watsat_sf (:) = nan - - ! SurfaceResistanceMod - allocate(this%d_max (begg:endc)) ; this%d_max (:) = nan - allocate(this%frac_sat_soil_dsl_init (begg:endc)) ; this%frac_sat_soil_dsl_init (:) = nan - - ! SurfaceWaterMod - allocate(this%pc (begg:endc)) ; this%pc (:) = nan - allocate(this%mu (begg:endc)) ; this%mu (:) = nan - - ! WaterDiagnosticBulkType - allocate(this%zlnd (begg:endc)) ; this%zlnd (:) = nan - allocate(this%snw_rds_min (begg:endc)) ; this%snw_rds_min (:) = nan - - ! atm2lndType - allocate(this%precip_repartition_nonglc_all_rain_t (begg:endc)) ; this%precip_repartition_nonglc_all_rain_t (:) = nan - allocate(this%precip_repartition_nonglc_all_snow_t (begg:endc)) ; this%precip_repartition_nonglc_all_snow_t (:) = nan + if ( this%is_distributed )then + param_val = this%val(g) + else + param_val = this%val(1) + end if + end function param_val - ! allocate(this% (begg:endc)) ; this% (:) = nan + !------------------------------------------------------------------------ + subroutine ReadScalarParameter(this,ncid) + ! + ! !USES: + use paramUtilMod , only : readNcdioScalar + ! !ARGUMENTS: + implicit none + class(distparam_class), intent(inout) :: this + type(file_desc_t) , intent(inout) :: ncid ! pio netCDF file id + ! + ! !LOCAL VARIABLES: + real(r8) :: fscalar_in ! read in - scalar - float + character(len=*), parameter :: subname = 'ReadScalarParameter' + !-------------------------------------------------------------------- - end subroutine Init + if ( .not. this%is_distributed ) then + allocate(this%val(1)) + call readNcdioScalar(ncid, this%name, subname, fscalar_in) + this%val(1) = fscalar_in + endif - !----------------------------------------------------------------------- - subroutine readDistributedParams(this,bounds) + end subroutine ReadScalarParameter + + !------------------------------------------------------------------------ + subroutine InitGlobalParameters(bounds) ! ! !USES: use ncdio_pio , only : check_var, ncd_io @@ -177,721 +155,288 @@ subroutine readDistributedParams(this,bounds) ! ! !ARGUMENTS: implicit none - class(distparam_type) :: this type(bounds_type), intent(in) :: bounds ! ! !LOCAL VARIABLES: logical :: readvar ! whether the variable was found integer :: c, g - integer :: dimid ! netCDF dimension id - integer :: dim1, dim2 ! dimension to compare - ! namelist variable names must match file - real(r8) :: baseflow_scalar ! read in namelist - real(r8) :: precip_repartition_nonglc_all_rain_t ! read in namelist - real(r8) :: precip_repartition_nonglc_all_snow_t ! read in namelist real(r8) :: fscalar_in ! read in - scalar - float real(r8), pointer :: fparam_in(:) ! read in - 1D - float - type(file_desc_t) :: ncids, ncidd ! pio netCDF file ids - character(len=256) :: locfns, locfnd ! local filenames - integer :: ierr ! error code - integer :: unitn ! unit for namelist file - character(len=256) :: nmlname ! namelist - character(len=*), parameter :: subname = 'readDistributedParams' + type(file_desc_t) :: ncid ! pio netCDF file id + character(len=256) :: locfn ! local filename + character(len=*), parameter :: subname = 'InitGlobalParameters' !-------------------------------------------------------------------- if (masterproc) then - write(iulog,*) trim(subname)//' :: reading CLM spatially distributed parameters' + write(iulog,*) trim(subname)//' :: reading CLM parameters' end if - ! global (scalar) parameters - call getfil (paramfile, locfns, 0) - call ncd_pio_openfile (ncids, trim(locfns), 0) - ! spatially distributed parameters - call getfil (distributed_paramfile, locfnd, 0) - call ncd_pio_openfile (ncidd, trim(locfnd), 0) - - ! should a check be added to make sure dimensions match surface data / mesh files? - call ncd_inqdid(ncidd,'lon',dimid) - call ncd_inqdlen(ncidd,dimid,dim1) - call ncd_inqdid(ncidd,'lat',dimid) - call ncd_inqdlen(ncidd,dimid,dim2) - if (dim1 /= ldomain%ni .or. dim2 /= ldomain%nj) then - call endrun(msg="ERROR dimensions of distributed parameter file incompatible with domain dimensions"//errmsg(sourcefile, __LINE__)) - endif - - !----------------------------------------------------------- - ! CanopyHydrology ! - !----------------------------------------------------------- - - ! these should be patch/pft not column... leaving in for now. - - ! Canopy-storage-of-liquid-water parameter (kg/m2) - call check_var(ncidd, 'liq_canopy_storage_scalar', readvar) - if (.false.) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%liq_canopy_storage_scalar(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'liq_canopy_storage_scalar', subname, fscalar_in) - this%liq_canopy_storage_scalar(:) = fscalar_in - endif - - ! Canopy-storage-of-snow parameter (kg/m2) - call check_var(ncidd, 'snow_canopy_storage_scalar', readvar) - if (.false.) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%snow_canopy_storage_scalar(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'snow_canopy_storage_scalar', subname, fscalar_in) - this%snow_canopy_storage_scalar(:) = fscalar_in - endif - - ! Temperature canopy snow unload scaling (C2 in Eq. 14, Roesch et al. 2001) (K*s) - call check_var(ncidd, 'snowcan_unload_temp_fact', readvar) - if (.false.) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='snowcan_unload_temp_fact', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%snowcan_unload_temp_fact(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'snowcan_unload_temp_fact', subname, fscalar_in) - this%snowcan_unload_temp_fact(:) = fscalar_in - endif + !distributed_parameter_stream%Init must be called before this routine is called. - ! Wind canopy snow unload scaling (modifies 1.56e5, where 1.56e5 is C3 in Eq. 15, Roesch et al. 2001) (-) - call check_var(ncidd, 'snowcan_unload_wind_fact', readvar) - if (.false.) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='snowcan_unload_wind_fact', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%snowcan_unload_wind_fact(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'snowcan_unload_wind_fact', subname, fscalar_in) - this%snowcan_unload_wind_fact(:) = fscalar_in - endif - - ! Fraction of intercepted precipitation (-) - call check_var(ncidd, 'interception_fraction', readvar) - if (.false.) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='interception_fraction', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%interception_fraction(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'interception_fraction', subname, fscalar_in) - this%interception_fraction(:) = fscalar_in - endif - - ! Maximum fraction of leaf that may be wet (-) - call check_var(ncidd, 'maximum_leaf_wetted_fraction', readvar) - if (.false.) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='maximum_leaf_wetted_fraction', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%maximum_leaf_wetted_fraction(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'maximum_leaf_wetted_fraction', subname, fscalar_in) - this%maximum_leaf_wetted_fraction(:) = fscalar_in - endif + ! global (scalar) parameters + call getfil (paramfile, locfn, 0) + call ncd_pio_openfile (ncid, trim(locfn), 0) + ! Any parameters not read from stream will be read here from parameter file or namelist !----------------------------------------------------------- ! SoilHydrology ! !----------------------------------------------------------- ! Minimum aquifer specific yield (unitless) - call check_var(ncidd, 'aq_sp_yield_min', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='aq_sp_yield_min', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%aq_sp_yield_min(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'aq_sp_yield_min', subname, fscalar_in) - this%aq_sp_yield_min(:) = fscalar_in - endif - + call ReadScalarParameter(distributed_parameters%aq_sp_yield_min,ncid) + ! Drainage power law exponent (unitless) - call check_var(ncidd, 'n_baseflow', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='n_baseflow', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%n_baseflow(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'n_baseflow', subname, fscalar_in) - this%n_baseflow(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%n_baseflow,ncid) + + ! Scalar multiplier for perched base flow rate () + call ReadScalarParameter(distributed_parameters%baseflow_scalar,ncid) ! Scalar multiplier for perched base flow rate (kg/m2/s) - call check_var(ncidd, 'perched_baseflow_scalar', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='perched_baseflow_scalar', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%perched_baseflow_scalar(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'perched_baseflow_scalar', subname, fscalar_in) - this%perched_baseflow_scalar(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%perched_baseflow_scalar,ncid) ! Soil ice impedance factor (unitless) - call check_var(ncidd, 'e_ice', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='e_ice', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%e_ice(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'e_ice', subname, fscalar_in) - this%e_ice(:) = fscalar_in - endif - - ! Scalar multiplier for base flow rate () - call check_var(ncidd, 'baseflow_scalar', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='baseflow_scalar', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%baseflow_scalar(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! global parameter value from namelist - nmlname = 'soilhydrology_inparm' - namelist /soilhydrology_inparm/ baseflow_scalar - - if (masterproc) then - unitn = getavu() - write(iulog,*) 'Read in '//nmlname//' namelist' - call opnfil (NLFilename, unitn, 'F') - call shr_nl_find_group_name(unitn, nmlname, status=ierr) - if (ierr == 0) then - read(unitn, nml=soilhydrology_inparm, iostat=ierr) - if (ierr /= 0) then - call endrun(msg="ERROR reading "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) - end if - else - call endrun(msg="ERROR could NOT find "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) - end if - call relavu( unitn ) - end if - - call shr_mpi_bcast (baseflow_scalar, mpicom) - this%baseflow_scalar(:) = baseflow_scalar - endif + call ReadScalarParameter(distributed_parameters%e_ice,ncid) !----------------------------------------------------------- ! SaturatedExcessRunoff ! !----------------------------------------------------------- + ! Decay factor for fractional saturated area (1/m) - call check_var(ncidd, 'fff', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='fff', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%fff(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'fff', subname, fscalar_in) - this%fff(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%fff,ncid) - !----------------------------------------------------------- ! initVertical ! !----------------------------------------------------------- ! exponent for microtopography pdf sigma (unitless) - call check_var(ncidd, 'slopebeta', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='slopebeta', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%slopebeta(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'slopebeta', subname, fscalar_in) - this%slopebeta(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%slopebeta,ncid) ! max topographic slope for microtopography pdf sigma (unitless) - call check_var(ncidd, 'slopemax', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='slopemax', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%slopemax(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'slopemax', subname, fscalar_in) - this%slopemax(:) = fscalar_in - endif - + call ReadScalarParameter(distributed_parameters%slopemax,ncid) + !----------------------------------------------------------- ! SnowCoverFractionSwensonLawrence2012 ! !----------------------------------------------------------- ! SCA shape parameter - call check_var(ncidd, 'n_melt_coef', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='n_melt_coef', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%n_melt_coef(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'n_melt_coef', subname, fscalar_in) - this%n_melt_coef(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%n_melt_coef,ncid) ! Accumulation constant for fractional snow covered area (unitless) - call check_var(ncidd, 'accum_factor', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='accum_factor', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%accum_factor(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'accum_factor', subname, fscalar_in) - this%accum_factor(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%accum_factor,ncid) !----------------------------------------------------------- ! SnowHydrologyMod ! !----------------------------------------------------------- ! Upper limit on destructive metamorphism compaction (kg/m3) - call check_var(ncidd, 'upplim_destruct_metamorph', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='upplim_destruct_metamorph', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%upplim_destruct_metamorph(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'upplim_destruct_metamorph', subname, fscalar_in) - this%upplim_destruct_metamorph(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%upplim_destruct_metamorph,ncid) !----------------------------------------------------------- ! SoilStateInitTimeConstMod ! !----------------------------------------------------------- - call check_var(ncidd, 'bsw_sf', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='bsw_sf', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%bsw_sf(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'bsw_sf', subname, fscalar_in) - this%bsw_sf(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%bsw_sf,ncid) - call check_var(ncidd, 'hksat_sf', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='hksat_sf', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%hksat_sf(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'hksat_sf', subname, fscalar_in) - this%hksat_sf(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%hksat_sf,ncid) - call check_var(ncidd, 'sucsat_sf', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='sucsat_sf', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%sucsat_sf(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'sucsat_sf', subname, fscalar_in) - this%sucsat_sf(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%sucsat_sf,ncid) - call check_var(ncidd, 'watsat_sf', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='watsat_sf', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%watsat_sf(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'watsat_sf', subname, fscalar_in) - this%watsat_sf(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%watsat_sf,ncid) !----------------------------------------------------------- ! SurfaceResistanceMod ! !----------------------------------------------------------- ! Dry surface layer parameter (mm) - call check_var(ncidd, 'd_max', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='d_max', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%d_max(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'd_max', subname, fscalar_in) - this%d_max(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%d_max,ncid) ! Fraction of saturated soil for moisture value at which DSL initiates (unitless) - call check_var(ncidd, 'frac_sat_soil_dsl_init', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='frac_sat_soil_dsl_init', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%frac_sat_soil_dsl_init(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'frac_sat_soil_dsl_init', subname, fscalar_in) - this%frac_sat_soil_dsl_init(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%frac_sat_soil_dsl_init,ncid) !----------------------------------------------------------- ! SurfaceWaterMod ! !----------------------------------------------------------- - ! Threshold probability for surface water (unitless) - call check_var(ncidd, 'pc', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='pc', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%pc(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'pc', subname, fscalar_in) - this%pc(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%pc,ncid) ! Connectivity exponent for surface water (unitless) - call check_var(ncidd, 'mu', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='mu', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%mu(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'mu', subname, fscalar_in) - this%mu(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%mu,ncid) !----------------------------------------------------------- ! WaterDiagnosticBulkType ! !----------------------------------------------------------- ! Momentum roughness length for soil, glacier, wetland (m) - call check_var(ncidd, 'zlnd', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='zlnd', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%zlnd(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'zlnd', subname, fscalar_in) - this%zlnd(:) = fscalar_in - endif - - ! minimum allowed snow effective radius (also cold "fresh snow" value) [microns] - call check_var(ncidd, 'snw_rds_min', readvar) - if (.false.) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='snw_rds_min', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%snw_rds_min(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! use global parameter value - call readNcdioScalar(ncids, 'snw_rds_min', subname, fscalar_in) - this%snw_rds_min(:) = fscalar_in - endif + call ReadScalarParameter(distributed_parameters%zlnd,ncid) + ! Minimum allowed snow effective radius (also cold "fresh snow" value) [microns] + call ReadScalarParameter(distributed_parameters%snw_rds_min,ncid) + !----------------------------------------------------------- ! atm2lndType ! !----------------------------------------------------------- - call check_var(ncidd, 'precip_repartition_nonglc_all_snow_t', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='precip_repartition_nonglc_all_snow_t', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%precip_repartition_nonglc_all_snow_t(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! global parameter value from namelist - nmlname = 'atm2lnd_inparm' - namelist /atm2lnd_inparm/ precip_repartition_nonglc_all_snow_t - - if (masterproc) then - unitn = getavu() - write(iulog,*) 'Read in '//nmlname//' namelist' - call opnfil (NLFilename, unitn, 'F') - call shr_nl_find_group_name(unitn, nmlname, status=ierr) - if (ierr == 0) then - read(unitn, nml=atm2lnd_inparm, iostat=ierr) - if (ierr /= 0) then - call endrun(msg="ERROR reading "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) - end if - else - call endrun(msg="ERROR could NOT find "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) - end if - call relavu( unitn ) - end if - - call shr_mpi_bcast (precip_repartition_nonglc_all_snow_t, mpicom) - this%precip_repartition_nonglc_all_snow_t(:) = precip_repartition_nonglc_all_snow_t - endif - - call check_var(ncidd, 'precip_repartition_nonglc_all_rain_t', readvar) - if (readvar) then - ! if distributed parameter found, overwrite global value - allocate(fparam_in(bounds%begg:bounds%endg)) - call ncd_io(ncid=ncidd, varname='precip_repartition_nonglc_all_rain_t', flag='read', data=fparam_in, dim1name=grlnd, readvar=readvar) - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - this%precip_repartition_nonglc_all_rain_t(c) = fparam_in(g) - enddo - - deallocate(fparam_in) - else ! global parameter value from namelist - nmlname = 'atm2lnd_inparm' - namelist /atm2lnd_inparm/ precip_repartition_nonglc_all_rain_t - - if (masterproc) then - unitn = getavu() - write(iulog,*) 'Read in '//nmlname//' namelist' - call opnfil (NLFilename, unitn, 'F') - call shr_nl_find_group_name(unitn, nmlname, status=ierr) - if (ierr == 0) then - read(unitn, nml=atm2lnd_inparm, iostat=ierr) - if (ierr /= 0) then - call endrun(msg="ERROR reading "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) - end if - else - call endrun(msg="ERROR could NOT find "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) - end if - call relavu( unitn ) - end if - - call shr_mpi_bcast (precip_repartition_nonglc_all_rain_t, mpicom) - this%precip_repartition_nonglc_all_rain_t(:) = precip_repartition_nonglc_all_rain_t - endif - ! close parameter files - call ncd_pio_closefile(ncids) - call ncd_pio_closefile(ncidd) - - end subroutine readDistributedParams + ! non-glacier all rain temperature (degrees C) + call ReadScalarParameter(distributed_parameters%precip_repartition_nonglc_all_rain_t_celsius,ncid) + ! non-glacier all snow temperature (degrees C) + call ReadScalarParameter(distributed_parameters%precip_repartition_nonglc_all_snow_t_celsius,ncid) + ! glacier all rain temperature (degrees C) + call ReadScalarParameter(distributed_parameters%precip_repartition_glc_all_rain_t_celsius,ncid) + ! glacier all snow temperature (degrees C) + call ReadScalarParameter(distributed_parameters%precip_repartition_glc_all_snow_t_celsius,ncid) + + ! close parameter file + call ncd_pio_closefile(ncid) + + end subroutine InitGlobalParameters !------------------------------------------------------------------------ - subroutine Clean(this) + subroutine Init(this) ! ! !ARGUMENTS: - class(distparam_type) :: this + class(distributed_parameter_type) :: this !------------------------------------------------------------------------ - deallocate(this%liq_canopy_storage_scalar) - deallocate(this%snow_canopy_storage_scalar) - deallocate(this%snowcan_unload_temp_fact) - deallocate(this%snowcan_unload_wind_fact) - deallocate(this%interception_fraction) - deallocate(this%maximum_leaf_wetted_fraction) + allocate(this%aq_sp_yield_min) + allocate(this%n_baseflow) + allocate(this%perched_baseflow_scalar) + allocate(this%e_ice) + allocate(this%baseflow_scalar) + allocate(this%fff) + allocate(this%slopebeta) + allocate(this%slopemax) + allocate(this%n_melt_coef) + allocate(this%accum_factor) + allocate(this%upplim_destruct_metamorph) + allocate(this%pc) + allocate(this%mu) + allocate(this%bsw_sf) + allocate(this%hksat_sf) + allocate(this%sucsat_sf) + allocate(this%watsat_sf) + allocate(this%d_max) + allocate(this%frac_sat_soil_dsl_init) + allocate(this%zlnd) + allocate(this%snw_rds_min) + allocate(this%precip_repartition_nonglc_all_rain_t_celsius) + allocate(this%precip_repartition_nonglc_all_snow_t_celsius) + allocate(this%precip_repartition_glc_all_rain_t_celsius) + allocate(this%precip_repartition_glc_all_snow_t_celsius) - deallocate(this%aq_sp_yield_min) - deallocate(this%n_baseflow) - deallocate(this%perched_baseflow_scalar) - deallocate(this%e_ice) + this%aq_sp_yield_min%name = 'aq_sp_yield_min' + this%aq_sp_yield_min%is_distributed = .false. - deallocate(this%fff) + this%n_baseflow%name = 'n_baseflow' + this%n_baseflow%is_distributed = .false. - deallocate(this%slopebeta) - deallocate(this%slopemax) + this%perched_baseflow_scalar%name = 'perched_baseflow_scalar' + this%perched_baseflow_scalar%is_distributed = .false. - deallocate(this%n_melt_coef) - deallocate(this%accum_factor) + this%e_ice%name = 'e_ice' + this%e_ice%is_distributed = .false. - deallocate(this%upplim_destruct_metamorph) + this%baseflow_scalar%name = 'baseflow_scalar' + this%baseflow_scalar%is_distributed = .false. - deallocate(this%pc) - deallocate(this%mu) + this%fff%name = 'fff' + this%fff%is_distributed = .false. - deallocate(this%bsw_sf) - deallocate(this%hksat_sf) - deallocate(this%sucsat_sf) - deallocate(this%watsat_sf) + this%slopebeta%name = 'slopebeta' + this%slopebeta%is_distributed = .false. - deallocate(this%d_max) - deallocate(this%frac_sat_soil_dsl_init) - - deallocate(this%zlnd) - deallocate(this%snw_rds_min) + this%slopemax%name = 'slopemax' + this%slopemax%is_distributed = .false. + + this%n_melt_coef%name = 'n_melt_coef' + this%n_melt_coef%is_distributed = .false. + + this%accum_factor%name = 'accum_factor' + this%accum_factor%is_distributed = .false. + + this%upplim_destruct_metamorph%name = 'upplim_destruct_metamorph' + this%upplim_destruct_metamorph%is_distributed = .false. + + this%pc%name = 'pc' + this%pc%is_distributed = .false. + + this%mu%name = 'mu' + this%mu%is_distributed = .false. + + this%bsw_sf%name = 'bsw_sf' + this%bsw_sf%is_distributed = .false. + + this%hksat_sf%name = 'hksat_sf' + this%hksat_sf%is_distributed = .false. + + this%sucsat_sf%name = 'sucsat_sf' + this%sucsat_sf%is_distributed = .false. - deallocate(this%precip_repartition_nonglc_all_rain_t) - deallocate(this%precip_repartition_nonglc_all_snow_t) + this%watsat_sf%name = 'watsat_sf' + this%watsat_sf%is_distributed = .false. - ! deallocate(this%) + this%d_max%name = 'd_max' + this%d_max%is_distributed = .false. + + this%frac_sat_soil_dsl_init%name = 'frac_sat_soil_dsl_init' + this%frac_sat_soil_dsl_init%is_distributed = .false. + + this%zlnd%name = 'zlnd' + this%zlnd%is_distributed = .false. + + this%snw_rds_min%name = 'snw_rds_min' + this%snw_rds_min%is_distributed = .false. + + this%precip_repartition_nonglc_all_rain_t_celsius%name = & + 'precip_repartition_nonglc_all_rain_t' + this%precip_repartition_nonglc_all_rain_t_celsius%is_distributed = .false. + + this%precip_repartition_nonglc_all_snow_t_celsius%name = & + 'precip_repartition_nonglc_all_snow_t' + this%precip_repartition_nonglc_all_snow_t_celsius%is_distributed = .false. + + this%precip_repartition_glc_all_rain_t_celsius%name = & + 'precip_repartition_glc_all_rain_t' + this%precip_repartition_glc_all_rain_t_celsius%is_distributed = .false. + + this%precip_repartition_glc_all_snow_t_celsius%name = & + 'precip_repartition_glc_all_snow_t' + this%precip_repartition_glc_all_snow_t_celsius%is_distributed = .false. + + end subroutine Init + + !------------------------------------------------------------------------ + subroutine Clean(this) + ! + ! !ARGUMENTS: + class(distributed_parameter_type) :: this + !------------------------------------------------------------------------ + deallocate(this%aq_sp_yield_min%val) + deallocate(this%n_baseflow%val) + deallocate(this%perched_baseflow_scalar%val) + deallocate(this%e_ice%val) + deallocate(this%baseflow_scalar%val) + deallocate(this%fff%val) + deallocate(this%slopebeta%val) + deallocate(this%slopemax%val) + deallocate(this%n_melt_coef%val) + deallocate(this%accum_factor%val) + deallocate(this%upplim_destruct_metamorph%val) + deallocate(this%pc%val) + deallocate(this%mu%val) + deallocate(this%bsw_sf%val) + deallocate(this%hksat_sf%val) + deallocate(this%sucsat_sf%val) + deallocate(this%watsat_sf%val) + deallocate(this%d_max%val) + deallocate(this%frac_sat_soil_dsl_init%val) + deallocate(this%zlnd%val) + deallocate(this%snw_rds_min%val) + deallocate(this%precip_repartition_nonglc_all_rain_t_celsius%val) + deallocate(this%precip_repartition_nonglc_all_snow_t_celsius%val) + deallocate(this%precip_repartition_glc_all_rain_t_celsius%val) + deallocate(this%precip_repartition_glc_all_snow_t_celsius%val) end subroutine Clean diff --git a/src/main/atm2lndMod.F90 b/src/main/atm2lndMod.F90 index 5da4ff6333..cfd80264f4 100644 --- a/src/main/atm2lndMod.F90 +++ b/src/main/atm2lndMod.F90 @@ -25,7 +25,6 @@ module atm2lndMod use landunit_varcon, only : istice use WaterType , only : water_type use Wateratm2lndBulkType, only : wateratm2lndbulk_type - ! ! !PUBLIC TYPES: implicit none @@ -266,7 +265,7 @@ subroutine downscale_forcings(bounds, & end do - ! adjust hillslope precpitation before repartitioning rain/snow + ! adjust hillslope precipitation before repartitioning rain/snow if (use_hillslope .and. downscale_hillslope_meteorology) then call downscale_hillslope_solar(bounds, atm2lnd_inst, surfalb_inst) call downscale_hillslope_precipitation(bounds, topo_inst, atm2lnd_inst, wateratm2lndbulk_inst) @@ -359,15 +358,16 @@ subroutine partition_precip(bounds, atm2lnd_inst, wateratm2lndbulk_inst, eflx_sh if (atm2lnd_inst%params%repartition_rain_snow) then do c = bounds%begc, bounds%endc if (col%active(c)) then + g = col%gridcell(c) l = col%landunit(c) rain_orig = forc_rain_c(c) snow_orig = forc_snow_c(c) if (lun%itype(l) == istice) then - all_snow_t = atm2lnd_inst%params%precip_repartition_glc_all_snow_t - frac_rain_slope = atm2lnd_inst%params%precip_repartition_glc_frac_rain_slope + all_snow_t = atm2lnd_inst%params%precip_repartition_glc_all_snow_t%param_val(g) + frac_rain_slope = atm2lnd_inst%params%precip_repartition_glc_frac_rain_slope%param_val(g) else - all_snow_t = atm2lnd_inst%params%precip_repartition_nonglc_all_snow_t - frac_rain_slope = atm2lnd_inst%params%precip_repartition_nonglc_frac_rain_slope + all_snow_t = atm2lnd_inst%params%precip_repartition_nonglc_all_snow_t%param_val(g) + frac_rain_slope = atm2lnd_inst%params%precip_repartition_nonglc_frac_rain_slope%param_val(g) end if call repartition_rain_snow_one_col(& temperature = forc_t_c(c), & @@ -387,6 +387,7 @@ subroutine partition_precip(bounds, atm2lnd_inst, wateratm2lndbulk_inst, eflx_sh sens_heat_flux = eflx_sh_precip_conversion(c)) end if end do + end if end associate diff --git a/src/main/atm2lndType.F90 b/src/main/atm2lndType.F90 index 298ca4a41d..e2df497379 100644 --- a/src/main/atm2lndType.F90 +++ b/src/main/atm2lndType.F90 @@ -14,6 +14,10 @@ module atm2lndType use decompMod , only : bounds_type use abortutils , only : endrun use PatchType , only : patch + use DistParamType , only : distparam_class, distparams => distributed_parameters + + + use GridcellType , only : grc ! debug ! ! !PUBLIC TYPES: implicit none @@ -43,14 +47,19 @@ module atm2lndType ! Rain-snow ramp for glacier landunits ! frac_rain = (temp - all_snow_t) * frac_rain_slope ! (all_snow_t is in K) - real(r8) :: precip_repartition_glc_all_snow_t - real(r8) :: precip_repartition_glc_frac_rain_slope + class(distparam_class), pointer :: precip_repartition_glc_all_snow_t => NULL() + class(distparam_class), pointer :: precip_repartition_glc_frac_rain_slope => NULL() +! real(r8) :: precip_repartition_glc_all_snow_t +! real(r8) :: precip_repartition_glc_frac_rain_slope ! Rain-snow ramp for non-glacier landunits ! frac_rain = (temp - all_snow_t) * frac_rain_slope ! (all_snow_t is in K) - real(r8) :: precip_repartition_nonglc_all_snow_t - real(r8) :: precip_repartition_nonglc_frac_rain_slope + class(distparam_class), pointer :: precip_repartition_nonglc_all_snow_t => NULL() + class(distparam_class), pointer :: precip_repartition_nonglc_frac_rain_slope => NULL() + +! real(r8) :: precip_repartition_nonglc_all_snow_t +! real(r8) :: precip_repartition_nonglc_frac_rain_slope end type atm2lnd_params_type !---------------------------------------------------- @@ -143,9 +152,7 @@ module atm2lndType !----------------------------------------------------------------------- function atm2lnd_params_constructor(repartition_rain_snow, glcmec_downscale_longwave, & - lapse_rate, lapse_rate_longwave, longwave_downscaling_limit, & - precip_repartition_glc_all_snow_t, precip_repartition_glc_all_rain_t, & - precip_repartition_nonglc_all_snow_t, precip_repartition_nonglc_all_rain_t) & + lapse_rate, lapse_rate_longwave, longwave_downscaling_limit) & result(params) ! ! !DESCRIPTION: @@ -169,18 +176,9 @@ function atm2lnd_params_constructor(repartition_rain_snow, glcmec_downscale_long ! Must be present if glcmec_downscale_longwave is true; ignored otherwise real(r8), intent(in), optional :: longwave_downscaling_limit - ! End-points of the rain-snow ramp for glacier landunits (degrees C) - ! Must be present if repartition_rain_snow is true; ignored otherwise - real(r8), intent(in), optional :: precip_repartition_glc_all_snow_t - real(r8), intent(in), optional :: precip_repartition_glc_all_rain_t - - ! End-points of the rain-snow ramp for non-glacier landunits (degrees C) - ! Must be present if repartition_rain_snow is true; ignored otherwise - real(r8), intent(in), optional :: precip_repartition_nonglc_all_snow_t - real(r8), intent(in), optional :: precip_repartition_nonglc_all_rain_t ! ! !LOCAL VARIABLES: - + integer :: beg_index, end_index character(len=*), parameter :: subname = 'atm2lnd_params_constructor' !----------------------------------------------------------------------- @@ -212,70 +210,102 @@ function atm2lnd_params_constructor(repartition_rain_snow, glcmec_downscale_long params%longwave_downscaling_limit = nan end if - if (repartition_rain_snow) then - ! Make sure all of the repartitioning-related parameters are present + ! allocate derived type variables + allocate(params%precip_repartition_glc_all_snow_t) + allocate(params%precip_repartition_glc_frac_rain_slope) + allocate(params%precip_repartition_nonglc_all_snow_t) + allocate(params%precip_repartition_nonglc_frac_rain_slope) + + ! allocate derived type arrays based on distributed parameter arrays + if (.not. allocated(distparams%precip_repartition_glc_all_rain_t_celsius%val)) then + call endrun(subname //' ERROR: Must initialize precip_repartition_glc_all_rain_t_celsius') + endif + if (.not. allocated(distparams%precip_repartition_glc_all_snow_t_celsius%val)) then + call endrun(subname //' ERROR: Must initialize precip_repartition_glc_all_snow_t_celsius') + endif + if (.not. allocated(distparams%precip_repartition_nonglc_all_rain_t_celsius%val)) then + call endrun(subname //' ERROR: Must initialize precip_repartition_nonglc_all_rain_t_celsius') + endif + if (.not. allocated(distparams%precip_repartition_nonglc_all_snow_t_celsius%val)) then + call endrun(subname //' ERROR: Must initialize precip_repartition_nonglc_all_snow_t_celsius') + endif - if (.not. present(precip_repartition_glc_all_snow_t)) then - call endrun(subname // & - ' ERROR: For repartition_rain_snow true, precip_repartition_glc_all_snow_t must be provided') - end if - if (.not. present(precip_repartition_glc_all_rain_t)) then - call endrun(subname // & - ' ERROR: For repartition_rain_snow true, precip_repartition_glc_all_rain_t must be provided') - end if - if (.not. present(precip_repartition_nonglc_all_snow_t)) then - call endrun(subname // & - ' ERROR: For repartition_rain_snow true, precip_repartition_nonglc_all_snow_t must be provided') - end if - if (.not. present(precip_repartition_nonglc_all_rain_t)) then - call endrun(subname // & - ' ERROR: For repartition_rain_snow true, precip_repartition_nonglc_all_rain_t must be provided') - end if + ! set %is_distributed based on distparams values + params%precip_repartition_glc_all_snow_t%is_distributed = & + distparams%precip_repartition_glc_all_rain_t_celsius%is_distributed + params%precip_repartition_glc_frac_rain_slope%is_distributed = & + distparams%precip_repartition_glc_all_rain_t_celsius%is_distributed + params%precip_repartition_nonglc_all_snow_t%is_distributed = & + distparams%precip_repartition_nonglc_all_rain_t_celsius%is_distributed + params%precip_repartition_nonglc_frac_rain_slope%is_distributed = & + distparams%precip_repartition_nonglc_all_rain_t_celsius%is_distributed + + ! allocate glacier parameter arrays + beg_index = lbound(distparams%precip_repartition_glc_all_rain_t_celsius%val,1) + end_index = ubound(distparams%precip_repartition_glc_all_rain_t_celsius%val,1) + allocate(params%precip_repartition_glc_all_snow_t%val(beg_index:end_index)) + allocate(params%precip_repartition_glc_frac_rain_slope%val(beg_index:end_index)) + params%precip_repartition_glc_all_snow_t%val(:) = nan + params%precip_repartition_glc_frac_rain_slope%val(:) = nan + + ! allocate non-glacier parameter arrays + beg_index = lbound(distparams%precip_repartition_nonglc_all_rain_t_celsius%val,1) + end_index = ubound(distparams%precip_repartition_nonglc_all_rain_t_celsius%val,1) + allocate(params%precip_repartition_nonglc_all_snow_t%val(beg_index:end_index)) + allocate(params%precip_repartition_nonglc_frac_rain_slope%val(beg_index:end_index)) + params%precip_repartition_nonglc_all_snow_t%val(:) = nan + params%precip_repartition_nonglc_frac_rain_slope%val(:) = nan - ! Do some other error checking + if (repartition_rain_snow) then - if (precip_repartition_glc_all_rain_t <= precip_repartition_glc_all_snow_t) then - call endrun(subname // & - ' ERROR: Must have precip_repartition_glc_all_snow_t < precip_repartition_glc_all_rain_t') - end if + ! Do some other error checking - if (precip_repartition_nonglc_all_rain_t <= precip_repartition_nonglc_all_snow_t) then - call endrun(subname // & - ' ERROR: Must have precip_repartition_nonglc_all_snow_t < precip_repartition_nonglc_all_rain_t') - end if + !todo: should this be done on a gridcell basis, or just removed? + +!!$ if (precip_repartition_glc_all_rain_t <= precip_repartition_glc_all_snow_t) then +!!$ call endrun(subname // & +!!$ ' ERROR: Must have precip_repartition_glc_all_snow_t < precip_repartition_glc_all_rain_t') +!!$ end if +!!$ +!!$ if (precip_repartition_nonglc_all_rain_t <= precip_repartition_nonglc_all_snow_t) then +!!$ call endrun(subname // & +!!$ ' ERROR: Must have precip_repartition_nonglc_all_snow_t < precip_repartition_nonglc_all_rain_t') +!!$ end if ! Convert to the form of the parameters we want for the main code call compute_ramp_params( & - all_snow_t_c = precip_repartition_glc_all_snow_t, & - all_rain_t_c = precip_repartition_glc_all_rain_t, & - all_snow_t_k = params%precip_repartition_glc_all_snow_t, & - frac_rain_slope = params%precip_repartition_glc_frac_rain_slope) + all_snow_t_c = distparams%precip_repartition_glc_all_snow_t_celsius%val, & + all_rain_t_c = distparams%precip_repartition_glc_all_rain_t_celsius%val, & + all_snow_t_k = params%precip_repartition_glc_all_snow_t%val, & + frac_rain_slope = params%precip_repartition_glc_frac_rain_slope%val) call compute_ramp_params( & - all_snow_t_c = precip_repartition_nonglc_all_snow_t, & - all_rain_t_c = precip_repartition_nonglc_all_rain_t, & - all_snow_t_k = params%precip_repartition_nonglc_all_snow_t, & - frac_rain_slope = params%precip_repartition_nonglc_frac_rain_slope) - - else ! .not. repartition_rain_snow - params%precip_repartition_glc_all_snow_t = nan - params%precip_repartition_glc_frac_rain_slope = nan - params%precip_repartition_nonglc_all_snow_t = nan - params%precip_repartition_nonglc_frac_rain_slope = nan + all_snow_t_c = distparams%precip_repartition_nonglc_all_snow_t_celsius%val, & + all_rain_t_c = distparams%precip_repartition_nonglc_all_rain_t_celsius%val, & + all_snow_t_k = params%precip_repartition_nonglc_all_snow_t%val, & + frac_rain_slope = params%precip_repartition_nonglc_frac_rain_slope%val) + end if contains subroutine compute_ramp_params(all_snow_t_c, all_rain_t_c, & all_snow_t_k, frac_rain_slope) - real(r8), intent(in) :: all_snow_t_c ! Temperature at which precip falls entirely as rain (deg C) - real(r8), intent(in) :: all_rain_t_c ! Temperature at which precip falls entirely as snow (deg C) - real(r8), intent(out) :: all_snow_t_k ! Temperature at which precip falls entirely as snow (K) - real(r8), intent(out) :: frac_rain_slope ! Slope of the frac_rain vs. T relationship - - frac_rain_slope = 1._r8 / (all_rain_t_c - all_snow_t_c) - all_snow_t_k = all_snow_t_c + tfrz + real(r8), intent(in) :: all_snow_t_c(:) ! Temperature at which precip falls entirely as rain (deg C) + real(r8), intent(in) :: all_rain_t_c(:) ! Temperature at which precip falls entirely as snow (deg C) + real(r8), intent(out) :: all_snow_t_k(:) ! Temperature at which precip falls entirely as snow (K) + real(r8), intent(out) :: frac_rain_slope(:) ! Slope of the frac_rain vs. T relationship + integer :: g + integer :: beg_index, end_index + + beg_index = lbound(all_rain_t_c,1) + end_index = ubound(all_rain_t_c,1) + + do g = beg_index,end_index + frac_rain_slope(g) = 1._r8 / (all_rain_t_c(g) - all_snow_t_c(g)) + all_snow_t_k(g) = all_snow_t_c(g) + tfrz + enddo end subroutine compute_ramp_params end function atm2lnd_params_constructor @@ -356,10 +386,6 @@ subroutine ReadNamelist(this, NLFilename) real(r8) :: lapse_rate real(r8) :: lapse_rate_longwave real(r8) :: longwave_downscaling_limit - real(r8) :: precip_repartition_glc_all_snow_t - real(r8) :: precip_repartition_glc_all_rain_t - real(r8) :: precip_repartition_nonglc_all_snow_t - real(r8) :: precip_repartition_nonglc_all_rain_t integer :: ierr ! error code integer :: unitn ! unit for namelist file @@ -369,9 +395,7 @@ subroutine ReadNamelist(this, NLFilename) !----------------------------------------------------------------------- namelist /atm2lnd_inparm/ repartition_rain_snow, glcmec_downscale_longwave, & - lapse_rate, lapse_rate_longwave, longwave_downscaling_limit, & - precip_repartition_glc_all_snow_t, precip_repartition_glc_all_rain_t, & - precip_repartition_nonglc_all_snow_t, precip_repartition_nonglc_all_rain_t + lapse_rate, lapse_rate_longwave, longwave_downscaling_limit ! Initialize namelist variables to defaults repartition_rain_snow = .false. @@ -379,10 +403,6 @@ subroutine ReadNamelist(this, NLFilename) lapse_rate = nan lapse_rate_longwave = nan longwave_downscaling_limit = nan - precip_repartition_glc_all_snow_t = nan - precip_repartition_glc_all_rain_t = nan - precip_repartition_nonglc_all_snow_t = nan - precip_repartition_nonglc_all_rain_t = nan if (masterproc) then unitn = getavu() @@ -404,10 +424,6 @@ subroutine ReadNamelist(this, NLFilename) call shr_mpi_bcast(lapse_rate, mpicom) call shr_mpi_bcast(lapse_rate_longwave, mpicom) call shr_mpi_bcast(longwave_downscaling_limit, mpicom) - call shr_mpi_bcast(precip_repartition_glc_all_snow_t, mpicom) - call shr_mpi_bcast(precip_repartition_glc_all_rain_t, mpicom) - call shr_mpi_bcast(precip_repartition_nonglc_all_snow_t, mpicom) - call shr_mpi_bcast(precip_repartition_nonglc_all_rain_t, mpicom) if (masterproc) then write(iulog,*) ' ' @@ -421,12 +437,6 @@ subroutine ReadNamelist(this, NLFilename) write(iulog,*) 'lapse_rate_longwave = ', lapse_rate_longwave write(iulog,*) 'longwave_downscaling_limit = ', longwave_downscaling_limit end if - if (repartition_rain_snow) then - write(iulog,*) 'precip_repartition_glc_all_snow_t = ', precip_repartition_glc_all_snow_t - write(iulog,*) 'precip_repartition_glc_all_rain_t = ', precip_repartition_glc_all_rain_t - write(iulog,*) 'precip_repartition_nonglc_all_snow_t = ', precip_repartition_nonglc_all_snow_t - write(iulog,*) 'precip_repartition_nonglc_all_rain_t = ', precip_repartition_nonglc_all_rain_t - end if write(iulog,*) ' ' end if @@ -435,15 +445,10 @@ subroutine ReadNamelist(this, NLFilename) glcmec_downscale_longwave = glcmec_downscale_longwave, & lapse_rate = lapse_rate, & lapse_rate_longwave = lapse_rate_longwave, & - longwave_downscaling_limit = longwave_downscaling_limit, & - precip_repartition_glc_all_snow_t = precip_repartition_glc_all_snow_t, & - precip_repartition_glc_all_rain_t = precip_repartition_glc_all_rain_t, & - precip_repartition_nonglc_all_snow_t = precip_repartition_nonglc_all_snow_t, & - precip_repartition_nonglc_all_rain_t = precip_repartition_nonglc_all_rain_t) + longwave_downscaling_limit = longwave_downscaling_limit) end subroutine ReadNamelist - !------------------------------------------------------------------------ subroutine InitAllocate(this, bounds) ! diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90 index a703bfeb00..543adc4c27 100644 --- a/src/main/clm_initializeMod.F90 +++ b/src/main/clm_initializeMod.F90 @@ -25,7 +25,8 @@ module clm_initializeMod use LandunitType , only : lun ! instance use ColumnType , only : col ! instance use PatchType , only : patch ! instance - use DistParamType , only : distparams +! use DistParamType , only : distparams, distparamsstreams + use DistParamMod , only : InitDistributedParameters use reweightMod , only : reweight_wrapup use filterMod , only : allocFilters, filter, filter_inactive_and_active use CLMFatesInterfaceMod , only : CLMFatesGlobals1,CLMFatesGlobals2 @@ -350,8 +351,9 @@ subroutine initialize2(ni,nj, currtime) call readParameters(photosyns_inst) ! Read in spatially distributed parameters - call distparams%Init(bounds_proc) - call distparams%readDistributedParams(bounds_proc) +! call distparams%Init(bounds_proc) +!scs call distparams%readDistributedParams(bounds_proc) + call InitDistributedParameters(bounds_proc) ! Initialize time manager if (nsrest == nsrStartup) then diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90 index a4d7b59f87..1b2005cb3f 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -398,6 +398,12 @@ module clm_varctl integer, public :: CN_evergreen_phenology_opt = 0 integer, public :: carbon_resp_opt = 0 + !---------------------------------------------------------- + ! spatially distributed parameters switch + !---------------------------------------------------------- + + logical, public :: use_distributed_parameters = .false. ! true => use spatially distributed parameter stream data + !---------------------------------------------------------- ! prescribed soil moisture streams switch !---------------------------------------------------------- diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90 index 6e0be9055d..983cd231cd 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -46,7 +46,7 @@ module controlMod use SoilBiogeochemLittVertTranspMod , only: som_adv_flux, max_depth_cryoturb use SoilBiogeochemVerticalProfileMod , only: surfprof_exp use SoilBiogeochemNitrifDenitrifMod , only: no_frozen_nitrif_denitrif - use SoilHydrologyMod , only: soilHydReadNML, hillslope_hydrology_ReadNML + use SoilHydrologyMod , only: hillslope_hydrology_ReadNML use CNFireFactoryMod , only: CNFireReadNML use CanopyFluxesMod , only: CanopyFluxesReadNML use shr_drydep_mod , only: n_drydep @@ -616,7 +616,6 @@ subroutine control_init(dtime) call CNFUNReadNML( NLFilename ) end if - call soilHydReadNML( NLFilename ) if ( use_hillslope ) then call hillslope_hydrology_ReadNML( NLFilename ) endif diff --git a/src/main/initVerticalMod.F90 b/src/main/initVerticalMod.F90 index dd67ff8169..87a540e668 100644 --- a/src/main/initVerticalMod.F90 +++ b/src/main/initVerticalMod.F90 @@ -27,7 +27,7 @@ module initVerticalMod use LandunitType , only : lun use GridcellType , only : grc use ColumnType , only : col - use DistParamType , only : distparams + use DistParamType , only : distparams => distributed_parameters use glcBehaviorMod , only : glc_behavior_type use abortUtils , only : endrun use ncdio_pio @@ -651,12 +651,12 @@ subroutine initVertical(bounds, glc_behavior, thick_wall, thick_roof) do c = begc,endc ! microtopographic parameter, units are meters (try smooth function of slope) - slope0 = distparams%slopemax(c)**(1._r8/distparams%slopebeta(c)) + slope0 = distparams%slopemax%param_val(col%gridcell(c))**(1._r8/distparams%slopebeta%param_val(col%gridcell(c))) if (col%is_hillslope_column(c)) then - col%micro_sigma(c) = (atan(col%hill_slope(c)) + slope0)**(distparams%slopebeta(c)) + col%micro_sigma(c) = (atan(col%hill_slope(c)) + slope0)**(distparams%slopebeta%param_val(col%gridcell(c))) else - col%micro_sigma(c) = (col%topo_slope(c) + slope0)**(distparams%slopebeta(c)) + col%micro_sigma(c) = (col%topo_slope(c) + slope0)**(distparams%slopebeta%param_val(col%gridcell(c))) endif end do diff --git a/src/main/ncdio_pio.F90.in b/src/main/ncdio_pio.F90.in index 991823fb67..7f0a5931af 100644 --- a/src/main/ncdio_pio.F90.in +++ b/src/main/ncdio_pio.F90.in @@ -62,6 +62,7 @@ module ncdio_pio public :: ncd_inqdname ! inquire dimension name public :: ncd_inqdlen ! inquire dimension length public :: ncd_inqfdims ! inquire file dimnesions + public :: ncd_inqnvars ! inquire number of variables on file public :: ncd_defvar ! define variables public :: ncd_inqvid ! inquire variable id public :: ncd_inqvname ! inquire variable name @@ -678,7 +679,7 @@ contains !----------------------------------------------------------------------- vname = '' - status = PIO_inq_varname(ncid,vardesc,vname) + status = PIO_inq_varname(ncid,varid,vname) end subroutine ncd_inqvname @@ -897,6 +898,26 @@ contains end subroutine ncd_inqnatts + !----------------------------------------------------------------------- + subroutine ncd_inqnvars(ncid, nvariables) + ! + ! !DESCRIPTION: + ! Inquire number of variables + ! + ! !ARGUMENTS: + class(file_desc_t), intent(inout) :: ncid ! netcdf file id + integer , intent(out) :: nvariables ! number of variables + ! + ! !LOCAL VARIABLES: + integer :: status + + character(len=*), parameter :: subname = 'ncd_inqnvars' + !----------------------------------------------------------------------- + + status = PIO_inquire(ncid, nvariables=nvariables) + + end subroutine ncd_inqnvars + !----------------------------------------------------------------------- subroutine ncd_inqattname(ncid, varid, attnum, attname) ! From 16442eecb7fc3333e6c3f2b4599bde43ff6119dd Mon Sep 17 00:00:00 2001 From: Sean Swenson Date: Mon, 11 Aug 2025 08:28:33 -0600 Subject: [PATCH 6/6] remove distributed_parameter file --- bld/CLMBuildNamelist.pm | 1 - bld/namelist_files/namelist_defaults_ctsm.xml | 6 +----- bld/namelist_files/namelist_definition_ctsm.xml | 5 ----- 3 files changed, 1 insertion(+), 11 deletions(-) diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index ca79a93b04..7fd1bbac51 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -2255,7 +2255,6 @@ sub setup_logic_params_file { add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'paramfile', 'phys'=>$nl_flags->{'phys'}, 'use_flexibleCN'=>$nl_flags->{'use_flexibleCN'} ); - add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'distributed_paramfile'); } #------------------------------------------------------------------------------- diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index 651835dd0a..0eeeb73b74 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -527,13 +527,9 @@ attributes from the config_cache.xml file (with keys converted to upper-case). -lnd/clm2/paramdata/ctsm5.3.041.Nfix_params.v13.c250221_upplim250.nc +lnd/clm2/paramdata/ctsm5.3.041.Nfix_params.v13.c250221_upplim250_add_nml.nc lnd/clm2/paramdata/clm50_params.c250311.nc lnd/clm2/paramdata/clm45_params.c250311.nc - -lnd/clm2/paramdata/ctsm5.3_distributed_params_placeholder.nc -lnd/clm2/paramdata/ctsm5.3_distributed_params_placeholder.nc - diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 227e8e10a7..457a87293a 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -998,11 +998,6 @@ Full pathname datafile with plant function type (PFT) constants combined with constants for biogeochem modules - -Full pathname spatially distributed parameter file - - Full pathname datafile with fates parameters