diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 06ea82d99b..7fd1bbac51 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 # ########################################## @@ -4420,6 +4425,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) = @_; @@ -4658,19 +4682,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"); - } - } - } } #------------------------------------------------------------------------------- @@ -5206,7 +5217,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 50e3cf68ad..0eeeb73b74 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. @@ -541,10 +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 - @@ -2082,6 +2067,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 820975655d..457a87293a 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -1694,30 +1694,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. - - @@ -2055,6 +2031,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/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..d5e4acc3b4 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 => distributed_parameters 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 @@ -267,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 ! @@ -296,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**(-params_inst%e_ice*(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 9956a7dfb8..cdee4c7ec2 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 => distributed_parameters use SoilHydrologyType, only : soilhydrology_type - use SoilStateType, only : soilstate_type + use SoilStateType , only : soilstate_type use WaterFluxBulkType, only : waterfluxbulk_type implicit none @@ -50,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: @@ -175,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 ! ======================================================================== @@ -315,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 ! @@ -349,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*params_inst%fff*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*params_inst%fff*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 f6a5ff41ee..a791b7ffe1 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 => distributed_parameters use glcBehaviorMod , only : glc_behavior_type use landunit_varcon, only : istice use paramUtilMod , only : readNcdioScalar @@ -44,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 @@ -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(this%accum_factor * 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(this%accum_factor * 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 @@ -291,7 +293,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' @@ -301,10 +302,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) @@ -314,7 +311,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 @@ -376,30 +372,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) ! @@ -430,7 +402,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 @@ -440,7 +412,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: @@ -461,7 +432,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%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 578769d9ea..870d846e25 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 => 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 @@ -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%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 5f38830891..7289987fff 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 => distributed_parameters use SoilStateInitTimeConstMod, only: organic_max ! implicit none @@ -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%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 5a4aa50f6e..a13d971977 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 => distributed_parameters ! ! !PUBLIC TYPES: @@ -40,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 @@ -57,18 +57,8 @@ 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 integer, private :: head_gradient_method ! Method for calculating hillslope saturated head gradient @@ -85,7 +75,7 @@ module SoilHydrologyMod __FILE__ contains - + !----------------------------------------------------------------------- subroutine hillslope_hydrology_ReadNML(NLFilename) ! @@ -171,91 +161,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 ) - ! - ! !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) @@ -791,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, params_inst%aq_sp_yield_min) + 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 @@ -806,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, params_inst%aq_sp_yield_min) + 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) @@ -821,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, params_inst%aq_sp_yield_min) + 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) @@ -1088,7 +993,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%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 @@ -1121,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**(-params_inst%e_ice*(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 @@ -1197,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**(-params_inst%e_ice*(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 @@ -1242,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**(-params_inst%e_ice*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**(-params_inst%e_ice*(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 @@ -1271,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, params_inst%aq_sp_yield_min) + 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 @@ -1308,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, params_inst%aq_sp_yield_min) + 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) @@ -1916,7 +1821,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%param_val(col%gridcell(c)) & * sin(col%topo_slope(c) * (rpi/180._r8)) wtsub = 0._r8 @@ -1963,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,params_inst%aq_sp_yield_min) + 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 @@ -2177,6 +2082,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 +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**(-params_inst%e_ice*icefrac(c,j)) + ice_imped(c,j)=10._r8**(-distparams%e_ice%param_val(col%gridcell(c))*icefrac(c,j)) end do end do @@ -2243,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**(-params_inst%e_ice*(icefracsum/dzsum)) + ice_imped_col(c)=10._r8**(-distparams%e_ice%param_val(col%gridcell(c))*(icefracsum/dzsum)) enddo do fc = 1, num_hydrologyc @@ -2391,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) * baseflow_scalar & + 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))**(params_inst%n_baseflow) + (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)) @@ -2456,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,params_inst%aq_sp_yield_min) + 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)) @@ -2483,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,params_inst%aq_sp_yield_min) + 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) @@ -2822,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, params_inst%aq_sp_yield_min) + 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 a730417315..3f5311f8ca 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 => distributed_parameters use abortUtils , only : endrun use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) ! @@ -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%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) = params_inst%bsw_sf * ( (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) = params_inst%sucsat_sf * ( (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 @@ -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%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)) @@ -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%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) = params_inst%bsw_sf * ( (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) = params_inst%sucsat_sf * ( (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 @@ -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%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 85bcf42c5e..3c168cd8a2 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 => distributed_parameters ! 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%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* & @@ -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,8 +2123,10 @@ end subroutine compute_qcharge !#13 !----------------------------------------------------------------------- - subroutine IceImpedance(icefrac, imped) + subroutine IceImpedance(c, icefrac, imped) ! + ! USES + use ColumnType , only : col !DESCRIPTION ! compute soil suction potential ! @@ -2133,15 +2136,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%param_val(col%gridcell(c))*icefrac) end subroutine IceImpedance diff --git a/src/biogeophys/SurfaceResistanceMod.F90 b/src/biogeophys/SurfaceResistanceMod.F90 index 9b04327252..64d2bdba79 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 => distributed_parameters + use WaterStateBulkType, only : waterstatebulk_type + use WaterDiagnosticBulkType, only : waterdiagnosticbulk_type + use TemperatureType , only : temperature_type implicit none save private @@ -29,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__ @@ -165,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, & @@ -401,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) = 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%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, (params_inst%frac_sat_soil_dsl_init * 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 562c64cc18..e59471d514 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 => distributed_parameters use NumericsMod , only : truncate_small_values use InfiltrationExcessRunoffMod , only : infiltration_excess_runoff_type use EnergyFluxType , only : energyflux_type @@ -32,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) @@ -473,10 +446,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%param_val(col%gridcell(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%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 f91aaca761..b8e9ed4250 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 => distributed_parameters 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%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 new file mode 100644 index 0000000000..ea57cb956f --- /dev/null +++ b/src/cpl/share_esmf/DistParamStreamMod.F90 @@ -0,0 +1,364 @@ +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_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) + 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 diff --git a/src/main/DistParamType.F90 b/src/main/DistParamType.F90 new file mode 100644 index 0000000000..e17a89cd04 --- /dev/null +++ b/src/main/DistParamType.F90 @@ -0,0 +1,443 @@ +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 + 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_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 + 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 + class(distparam_class), pointer :: fff => NULL() ! Decay factor for fractional saturated area (1/m) + + ! initVerticalMod + 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 + 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 + class(distparam_class), pointer :: upplim_destruct_metamorph => NULL() ! Upper limit on destructive metamorphism compaction (kg/m3) + + ! SoilStateInitTimeConstMod + 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 + 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 + 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 + 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 + 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 :: Clean + + end type distributed_parameter_type + + type(distributed_parameter_type), public, target :: distributed_parameters + + public :: InitGlobalParameters + + !------------------------------------------------------------------------ + +contains + + !------------------------------------------------------------------------ + function param_val(this,g) + ! + ! !DESCRIPTION: + ! Return distributed parameter value if allocated, otherwise return scalar value + ! + ! !ARGUMENTS: + implicit none + class(distparam_class) :: this + real(r8) :: param_val + integer :: g + ! + if ( this%is_distributed )then + param_val = this%val(g) + else + param_val = this%val(1) + end if + end function param_val + + !------------------------------------------------------------------------ + 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' + !-------------------------------------------------------------------- + + if ( .not. this%is_distributed ) then + allocate(this%val(1)) + call readNcdioScalar(ncid, this%name, subname, fscalar_in) + this%val(1) = fscalar_in + endif + + end subroutine ReadScalarParameter + + !------------------------------------------------------------------------ + subroutine InitGlobalParameters(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 + type(bounds_type), intent(in) :: bounds + ! + ! !LOCAL VARIABLES: + logical :: readvar ! whether the variable was found + integer :: c, g + real(r8) :: fscalar_in ! read in - scalar - float + real(r8), pointer :: fparam_in(:) ! read in - 1D - float + 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 parameters' + end if + + !distributed_parameter_stream%Init must be called before this routine is called. + + ! 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 ReadScalarParameter(distributed_parameters%aq_sp_yield_min,ncid) + + ! Drainage power law exponent (unitless) + 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 ReadScalarParameter(distributed_parameters%perched_baseflow_scalar,ncid) + + ! Soil ice impedance factor (unitless) + call ReadScalarParameter(distributed_parameters%e_ice,ncid) + + !----------------------------------------------------------- + ! SaturatedExcessRunoff ! + !----------------------------------------------------------- + + ! Decay factor for fractional saturated area (1/m) + call ReadScalarParameter(distributed_parameters%fff,ncid) + + !----------------------------------------------------------- + ! initVertical ! + !----------------------------------------------------------- + + ! exponent for microtopography pdf sigma (unitless) + call ReadScalarParameter(distributed_parameters%slopebeta,ncid) + + ! max topographic slope for microtopography pdf sigma (unitless) + call ReadScalarParameter(distributed_parameters%slopemax,ncid) + + !----------------------------------------------------------- + ! SnowCoverFractionSwensonLawrence2012 ! + !----------------------------------------------------------- + + ! SCA shape parameter + call ReadScalarParameter(distributed_parameters%n_melt_coef,ncid) + + ! Accumulation constant for fractional snow covered area (unitless) + call ReadScalarParameter(distributed_parameters%accum_factor,ncid) + + !----------------------------------------------------------- + ! SnowHydrologyMod ! + !----------------------------------------------------------- + + ! Upper limit on destructive metamorphism compaction (kg/m3) + call ReadScalarParameter(distributed_parameters%upplim_destruct_metamorph,ncid) + + !----------------------------------------------------------- + ! SoilStateInitTimeConstMod ! + !----------------------------------------------------------- + + call ReadScalarParameter(distributed_parameters%bsw_sf,ncid) + + call ReadScalarParameter(distributed_parameters%hksat_sf,ncid) + + call ReadScalarParameter(distributed_parameters%sucsat_sf,ncid) + + call ReadScalarParameter(distributed_parameters%watsat_sf,ncid) + + !----------------------------------------------------------- + ! SurfaceResistanceMod ! + !----------------------------------------------------------- + + ! Dry surface layer parameter (mm) + call ReadScalarParameter(distributed_parameters%d_max,ncid) + + ! Fraction of saturated soil for moisture value at which DSL initiates (unitless) + call ReadScalarParameter(distributed_parameters%frac_sat_soil_dsl_init,ncid) + + !----------------------------------------------------------- + ! SurfaceWaterMod ! + !----------------------------------------------------------- + + ! Threshold probability for surface water (unitless) + call ReadScalarParameter(distributed_parameters%pc,ncid) + + ! Connectivity exponent for surface water (unitless) + call ReadScalarParameter(distributed_parameters%mu,ncid) + + !----------------------------------------------------------- + ! WaterDiagnosticBulkType ! + !----------------------------------------------------------- + + ! Momentum roughness length for soil, glacier, wetland (m) + 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 ! + !----------------------------------------------------------- + + ! 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 Init(this) + ! + ! !ARGUMENTS: + class(distributed_parameter_type) :: this + !------------------------------------------------------------------------ + + 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) + + this%aq_sp_yield_min%name = 'aq_sp_yield_min' + this%aq_sp_yield_min%is_distributed = .false. + + this%n_baseflow%name = 'n_baseflow' + this%n_baseflow%is_distributed = .false. + + this%perched_baseflow_scalar%name = 'perched_baseflow_scalar' + this%perched_baseflow_scalar%is_distributed = .false. + + this%e_ice%name = 'e_ice' + this%e_ice%is_distributed = .false. + + this%baseflow_scalar%name = 'baseflow_scalar' + this%baseflow_scalar%is_distributed = .false. + + this%fff%name = 'fff' + this%fff%is_distributed = .false. + + this%slopebeta%name = 'slopebeta' + this%slopebeta%is_distributed = .false. + + 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. + + this%watsat_sf%name = 'watsat_sf' + this%watsat_sf%is_distributed = .false. + + 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 + +end module Distparamtype 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 8c0b50230b..543adc4c27 100644 --- a/src/main/clm_initializeMod.F90 +++ b/src/main/clm_initializeMod.F90 @@ -25,6 +25,8 @@ module clm_initializeMod use LandunitType , only : lun ! instance use ColumnType , only : col ! instance use PatchType , only : patch ! instance +! 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 @@ -342,12 +344,16 @@ 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) +!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 41978ae695..1b2005cb3f 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 @@ -397,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 6d363a9a6e..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 @@ -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 @@ -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 @@ -755,6 +754,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 +1045,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..87a540e668 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 => distributed_parameters use glcBehaviorMod , only : glc_behavior_type use abortUtils , only : endrun use ncdio_pio @@ -45,8 +46,6 @@ 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 @@ -76,11 +75,6 @@ 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) @@ -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%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)**(params_inst%slopebeta) + 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)**(params_inst%slopebeta) + 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) ! diff --git a/src/main/readParamsMod.F90 b/src/main/readParamsMod.F90 index c11a741198..ffe0b6a123 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 @@ -49,20 +52,17 @@ 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 + ! ! !ARGUMENTS: type(photosyns_type) , intent(in) :: photosyns_inst @@ -126,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)