From 43241150ec1ed9ca20dff3b01f96d6bf1f477f38 Mon Sep 17 00:00:00 2001 From: Rich Fiorella Date: Sat, 30 Mar 2024 19:31:42 -0600 Subject: [PATCH 1/7] Add maxalt across overall simulation. Adds another variable of maxalt_ever and its index following existing logic for calculating the maximum active layer thickness (really thaw depth as calculated) to date in the simulation. --- .../elm/src/biogeophys/ActiveLayerMod.F90 | 42 +++++++++++-------- .../elm/src/biogeophys/CanopyStateType.F90 | 22 ++++++++++ 2 files changed, 46 insertions(+), 18 deletions(-) diff --git a/components/elm/src/biogeophys/ActiveLayerMod.F90 b/components/elm/src/biogeophys/ActiveLayerMod.F90 index 940fe5a3ba8d..e3513fd7a71e 100644 --- a/components/elm/src/biogeophys/ActiveLayerMod.F90 +++ b/components/elm/src/biogeophys/ActiveLayerMod.F90 @@ -1,5 +1,5 @@ module ActiveLayerMod - + !----------------------------------------------------------------------- ! !DESCRIPTION: ! Module holding routines for calculation of active layer dynamics @@ -10,9 +10,9 @@ module ActiveLayerMod use elm_varctl , only : iulog use TemperatureType , only : temperature_type use CanopyStateType , only : canopystate_type - use GridcellType , only : grc_pp + use GridcellType , only : grc_pp use ColumnType , only : col_pp - use ColumnDataType , only : col_es + use ColumnDataType , only : col_es ! implicit none save @@ -21,12 +21,12 @@ module ActiveLayerMod ! !PUBLIC MEMBER FUNCTIONS: public:: alt_calc !----------------------------------------------------------------------- - + contains !----------------------------------------------------------------------- subroutine alt_calc(num_soilc, filter_soilc, & - temperature_vars, canopystate_vars) + temperature_vars, canopystate_vars) ! ! !DESCRIPTION: ! define active layer thickness similarly to frost_table, except set as deepest thawed layer and define on nlevgrnd @@ -68,15 +68,17 @@ subroutine alt_calc(num_soilc, filter_soilc, & real(r8) :: t1, t2, z1, z2 ! temporary variables !----------------------------------------------------------------------- - associate( & - t_soisno => col_es%t_soisno , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) (-nlevsno+1:nlevgrnd) - - alt => canopystate_vars%alt_col , & ! Output: [real(r8) (:) ] current depth of thaw - altmax => canopystate_vars%altmax_col , & ! Output: [real(r8) (:) ] maximum annual depth of thaw - altmax_lastyear => canopystate_vars%altmax_lastyear_col , & ! Output: [real(r8) (:) ] prior year maximum annual depth of thaw - alt_indx => canopystate_vars%alt_indx_col , & ! Output: [integer (:) ] current depth of thaw - altmax_indx => canopystate_vars%altmax_indx_col , & ! Output: [integer (:) ] maximum annual depth of thaw - altmax_lastyear_indx => canopystate_vars%altmax_lastyear_indx_col & ! Output: [integer (:) ] prior year maximum annual depth of thaw + associate( & + t_soisno => col_es%t_soisno , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) (-nlevsno+1:nlevgrnd) + + alt => canopystate_vars%alt_col , & ! Output: [real(r8) (:) ] current depth of thaw + altmax => canopystate_vars%altmax_col , & ! Output: [real(r8) (:) ] maximum annual depth of thaw + altmax_lastyear => canopystate_vars%altmax_lastyear_col , & ! Output: [real(r8) (:) ] prior year maximum annual depth of thaw + altmax_ever => canopystate_vars%altmax_ever_col , & ! Output: [real(r8) (:) ] maximum thaw depth since initialization + alt_indx => canopystate_vars%alt_indx_col , & ! Output: [integer (:) ] current depth of thaw + altmax_indx => canopystate_vars%altmax_indx_col , & ! Output: [integer (:) ] maximum annual depth of thaw + altmax_lastyear_indx => canopystate_vars%altmax_lastyear_indx_col , & ! Output: [integer (:) ] prior year maximum annual depth of thaw + altmax_ever_indx => canopystate_vars%altmax_ever_indx_col & ! Output: [integer (:) ] maximum thaw depth since initialization ) ! on a set annual timestep, update annual maxima @@ -88,7 +90,7 @@ subroutine alt_calc(num_soilc, filter_soilc, & c = filter_soilc(fc) g = col_pp%gridcell(c) if ( grc_pp%lat(g) > 0. ) then - + altmax_lastyear(c) = altmax(c) altmax_lastyear_indx(c) = altmax_indx(c) altmax(c) = 0. @@ -100,7 +102,7 @@ subroutine alt_calc(num_soilc, filter_soilc, & do fc = 1,num_soilc c = filter_soilc(fc) g = col_pp%gridcell(c) - if ( grc_pp%lat(g) <= 0. ) then + if ( grc_pp%lat(g) <= 0. ) then altmax_lastyear(c) = altmax(c) altmax_lastyear_indx(c) = altmax_indx(c) altmax(c) = 0. @@ -150,11 +152,15 @@ subroutine alt_calc(num_soilc, filter_soilc, & altmax(c) = alt(c) altmax_indx(c) = alt_indx(c) endif + if (alt(c) > altmax_ever(c)) then + altmax_ever(c) = alt(c) + altmax_ever_indx(c) = alt_indx(c) + endif end do - end associate + end associate end subroutine alt_calc - + end module ActiveLayerMod diff --git a/components/elm/src/biogeophys/CanopyStateType.F90 b/components/elm/src/biogeophys/CanopyStateType.F90 index 2f4598782c73..098859139de7 100644 --- a/components/elm/src/biogeophys/CanopyStateType.F90 +++ b/components/elm/src/biogeophys/CanopyStateType.F90 @@ -62,8 +62,10 @@ module CanopyStateType integer , pointer :: alt_indx_col (:) ! col current depth of thaw real(r8) , pointer :: altmax_col (:) ! col maximum annual depth of thaw real(r8) , pointer :: altmax_lastyear_col (:) ! col prior year maximum annual depth of thaw + real(r8) , pointer :: altmax_ever_col (:) ! col maximum thaw depth from beginning of simulation integer , pointer :: altmax_indx_col (:) ! col maximum annual depth of thaw integer , pointer :: altmax_lastyear_indx_col (:) ! col prior year maximum annual depth of thaw + integer , pointer :: altmax_ever_indx_col (:) ! col maximum thaw depth from beginning of simulation real(r8) , pointer :: dewmx_patch (:) ! patch maximum allowed dew [mm] @@ -142,9 +144,11 @@ subroutine InitAllocate(this, bounds ) allocate(this%alt_col (begc:endc)) ; this%alt_col (:) = spval; allocate(this%altmax_col (begc:endc)) ; this%altmax_col (:) = spval allocate(this%altmax_lastyear_col (begc:endc)) ; this%altmax_lastyear_col (:) = spval + allocate(this%altmax_ever_col (begc:endc)) ; this%altmax_ever_col (:) = spval allocate(this%alt_indx_col (begc:endc)) ; this%alt_indx_col (:) = huge(1) allocate(this%altmax_indx_col (begc:endc)) ; this%altmax_indx_col (:) = huge(1) allocate(this%altmax_lastyear_indx_col (begc:endc)) ; this%altmax_lastyear_indx_col (:) = huge(1) + allocate(this%altmax_ever_indx_col (begc:endc)) ; this%altmax_ever_indx_col (:) = huge(1) allocate(this%dewmx_patch (begp:endp)) ; this%dewmx_patch (:) = spval allocate(this%dleaf_patch (begp:endp)) ; this%dleaf_patch (:) = spval @@ -277,6 +281,11 @@ subroutine InitHistory(this, bounds) call hist_addfld1d (fname='ALTMAX_LASTYEAR', units='m', & avgflag='A', long_name='maximum prior year active layer thickness', & ptr_col=this%altmax_lastyear_col) + + this%altmax_ever_col(begc:endc) = spval + call hist_addfld1d (fname='ALTMAX_EVER', units='m', & + avgflag='A', long_name='maximum ever active layer thickness throughout simulation', & + ptr_col=this%altmax_ever_col) end if ! Allow active layer fields to be optionally output even if not running CN @@ -296,6 +305,11 @@ subroutine InitHistory(this, bounds) call hist_addfld1d (fname='ALTMAX_LASTYEAR', units='m', & avgflag='A', long_name='maximum prior year active layer thickness', & ptr_col=this%altmax_lastyear_col, default='inactive') + + this%altmax_ever_col(begc:endc) = spval + call hist_addfld1d (fname='ALTMAX_EVER', units='m', & + avgflag='A', long_name='maximum ever active layer thickness throughout simulation', & + ptr_col=this%altmax_ever_col, default='inactive') end if ! Accumulated fields @@ -502,9 +516,11 @@ subroutine InitCold(this, bounds) this%alt_col(c) = 0._r8 this%altmax_col(c) = 0._r8 this%altmax_lastyear_col(c) = 0._r8 + this%altmax_ever_col(c) = 0._r8 this%alt_indx_col(c) = 0 this%altmax_indx_col(c) = 0 this%altmax_lastyear_indx_col(c) = 0 + this%altmax_ever_indx_col(c) = 0 end if end do @@ -572,12 +588,18 @@ subroutine Restart(this, bounds, ncid, flag) call restartvar(ncid=ncid, flag=flag, varname='altmax_lastyear', xtype=ncd_double, & dim1name='column', long_name='', units='', & interpinic_flag='interp', readvar=readvar, data=this%altmax_lastyear_col) + call restartvar(ncid=ncid, flag=flag, varname='altmax_ever', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=this%altmax_ever_col) call restartvar(ncid=ncid, flag=flag, varname='altmax_indx', xtype=ncd_int, & dim1name='column', long_name='', units='', & interpinic_flag='interp', readvar=readvar, data=this%altmax_indx_col) call restartvar(ncid=ncid, flag=flag, varname='altmax_lastyear_indx', xtype=ncd_int, & dim1name='column', long_name='', units='', & interpinic_flag='interp', readvar=readvar, data=this%altmax_lastyear_indx_col) + call restartvar(ncid=ncid, flag=flag, varname='altmax_ever_indx', xtype=ncd_int, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=this%altmax_ever_indx_col) end if if ( use_hydrstress ) then call restartvar(ncid=ncid, flag=flag, varname='vegwp', xtype=ncd_double, & From 6911dbf470d960e620e2adb0abda7e0a5d929a29 Mon Sep 17 00:00:00 2001 From: Rich Fiorella Date: Fri, 29 Mar 2024 08:16:37 -0600 Subject: [PATCH 2/7] Update ci-compile-test.yml --- .github/workflows/ci-compile-test.yml | 54 +++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) create mode 100644 .github/workflows/ci-compile-test.yml diff --git a/.github/workflows/ci-compile-test.yml b/.github/workflows/ci-compile-test.yml new file mode 100644 index 000000000000..241e2f7d9481 --- /dev/null +++ b/.github/workflows/ci-compile-test.yml @@ -0,0 +1,54 @@ +name: CI +on: + push: + branches: + - '**' +jobs: + container-test-job: + runs-on: ubuntu-latest + container: + image: rfiorella/e3sm-dev:latest + steps: + - name: Check out the NGEE-Arctic-E3SM repo + uses: actions/checkout@v4 + with: + ref: ${{ github.ref }} + submodules: true + - name: Setup cime and input datafiles + run: | + mkdir ~/.cime + cp /home/e3smuser/.cime/* /github/home/.cime/ + git clone https://github.com/rfiorella/pt-e3sm-inputdata /home/e3smuser/inputdata + cd /home/e3smuser/inputdata/lnd/clm2/firedata + tar -zxvf clmforc.Li_2012_hdm_0.5x0.5_AVHRR_simyr1850-2010_c130401.nc.tar.gz + cd /home/e3smuser/inputdata/atm/datm7/NASA_LIS/ + tar -zxvf clmforc.Li_2012_climo1995-2011.T62.lnfm_c130327.nc.tar.gz + tar -zxvf clmforc.Li_2012_climo1995-2011.T62.lnfm_Total_c140423.nc.tar.gz + ls -v /home/e3smuser/inputdata/atm/datm7/NASA_LIS/ + - name: create new case and setup + run: | + cd cime/scripts && ./create_newcase --mach docker --res ELM_USRDAT --compset ICB1850CNPRDCTCBC --case ~/output/ci_test \ + && cd ~/output/ci_test + ./xmlchange MOSART_MODE=NULL,DOUT_S=FALSE,DIN_LOC_ROOT=/home/e3smuser/inputdata + ./xmlchange DIN_LOC_ROOT_CLMFORC=/home/e3smuser/inputdata/atm/datm7 + ./xmlchange ELM_USRDAT_NAME=1x1pt_AK-TLG + ./xmlchange ATM_DOMAIN_PATH=/home/e3smuser/inputdata/share/domains/domain.clm + ./xmlchange LND_DOMAIN_PATH=/home/e3smuser/inputdata/share/domains/domain.clm + ./xmlchange ATM_DOMAIN_FILE=domain.lnd.1x1pt_teller-GRID_navy.nc + ./xmlchange LND_DOMAIN_FILE=domain.lnd.1x1pt_teller-GRID_navy.nc + ./xmlchange PIO_TYPENAME=netcdf,NTASKS=1,NTHRDS=1 + ./xmlchange STOP_OPTION=nyears,STOP_N=5 + ./xmlchange BATCH_SYSTEM=none + echo "&clm_inparm" >> user_nl_elm + echo " do_budgets = .false." >> user_nl_elm + echo " use_nofire = .true." >> user_nl_elm + echo " metdata_type = 'gswp3'" >> user_nl_elm + echo " metdata_bypass = '/home/e3smuser/inputdata/atm/datm7/atm_forcing.datm7.GSWP3.0.5d.v2.c180716_NGEE-Grid/cpl_bypass_teller-Grid'" >> user_nl_elm + echo " fsurdat = '/home/e3smuser/inputdata/lnd/clm2/surfdata_map/surfdata_1x1pt_teller-GRID_simyr1850_c360x720_c171002.nc'" >> user_nl_elm + echo " stream_fldfilename_popdens = '/home/e3smuser/inputdata/lnd/clm2/firedata/clmforc.Li_2012_hdm_0.5x0.5_AVHRR_simyr1850-2010_c130401.nc'" >> user_nl_elm + ./case.setup + # add coupler bypass flag to cmake macros. + echo 'string(APPEND CPPDEFS " -DCPL_BYPASS")' >> cmake_macros/universal.cmake + ./case.build + ./case.submit --no-batch + From 0f6496a2b8e3b293bccfcf9fa9bbac0ae52248ac Mon Sep 17 00:00:00 2001 From: Rich Fiorella Date: Wed, 3 Apr 2024 17:35:16 -0600 Subject: [PATCH 3/7] wip - correct initialization of altmax_ever Changes needed to get altmax_ever values to initialize correctly. Fixes an issue where altmax_ever would be >30 meters regardless of location/climate due to carrying forward the initial conditions. Now includes a spinup_state flag so that it is set as 0 during spinup, and then increases based on thaw depths over non-spinup cases. --- components/elm/src/biogeophys/ActiveLayerMod.F90 | 11 ++++++++--- components/elm/src/biogeophys/CanopyStateType.F90 | 2 +- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/components/elm/src/biogeophys/ActiveLayerMod.F90 b/components/elm/src/biogeophys/ActiveLayerMod.F90 index e3513fd7a71e..9756ee4aefce 100644 --- a/components/elm/src/biogeophys/ActiveLayerMod.F90 +++ b/components/elm/src/biogeophys/ActiveLayerMod.F90 @@ -7,7 +7,7 @@ module ActiveLayerMod ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 use shr_const_mod , only : SHR_CONST_TKFRZ - use elm_varctl , only : iulog + use elm_varctl , only : iulog, spinup_state use TemperatureType , only : temperature_type use CanopyStateType , only : canopystate_type use GridcellType , only : grc_pp @@ -153,8 +153,13 @@ subroutine alt_calc(num_soilc, filter_soilc, & altmax_indx(c) = alt_indx(c) endif if (alt(c) > altmax_ever(c)) then - altmax_ever(c) = alt(c) - altmax_ever_indx(c) = alt_indx(c) + if (spinup_state .eq. 0) then !overwrite if in spinup + altmax_ever(c) = alt(c) + altmax_ever_indx(c) = alt_indx(c) + else + altmax_ever(c) = 0._r8 + altmax_ever_indx(c) = 0 + endif endif end do diff --git a/components/elm/src/biogeophys/CanopyStateType.F90 b/components/elm/src/biogeophys/CanopyStateType.F90 index 098859139de7..3ce62408e999 100644 --- a/components/elm/src/biogeophys/CanopyStateType.F90 +++ b/components/elm/src/biogeophys/CanopyStateType.F90 @@ -148,7 +148,7 @@ subroutine InitAllocate(this, bounds ) allocate(this%alt_indx_col (begc:endc)) ; this%alt_indx_col (:) = huge(1) allocate(this%altmax_indx_col (begc:endc)) ; this%altmax_indx_col (:) = huge(1) allocate(this%altmax_lastyear_indx_col (begc:endc)) ; this%altmax_lastyear_indx_col (:) = huge(1) - allocate(this%altmax_ever_indx_col (begc:endc)) ; this%altmax_ever_indx_col (:) = huge(1) + allocate(this%altmax_ever_indx_col (begc:endc)) ; this%altmax_ever_indx_col (:) = huge(1) allocate(this%dewmx_patch (begp:endp)) ; this%dewmx_patch (:) = spval allocate(this%dleaf_patch (begp:endp)) ; this%dleaf_patch (:) = spval From 038dc0f399786115b9261b4e3291b2437fc3e0b9 Mon Sep 17 00:00:00 2001 From: Rich Fiorella Date: Mon, 15 Apr 2024 12:22:36 -0600 Subject: [PATCH 4/7] Add excess ground ice as a column variable Adds excess ground ice to the col_ws / column water state data structure as a varaible needed for NGEE-Arctic Phase 3 IM1 (improved inundation fraction). Not yet connected to any of the new polygonal ground subtypes of the vegetated columns. --- .../elm/src/data_types/ColumnDataType.F90 | 158 ++++++++++-------- components/elm/src/main/surfrdMod.F90 | 2 +- 2 files changed, 86 insertions(+), 74 deletions(-) diff --git a/components/elm/src/data_types/ColumnDataType.F90 b/components/elm/src/data_types/ColumnDataType.F90 index 32ebf37af051..66f167fbb4c0 100644 --- a/components/elm/src/data_types/ColumnDataType.F90 +++ b/components/elm/src/data_types/ColumnDataType.F90 @@ -47,7 +47,7 @@ module ColumnDataType use CNDecompCascadeConType , only : decomp_cascade_con use ColumnType , only : col_pp use LandunitType , only : lun_pp - use timeInfoMod , only : nstep_mod + use timeInfoMod , only : nstep_mod ! ! !PUBLIC TYPES: implicit none @@ -142,7 +142,8 @@ module ColumnDataType real(r8), pointer :: snow_persistence (:) => null() ! length of time that ground has had non-zero snow thickness (sec) real(r8), pointer :: snw_rds_top (:) => null() ! snow grain radius (top layer) (m^-6, microns) logical , pointer :: do_capsnow (:) => null() ! true => do snow capping - real(r8), pointer :: h2osoi_tend_tsl_col(:) => null() ! col moisture tendency due to vertical movement at topmost layer (m3/m3/s) + real(r8), pointer :: h2osoi_tend_tsl_col(:) => null() ! col moisture tendency due to vertical movement at topmost layer (m3/m3/s) + real(r8), pointer :: excess_ice (:) => null() ! NGEE-Arctic: tracking excess ground ice ! Area fractions real(r8), pointer :: frac_sno (:) => null() ! fraction of ground covered by snow (0 to 1) real(r8), pointer :: frac_sno_eff (:) => null() ! fraction of ground covered by snow (0 to 1) @@ -559,7 +560,7 @@ module ColumnDataType real(r8), pointer :: phr_vr (:,:) => null() ! potential hr (not N-limited) (gC/m3/s) real(r8), pointer :: fphr (:,:) => null() ! fraction of potential heterotrophic respiration real(r8), pointer :: som_c_leached (:) => null() ! total SOM C loss from vertical transport (gC/m^2/s) - real(r8), pointer :: som_c_runoff (:) => null() + real(r8), pointer :: som_c_runoff (:) => null() ! phenology: litterfall and crop fluxes real(r8), pointer :: phenology_c_to_litr_met_c (:,:) => null() ! C fluxes associated with phenology (litterfall and crop) to litter metabolic pool (gC/m3/s) real(r8), pointer :: phenology_c_to_litr_cel_c (:,:) => null() ! C fluxes associated with phenology (litterfall and crop) to litter cellulose pool (gC/m3/s) @@ -775,10 +776,10 @@ module ColumnDataType ! leaching fluxes real(r8), pointer :: smin_no3_leached_vr (:,:) => null() ! vertically-resolved soil mineral NO3 loss to leaching (gN/m3/s) real(r8), pointer :: smin_no3_leached (:) => null() ! soil mineral NO3 pool loss to leaching (gN/m2/s) - real(r8), pointer :: smin_nh4_leached (:) => null() + real(r8), pointer :: smin_nh4_leached (:) => null() real(r8), pointer :: smin_no3_runoff_vr (:,:) => null() ! vertically-resolved rate of mineral NO3 loss with runoff (gN/m3/s) real(r8), pointer :: smin_no3_runoff (:) => null() ! soil mineral NO3 pool loss to runoff (gN/m2/s) - real(r8), pointer :: smin_nh4_runoff (:) => null() + real(r8), pointer :: smin_nh4_runoff (:) => null() real(r8), pointer :: nh3_soi_flx (:) => null() ! nitrification /denitrification diagnostic quantities real(r8), pointer :: smin_no3_massdens_vr (:,:) => null() ! (ugN / g soil) soil nitrate concentration @@ -875,7 +876,7 @@ module ColumnDataType real(r8), pointer :: nh3_stores (:) => null() ! NH3 emission from manure storage, (gN/m2/s) real(r8), pointer :: nh3_grz (:) => null() ! NH3 emission from manure on pastures, (gN/m2/s) real(r8), pointer :: nh3_manure_app (:) => null() ! NH3 emission from manure applied on crops and grasslands, (gN/m - real(r8), pointer :: nh3_fert (:) => null() ! NH3 emission from fertilizers applied on crops and grasslands, + real(r8), pointer :: nh3_fert (:) => null() ! NH3 emission from fertilizers applied on crops and grasslands, real(r8), pointer :: nh3_otherfert (:) => null() ! NH3 emission from non-urea fertilizers applied on crops and gra real(r8), pointer :: nh3_total (:) => null() ! Total NH3 emission from agriculture real(r8), pointer :: manure_no3_to_soil (:) => null() ! Nitrification flux from manure (gN/m2/s) @@ -888,7 +889,7 @@ module ColumnDataType real(r8), pointer :: fan_totnin (:) => null() ! Total output N from FAN pools, gN/m2/s real(r8), pointer :: manure_n_to_sminn (:) => null() ! Manure N from FAN pools to soil mineral pools, gN/m2/s real(r8), pointer :: synthfert_n_to_sminn (:) => null() ! Fertilizer N from FAN pools to soil mineral pools, gN/m2/s - real(r8), pointer :: manure_n_total (:) => null() ! Total manure N produced, gN/m2/s + real(r8), pointer :: manure_n_total (:) => null() ! Total manure N produced, gN/m2/s contains procedure, public :: Init => col_nf_init @@ -941,7 +942,7 @@ module ColumnDataType real(r8), pointer :: actual_immob_p (:) => null() ! vert-int (diagnostic) actual P immobilization (gP/m2/s) real(r8), pointer :: sminp_to_plant_vr (:,:) => null() ! vertically-resolved plant uptake of soil mineral P (gP/m3/s) real(r8), pointer :: sminp_to_plant (:) => null() ! vert-int (diagnostic) plant uptake of soil mineral P (gP/m2/s) - real(r8), pointer :: net_mineralization_p_vr (:,:) => null() + real(r8), pointer :: net_mineralization_p_vr (:,:) => null() real(r8), pointer :: supplement_to_sminp_vr (:,:) =>null() ! vertically-resolved supplemental P supply (gP/m3/s) real(r8), pointer :: supplement_to_sminp (:) =>null() ! vert-int (diagnostic) supplemental P supply (gP/m2/s) real(r8), pointer :: gross_pmin_vr (:,:) =>null() ! vertically-resolved gross rate of P mineralization (gP/m3/s) @@ -1383,10 +1384,10 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ allocate(this%h2osoi_liq (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_liq (:,:) = spval allocate(this%h2osoi_ice (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_ice (:,:) = spval allocate(this%h2osoi_vol (begc:endc, 1:nlevgrnd)) ; this%h2osoi_vol (:,:) = spval - allocate(this%h2osfc (begc:endc)) ; this%h2osfc (:) = spval - allocate(this%h2ocan (begc:endc)) ; this%h2ocan (:) = spval + allocate(this%h2osfc (begc:endc)) ; this%h2osfc (:) = spval + allocate(this%h2ocan (begc:endc)) ; this%h2ocan (:) = spval allocate(this%wslake_col (begc:endc)) ; this%wslake_col (:) = spval - allocate(this%total_plant_stored_h2o(begc:endc)) ; this%total_plant_stored_h2o(:)= spval + allocate(this%total_plant_stored_h2o(begc:endc)) ; this%total_plant_stored_h2o(:)= spval allocate(this%h2osoi_liqvol (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_liqvol (:,:) = spval allocate(this%h2osoi_icevol (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_icevol (:,:) = spval allocate(this%h2osoi_liq_old (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_liq_old (:,:) = spval @@ -1416,6 +1417,7 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ if (use_fan) then allocate(this%h2osoi_tend_tsl_col(begc:endc)) ; this%h2osoi_tend_tsl_col(:) = spval end if + allocate(this%excess_ice (begc:endc)) ; this%excess_ice (:) = spval allocate(this%snw_rds_top (begc:endc)) ; this%snw_rds_top (:) = spval allocate(this%do_capsnow (begc:endc)) allocate(this%frac_sno (begc:endc)) ; this%frac_sno (:) = spval @@ -1562,7 +1564,12 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ ptr_col=this%h2osoi_tend_tsl_col, l2g_scale_type='veg', & default='inactive') end if - + + this%excess_ice(begc:endc) = spval + call hist_addfld1d (fname='EXCESS_ICE', units = 'kg/m2', & + avgflag='A', long_name='Excess ground ice', & + ptr_col=this$excess_ice, l2g_scale_type='veg') ! <- RPF: should this be natveg? + this%frac_sno(begc:endc) = spval call hist_addfld1d (fname='FSNO', units='1', & avgflag='A', long_name='fraction of ground covered by snow', & @@ -1584,9 +1591,9 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ call hist_addfld1d (fname='FH2OSFC', units='1', & avgflag='A', long_name='fraction of ground covered by surface water', & ptr_col=this%frac_h2osfc) - + this%frac_h2osfc_act(begc:endc) = spval - + if (use_cn) then this%wf(begc:endc) = spval call hist_addfld1d (fname='WF', units='proportion', & @@ -1819,12 +1826,12 @@ subroutine col_ws_restart(this, bounds, ncid, flag, watsat_input) this%h2osfc(bounds%begc:bounds%endc) = 0.0_r8 end if - if(do_budgets) then + if(do_budgets) then call restartvar(ncid=ncid, flag=flag, varname='ENDWB', xtype=ncd_double, & dim1name='column', & long_name='water balance at end of timestep', units='kg/m2', & interpinic_flag='interp', readvar=readvar, data=this%endwb) - end if + end if call restartvar(ncid=ncid, flag=flag, varname='H2OSOI_LIQ', xtype=ncd_double, & dim1name='column', dim2name='levtot', switchdim=.true., & @@ -1899,6 +1906,11 @@ subroutine col_ws_restart(this, bounds, ncid, flag, watsat_input) this%snow_persistence(:) = 0.0_r8 end if + call restartvar(ncid=ncid, flag=flag, varname='EXCESS_ICE', xtype=ncd_double, & + dim1name='column', & + long_name='excess ground ice', units='kg/m2', & + interpinic_flag='interp', readvar=readvar, data=this%excess_ice) + call restartvar(ncid=ncid, flag=flag, varname='frac_sno', xtype=ncd_double, & dim1name='column', & long_name='fraction of ground covered by snow (0 to 1)',units='1',& @@ -2055,8 +2067,8 @@ subroutine col_cs_init(this, begc, endc, carbon_type, ratio, c12_carbonstate_var allocate(this%totcolc (begc:endc)) ; this%totcolc (:) = spval allocate(this%totblgc (begc:endc)) ; this%totblgc (:) = spval allocate(this%totvegc_abg (begc:endc)) ; this%totvegc_abg (:) = spval - allocate(this%begcb (begc:endc)) ; this%begcb (:) = spval - allocate(this%endcb (begc:endc)) ; this%endcb (:) = spval + allocate(this%begcb (begc:endc)) ; this%begcb (:) = spval + allocate(this%endcb (begc:endc)) ; this%endcb (:) = spval allocate(this%errcb (begc:endc)) ; this%errcb (:) = spval allocate(this%totpftc_beg (begc:endc)) ; this%totpftc_beg (:) = spval allocate(this%cwdc_beg (begc:endc)) ; this%cwdc_beg (:) = spval @@ -2092,7 +2104,7 @@ subroutine col_cs_init(this, begc, endc, carbon_type, ratio, c12_carbonstate_var ! Do not define history variables for CWD when fates is active if( decomp_cascade_con%is_cwd(l) .and. use_fates ) cycle - + if ( nlevdecomp_full > 1 ) then data2dptr => this%decomp_cpools_vr(:,:,l) fieldname = trim(decomp_cascade_con%decomp_pool_name_history(l))//'C_vr' @@ -3267,18 +3279,18 @@ subroutine col_ns_init(this, begc, endc, col_cs) allocate(this%tan_f1 (begc:endc)) ; this%tan_f1 (:) = spval allocate(this%tan_f2 (begc:endc)) ; this%tan_f2 (:) = spval allocate(this%tan_f3 (begc:endc)) ; this%tan_f3 (:) = spval - allocate(this%tan_f4 (begc:endc)) ; this%tan_f4 (:) = spval - allocate(this%fert_u1 (begc:endc)) ; this%fert_u1 (:) = spval - allocate(this%fert_u2 (begc:endc)) ; this%fert_u2 (:) = spval - allocate(this%manure_u_grz (begc:endc)) ; this%manure_u_grz (:) = spval - allocate(this%manure_a_grz (begc:endc)) ; this%manure_a_grz (:) = spval - allocate(this%manure_r_grz (begc:endc)) ; this%manure_r_grz (:) = spval - allocate(this%manure_u_app (begc:endc)) ; this%manure_u_app (:) = spval - allocate(this%manure_a_app (begc:endc)) ; this%manure_a_app (:) = spval - allocate(this%manure_r_app (begc:endc)) ; this%manure_r_app (:) = spval - allocate(this%manure_n_stored (begc:endc)) ; this%manure_n_stored (:) = spval - allocate(this%manure_tan_stored (begc:endc)) ; this%manure_tan_stored (:) = spval - allocate(this%fan_grz_fract (begc:endc)) ; this%fan_grz_fract (:) = spval + allocate(this%tan_f4 (begc:endc)) ; this%tan_f4 (:) = spval + allocate(this%fert_u1 (begc:endc)) ; this%fert_u1 (:) = spval + allocate(this%fert_u2 (begc:endc)) ; this%fert_u2 (:) = spval + allocate(this%manure_u_grz (begc:endc)) ; this%manure_u_grz (:) = spval + allocate(this%manure_a_grz (begc:endc)) ; this%manure_a_grz (:) = spval + allocate(this%manure_r_grz (begc:endc)) ; this%manure_r_grz (:) = spval + allocate(this%manure_u_app (begc:endc)) ; this%manure_u_app (:) = spval + allocate(this%manure_a_app (begc:endc)) ; this%manure_a_app (:) = spval + allocate(this%manure_r_app (begc:endc)) ; this%manure_r_app (:) = spval + allocate(this%manure_n_stored (begc:endc)) ; this%manure_n_stored (:) = spval + allocate(this%manure_tan_stored (begc:endc)) ; this%manure_tan_stored (:) = spval + allocate(this%fan_grz_fract (begc:endc)) ; this%fan_grz_fract (:) = spval end if allocate(this%fan_totn (begc:endc)) ; this%fan_totn (:) = spval allocate(this%totpftn_beg (begc:endc)) ; this%totpftn_beg (:) = spval @@ -3316,7 +3328,7 @@ subroutine col_ns_init(this, begc, endc, col_cs) do l = 1, ndecomp_pools if( decomp_cascade_con%is_cwd(l) .and. use_fates ) cycle - + if ( nlevdecomp_full > 1 ) then data2dptr => this%decomp_npools_vr(:,:,l) fieldname = trim(decomp_cascade_con%decomp_pool_name_history(l))//'N_vr' @@ -5092,7 +5104,7 @@ subroutine col_ps_restart ( this, bounds, ncid, flag, cnstate_vars) ! assume soil below 50 cm has the same p pool concentration ! divide 0.5m when convert p pools from g/m2 to g/m3 ! assume p pools evenly distributed at dif layers - if (nu_com .eq. 'RD') then + if (nu_com .eq. 'RD') then smax_c = smax(isoilorder(c)) ks_sorption_c = ks_sorption(isoilorder(c)) this%solutionp_vr(c,j) = (cnstate_vars%labp_col(c)/0.5_r8*ks_sorption_c)/& @@ -6149,12 +6161,12 @@ subroutine col_cf_init(this, begc, endc, carbon_type) call hist_addfld1d (fname='NEP', units='gC/m^2/s', & avgflag='A', long_name='net ecosystem production, excludes fire, landuse, and harvest flux, positive for sink', & ptr_col=this%nep) - + this%nbp(begc:endc) = spval call hist_addfld1d (fname='NBP', units='gC/m^2/s', & avgflag='A', long_name='net biome production, includes fire, landuse, and harvest flux, positive for sink', & ptr_col=this%nbp) - + this%nee(begc:endc) = spval call hist_addfld1d (fname='NEE', units='gC/m^2/s', & avgflag='A', long_name='net ecosystem exchange of carbon, includes fire, landuse,'& @@ -6177,13 +6189,13 @@ subroutine col_cf_init(this, begc, endc, carbon_type) call hist_addfld1d (fname='CWDC_HR', units='gC/m^2/s', & avgflag='A', long_name='coarse woody debris C heterotrophic respiration', & ptr_col=this%cwdc_hr) - + this%cwdc_loss(begc:endc) = spval call hist_addfld1d (fname='CWDC_LOSS', units='gC/m^2/s', & avgflag='A', long_name='coarse woody debris C loss', & ptr_col=this%cwdc_loss) end if - + this%lithr(begc:endc) = spval call hist_addfld1d (fname='LITTERC_HR', units='gC/m^2/s', & avgflag='A', long_name='litter C heterotrophic respiration', & @@ -6585,20 +6597,20 @@ subroutine col_cf_init(this, begc, endc, carbon_type) avgflag='A', long_name='annual sum of column-level NPP', & ptr_col=this%annsum_npp, default='inactive') - if(.not.use_fates)then + if(.not.use_fates)then ! C4MIP output variable, plant carbon flux to cwd (a part of fVegLitter) this%plant_c_to_cwdc(begc:endc) = spval call hist_addfld1d (fname='VEGC_TO_CWDC', units='gC/m^2/s', & avgflag='A', long_name='plant carbon flux to cwd', & ptr_col=this%plant_c_to_cwdc, default='inactive') - + ! C4MIP output variable, plant phosphorus flux to cwd (a part of fVegLitter) this%plant_p_to_cwdp(begc:endc) = spval call hist_addfld1d (fname='VEGP_TO_CWDP', units='gP/m^2/s', & avgflag='A', long_name='plant phosphorus flux to cwd', & ptr_col=this%plant_p_to_cwdp, default='inactive') end if - + ! end of C12 block else if ( carbon_type == 'c13' ) then @@ -7712,7 +7724,7 @@ subroutine col_cf_setvalues ( this, num_column, filter_column, value_column) this%f_co2_soil_vr(i,j) = value_column end do end do - + do k = 1, ndecomp_pools do j = 1, nlevdecomp_full do fi = 1,num_column @@ -7722,7 +7734,7 @@ subroutine col_cf_setvalues ( this, num_column, filter_column, value_column) end do end do ! pflotran - if(nstep_mod == 0 .or. is_first_restart_step()) then + if(nstep_mod == 0 .or. is_first_restart_step()) then do k = 1, ndecomp_pools do j = 1, nlevdecomp_full do fi = 1,num_column @@ -7736,7 +7748,7 @@ subroutine col_cf_setvalues ( this, num_column, filter_column, value_column) this%f_co2_soil(i) = value_column this%externalc_to_decomp_delta(i) = value_column end do - end if + end if do k = 1, ndecomp_pools do fi = 1,num_column @@ -7744,7 +7756,7 @@ subroutine col_cf_setvalues ( this, num_column, filter_column, value_column) this%decomp_cpools_yield(i,k) = value_column !if ero_ccycle this%m_decomp_cpools_to_fire(i,k) = value_column end do - end do + end do do fi = 1,num_column i = filter_column(fi) @@ -7762,7 +7774,7 @@ subroutine col_cf_setvalues ( this, num_column, filter_column, value_column) this%coutputs(i) = value_column this%cwdc_hr(i) = value_column this%litterc_loss(i) = value_column - + this%nee(i) = value_column ! Zero p2c column fluxes @@ -7776,8 +7788,8 @@ subroutine col_cf_setvalues ( this, num_column, filter_column, value_column) this%wood_harvestc(i) = value_column this%hrv_xsmrpool_to_atm(i) = value_column end do - - if(use_crop) then + + if(use_crop) then do fi = 1,num_column i = filter_column(fi) this%somc_fire(i) = value_column @@ -7795,16 +7807,16 @@ subroutine col_cf_setvalues ( this, num_column, filter_column, value_column) this%somc_erode(i) = value_column this%somc_deposit(i) = value_column this%somc_yield(i) = value_column - enddo - end if + enddo + end if do k = 1, ndecomp_pools do fi = 1,num_column i = filter_column(fi) this%decomp_cpools_leached(i,k) = value_column this%decomp_cpools_erode(i,k) = value_column this%decomp_cpools_deposit(i,k) = value_column - end do - end do + end do + end do do k = 1, ndecomp_pools do j = 1, nlevdecomp_full do fi = 1,num_column @@ -7847,7 +7859,7 @@ subroutine col_cf_zero_forfates_veg(this, bounds, num_soilc, filter_soilc) this%hrv_xsmrpool_to_atm(c) = 0._r8 end do - + do k = 1, ndecomp_pools do j = 1, nlevdecomp_full do fc = 1,num_soilc @@ -8276,7 +8288,7 @@ subroutine col_nf_init(this, begc, endc) allocate(this%actual_immob_no3 (begc:endc)) ; this%actual_immob_no3 (:) = spval allocate(this%actual_immob_nh4 (begc:endc)) ; this%actual_immob_nh4 (:) = spval allocate(this%smin_no3_to_plant (begc:endc)) ; this%smin_no3_to_plant (:) = spval - allocate(this%smin_nh4_to_plant (begc:endc)) ; this%smin_nh4_to_plant (:) = spval + allocate(this%smin_nh4_to_plant (begc:endc)) ; this%smin_nh4_to_plant (:) = spval allocate(this%plant_to_litter_nflux (begc:endc)) ; this%plant_to_litter_nflux (:) = 0._r8 allocate(this%plant_to_cwd_nflux (begc:endc)) ; this%plant_to_cwd_nflux (:) = spval ! C4MIP output variable @@ -8334,7 +8346,7 @@ subroutine col_nf_init(this, begc, endc) allocate(this%fan_totnin (begc:endc)) ; this%fan_totnin (:) = spval allocate(this%fan_totnout (begc:endc)) ; this%fan_totnout (:) = spval end if - + !----------------------------------------------------------------------- ! initialize history fields for select members of col_nf !----------------------------------------------------------------------- @@ -8362,13 +8374,13 @@ subroutine col_nf_init(this, begc, endc) ptr_col=this%nfix_to_sminn) ! C4MIP output variable, plant nitrogen flux to cwd (a part of fVegLitter) - if(.not.use_fates)then + if(.not.use_fates)then this%plant_n_to_cwdn(begc:endc) = spval call hist_addfld1d (fname='VEGN_TO_CWDN', units='gN/m^2/s', & avgflag='A', long_name='plant nitrogen flux to cwd', & ptr_col=this%plant_n_to_cwdn, default='inactive') end if - + do k = 1, ndecomp_pools if ( decomp_cascade_con%is_litter(k) .or. decomp_cascade_con%is_cwd(k) ) then this%m_decomp_npools_to_fire(begc:endc,k) = spval @@ -9331,7 +9343,7 @@ subroutine col_nf_setvalues ( this, num_column, filter_column, value_column) end do end do - if( use_pflotran .and. pf_cmode) then + if( use_pflotran .and. pf_cmode) then do j = 1, nlevdecomp_full do fi = 1,num_column @@ -9339,13 +9351,13 @@ subroutine col_nf_setvalues ( this, num_column, filter_column, value_column) ! pflotran this%plant_ndemand_vr(i,j) = value_column !use_elm_interface.and.use_pflotran .and. pf_cmode this%f_ngas_decomp_vr(i,j) = value_column ! "" - this%f_ngas_nitri_vr(i,j) = value_column ! "" + this%f_ngas_nitri_vr(i,j) = value_column ! "" this%f_ngas_denit_vr(i,j) = value_column this%f_n2o_soil_vr(i,j) = value_column this%f_n2_soil_vr(i,j) = value_column - end do - end do - end if + end do + end do + end if do fi = 1,num_column i = filter_column(fi) @@ -9361,7 +9373,7 @@ subroutine col_nf_setvalues ( this, num_column, filter_column, value_column) this%smin_no3_to_plant(i) = value_column this%fire_nloss(i) = value_column this%som_n_leached(i) = value_column - + this%hrv_deadstemn_to_prod10n(i) = value_column this%hrv_deadstemn_to_prod100n(i) = value_column this%hrv_cropn_to_prod1n(i) = value_column @@ -9466,7 +9478,7 @@ subroutine col_nf_setvalues ( this, num_column, filter_column, value_column) ! pflotran !------------------------------------------------------------------------ - if(nstep_mod == 0 .or. is_first_restart_step()) then + if(nstep_mod == 0 .or. is_first_restart_step()) then do k = 1, ndecomp_pools do j = 1, nlevdecomp_full do fi = 1,num_column @@ -9494,7 +9506,7 @@ subroutine col_nf_setvalues ( this, num_column, filter_column, value_column) this%nh4_net_transport_vr(i,j) = value_column end do end do - + do fi = 1,num_column i = filter_column(fi) ! only initializing in the first time-step @@ -9502,15 +9514,15 @@ subroutine col_nf_setvalues ( this, num_column, filter_column, value_column) end do end if - if(use_crop) then + if(use_crop) then do j = 1, nlevdecomp_full do fi = 1,num_column i = filter_column(fi) this%f_nit_vr(i,j) = value_column this%f_denit_vr(i,j) = value_column end do - end do - end if + end do + end if end subroutine col_nf_setvalues !----------------------------------------------------------------------- @@ -10904,7 +10916,7 @@ subroutine col_pf_setvalues ( this, num_column, filter_column, value_column) this%plant_pdemand_vr(i,j) = value_column this%adsorb_to_labilep_vr(i,j) = value_column this%desorb_to_solutionp_vr(i,j) = value_column - + end do end do @@ -10995,7 +11007,7 @@ subroutine col_pf_setvalues ( this, num_column, filter_column, value_column) end do end do end do - + do k = 1, ndecomp_pools do j = 1, nlevdecomp_full do fi = 1,num_column @@ -11008,7 +11020,7 @@ subroutine col_pf_setvalues ( this, num_column, filter_column, value_column) end do ! pflotran - if(nstep_mod == 0 .or. is_first_restart_step() ) then + if(nstep_mod == 0 .or. is_first_restart_step() ) then do k = 1, ndecomp_pools do j = 1, nlevdecomp_full do fi = 1,num_column @@ -11030,7 +11042,7 @@ subroutine col_pf_setvalues ( this, num_column, filter_column, value_column) this%externalp_to_decomp_delta(i) = value_column this%sminp_net_transport_delta(i) = value_column end do - end if + end if end subroutine col_pf_setvalues @@ -11211,7 +11223,7 @@ subroutine col_pf_summary(this, bounds, num_soilc, filter_soilc) end if ! vertically integrate column-level fire P losses - if(.not. use_fates) then + if(.not. use_fates) then do k = 1, ndecomp_pools do j = 1, nlevdecomp do fc = 1,num_soilc @@ -11222,7 +11234,7 @@ subroutine col_pf_summary(this, bounds, num_soilc, filter_soilc) end do end do end do - end if + end if ! vertically integrate column-level P erosion flux if (ero_ccycle) then diff --git a/components/elm/src/main/surfrdMod.F90 b/components/elm/src/main/surfrdMod.F90 index 9cbcd4b2b8e0..e35f8320b30c 100755 --- a/components/elm/src/main/surfrdMod.F90 +++ b/components/elm/src/main/surfrdMod.F90 @@ -35,7 +35,7 @@ module surfrdMod ! ! !PUBLIC MEMBER FUNCTIONS: public :: surfrd_get_globmask ! Reads global land mask (needed for setting domain decomp) - public :: surfrd_get_grid ! Read grid/ladnfrac data into domain (after domain decomp) + public :: surfrd_get_grid ! Read grid/landfrac data into domain (after domain decomp) public :: surfrd_get_topo ! Read grid topography into domain (after domain decomp) public :: surfrd_get_data ! Read surface dataset and determine subgrid weights public :: surfrd_get_grid_conn ! Reads grid connectivity information from domain file From ca687d29c3deda394bbb37d4b212fd00596d704e Mon Sep 17 00:00:00 2001 From: Rich Fiorella Date: Mon, 15 Apr 2024 12:38:22 -0600 Subject: [PATCH 5/7] (rebase drop in future) add some more debugging tools to IM1 --- .github/workflows/ci-compile-test.yml | 8 +- components/elm/src/.vscode/settings.json | 3 + .../elm/src/biogeophys/CanopyHydrologyMod.F90 | 235 +++++++++--------- .../elm/src/biogeophys/WaterStateType.F90 | 138 +++++----- .../elm/src/data_types/ColumnDataType.F90 | 3 +- 5 files changed, 198 insertions(+), 189 deletions(-) create mode 100644 components/elm/src/.vscode/settings.json diff --git a/.github/workflows/ci-compile-test.yml b/.github/workflows/ci-compile-test.yml index 241e2f7d9481..0ad8b77c12fb 100644 --- a/.github/workflows/ci-compile-test.yml +++ b/.github/workflows/ci-compile-test.yml @@ -16,7 +16,7 @@ jobs: submodules: true - name: Setup cime and input datafiles run: | - mkdir ~/.cime + mkdir ~/.cime cp /home/e3smuser/.cime/* /github/home/.cime/ git clone https://github.com/rfiorella/pt-e3sm-inputdata /home/e3smuser/inputdata cd /home/e3smuser/inputdata/lnd/clm2/firedata @@ -49,6 +49,8 @@ jobs: ./case.setup # add coupler bypass flag to cmake macros. echo 'string(APPEND CPPDEFS " -DCPL_BYPASS")' >> cmake_macros/universal.cmake - ./case.build + ./case.build || cat /home/e3smuser/output/ci_test/bld/e3sm.bldlog.* ./case.submit --no-batch - + cat /home/e3smuser/output/ci_test/run/e3sm.log.* + cat /home/e3smuser/output/ci_test/run/lnd.log.* + diff --git a/components/elm/src/.vscode/settings.json b/components/elm/src/.vscode/settings.json new file mode 100644 index 000000000000..cad7657dfa54 --- /dev/null +++ b/components/elm/src/.vscode/settings.json @@ -0,0 +1,3 @@ +{ + "cmake.configureOnOpen": false +} \ No newline at end of file diff --git a/components/elm/src/biogeophys/CanopyHydrologyMod.F90 b/components/elm/src/biogeophys/CanopyHydrologyMod.F90 index 7726d03105e6..13c83ed010e5 100755 --- a/components/elm/src/biogeophys/CanopyHydrologyMod.F90 +++ b/components/elm/src/biogeophys/CanopyHydrologyMod.F90 @@ -12,7 +12,7 @@ module CanopyHydrologyMod ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 use shr_log_mod , only : errMsg => shr_log_errMsg - use shr_infnan_mod , only : isnan => shr_infnan_isnan + use shr_infnan_mod , only : isnan => shr_infnan_isnan use shr_sys_mod , only : shr_sys_flush use decompMod , only : bounds_type use abortutils , only : endrun @@ -22,10 +22,10 @@ module CanopyHydrologyMod use AerosolType , only : aerosol_type use CanopyStateType , only : canopystate_type use TopounitDataType , only : top_as, top_af ! Atmospheric state and flux variables - use ColumnType , only : col_pp - use ColumnDataType , only : col_es, col_ws, col_wf + use ColumnType , only : col_pp + use ColumnDataType , only : col_es, col_ws, col_wf use VegetationType , only : veg_pp - use VegetationDataType, only : veg_ws, veg_wf + use VegetationDataType, only : veg_ws, veg_wf use elm_varcon , only : snw_rds_min use pftvarcon , only : irrigated use GridcellType , only : grc_pp @@ -41,10 +41,10 @@ module CanopyHydrologyMod ! ! !PRIVATE MEMBER FUNCTIONS: private :: FracWet ! Determine fraction of vegetated surface that is wet - private :: FracH2oSfc ! Determine fraction of land surfaces which are submerged + private :: FracH2oSfc ! Determine fraction of land surfaces which are submerged ! ! !PRIVATE DATA MEMBERS: - integer :: oldfflag=0 ! use old fsno parameterization (N&Y07) + integer :: oldfflag=0 ! use old fsno parameterization (N&Y07) !----------------------------------------------------------------------- contains @@ -73,7 +73,7 @@ subroutine CanopyHydrology_readnl( NLFilename ) namelist / elm_canopyhydrology_inparm / oldfflag ! ---------------------------------------------------------------------- - ! Read namelist from standard input. + ! Read namelist from standard input. ! ---------------------------------------------------------------------- if ( masterproc )then @@ -128,7 +128,7 @@ subroutine CanopyHydrology(bounds, & use SnowHydrologyMod , only : NewSnowBulkDensity ! ! !ARGUMENTS: - type(bounds_type) , intent(in) :: bounds + type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_nolakec ! number of column non-lake points in column filter integer , intent(in) :: filter_nolakec(:) ! column filter for non-lake points integer , intent(in) :: num_nolakep ! number of pft non-lake points in pft filter @@ -171,93 +171,93 @@ subroutine CanopyHydrology(bounds, & real(r8) :: newsnow(bounds%begc:bounds%endc) real(r8) :: snowmelt(bounds%begc:bounds%endc) integer :: j - + !--------------------------------------------------- ! initializing variables used to adjust irrigation on local processer - real(r8) :: qflx_irrig_grid(bounds%begg:bounds%endg) ! irrigation at grid level [mm/s] + real(r8) :: qflx_irrig_grid(bounds%begg:bounds%endg) ! irrigation at grid level [mm/s] integer :: irrigated_ppg(bounds%begg:bounds%endg) ! irrigated pft per grid integer :: gg !----------------------------------------------------------------------- - associate( & + associate( & pgridcell => veg_pp%gridcell , & ! Input: [integer (:) ] pft's gridcell ptopounit => veg_pp%topounit , & ! Input: [integer (:) ] pft's topounit - plandunit => veg_pp%landunit , & ! Input: [integer (:) ] pft's landunit - pcolumn => veg_pp%column , & ! Input: [integer (:) ] pft's column + plandunit => veg_pp%landunit , & ! Input: [integer (:) ] pft's landunit + pcolumn => veg_pp%column , & ! Input: [integer (:) ] pft's column pgwgt => veg_pp%wtgcell , & ! Input: [integer (:) ] pft's weight in gridcell - ltype => lun_pp%itype , & ! Input: [integer (:) ] landunit type - urbpoi => lun_pp%urbpoi , & ! Input: [logical (:) ] true => landunit is an urban point - cgridcell => col_pp%gridcell , & ! Input: [integer (:) ] columns's gridcell - clandunit => col_pp%landunit , & ! Input: [integer (:) ] columns's landunit - ctype => col_pp%itype , & ! Input: [integer (:) ] column type - pfti => col_pp%pfti , & ! Input: [integer (:) ] column's beginning pft index - npfts => col_pp%npfts , & ! Input: [integer (:) ] number of patches in column - snl => col_pp%snl , & ! Input: [integer (:) ] number of snow layers - n_melt => col_pp%n_melt , & ! Input: [real(r8) (:) ] SCA shape parameter - zi => col_pp%zi , & ! Output: [real(r8) (:,:) ] interface level below a "z" level (m) - dz => col_pp%dz , & ! Output: [real(r8) (:,:) ] layer depth (m) - z => col_pp%z , & ! Output: [real(r8) (:,:) ] layer thickness (m) - - forc_rain => top_af%rain , & ! Input: [real(r8) (:) ] rain rate (kg H2O/m**2/s, or mm liquid H2O/s) - forc_snow => top_af%snow , & ! Input: [real(r8) (:) ] snow rate (kg H2O/m**2/s, or mm liquid H2O/s) - forc_t => top_as%tbot , & ! Input: [real(r8) (:) ] atmospheric temperature (Kelvin) + ltype => lun_pp%itype , & ! Input: [integer (:) ] landunit type + urbpoi => lun_pp%urbpoi , & ! Input: [logical (:) ] true => landunit is an urban point + cgridcell => col_pp%gridcell , & ! Input: [integer (:) ] columns's gridcell + clandunit => col_pp%landunit , & ! Input: [integer (:) ] columns's landunit + ctype => col_pp%itype , & ! Input: [integer (:) ] column type + pfti => col_pp%pfti , & ! Input: [integer (:) ] column's beginning pft index + npfts => col_pp%npfts , & ! Input: [integer (:) ] number of patches in column + snl => col_pp%snl , & ! Input: [integer (:) ] number of snow layers + n_melt => col_pp%n_melt , & ! Input: [real(r8) (:) ] SCA shape parameter + zi => col_pp%zi , & ! Output: [real(r8) (:,:) ] interface level below a "z" level (m) + dz => col_pp%dz , & ! Output: [real(r8) (:,:) ] layer depth (m) + z => col_pp%z , & ! Output: [real(r8) (:,:) ] layer thickness (m) + + forc_rain => top_af%rain , & ! Input: [real(r8) (:) ] rain rate (kg H2O/m**2/s, or mm liquid H2O/s) + forc_snow => top_af%snow , & ! Input: [real(r8) (:) ] snow rate (kg H2O/m**2/s, or mm liquid H2O/s) + forc_t => top_as%tbot , & ! Input: [real(r8) (:) ] atmospheric temperature (Kelvin) forc_wind => top_as%windbot , & ! Input: [real(r8) (:) ] atmospheric wind speed (m/s) - qflx_floodg => atm2lnd_vars%forc_flood_grc , & ! Input: [real(r8) (:) ] gridcell flux of flood water from RTM + qflx_floodg => atm2lnd_vars%forc_flood_grc , & ! Input: [real(r8) (:) ] gridcell flux of flood water from RTM - dewmx => canopystate_vars%dewmx_patch , & ! Input: [real(r8) (:) ] Maximum allowed dew [mm] + dewmx => canopystate_vars%dewmx_patch , & ! Input: [real(r8) (:) ] Maximum allowed dew [mm] frac_veg_nosno => canopystate_vars%frac_veg_nosno_patch , & ! Input: [integer (:) ] fraction of vegetation not covered by snow (0 OR 1) [-] elai => canopystate_vars%elai_patch , & ! Input: [real(r8) (:) ] one-sided leaf area index with burying by snow esai => canopystate_vars%esai_patch , & ! Input: [real(r8) (:) ] one-sided stem area index with burying by snow - t_grnd => col_es%t_grnd , & ! Input: [real(r8) (:) ] ground temperature (Kelvin) - t_soisno => col_es%t_soisno , & ! Output: [real(r8) (:,:) ] soil temperature (Kelvin) + t_grnd => col_es%t_grnd , & ! Input: [real(r8) (:) ] ground temperature (Kelvin) + t_soisno => col_es%t_soisno , & ! Output: [real(r8) (:,:) ] soil temperature (Kelvin) - do_capsnow => col_ws%do_capsnow , & ! Output: [logical (:) ] true => do snow capping - h2ocan => veg_ws%h2ocan , & ! Output: [real(r8) (:) ] total canopy water (mm H2O) - h2osfc => col_ws%h2osfc , & ! Output: [real(r8) (:) ] surface water (mm) - h2osno => col_ws%h2osno , & ! Output: [real(r8) (:) ] snow water (mm H2O) - snow_depth => col_ws%snow_depth , & ! Output: [real(r8) (:) ] snow height (m) - int_snow => col_ws%int_snow , & ! Output: [real(r8) (:) ] integrated snowfall [mm] + do_capsnow => col_ws%do_capsnow , & ! Output: [logical (:) ] true => do snow capping + h2ocan => veg_ws%h2ocan , & ! Output: [real(r8) (:) ] total canopy water (mm H2O) + h2osfc => col_ws%h2osfc , & ! Output: [real(r8) (:) ] surface water (mm) + h2osno => col_ws%h2osno , & ! Output: [real(r8) (:) ] snow water (mm H2O) + snow_depth => col_ws%snow_depth , & ! Output: [real(r8) (:) ] snow height (m) + int_snow => col_ws%int_snow , & ! Output: [real(r8) (:) ] integrated snowfall [mm] frac_sno_eff => col_ws%frac_sno_eff , & ! Output: [real(r8) (:) ] eff. fraction of ground covered by snow (0 to 1) frac_sno => col_ws%frac_sno , & ! Output: [real(r8) (:) ] fraction of ground covered by snow (0 to 1) frac_h2osfc => col_ws%frac_h2osfc , & ! Output: [real(r8) (:) ] fraction of ground covered by surface water (0 to 1) frac_iceold => col_ws%frac_iceold , & ! Output: [real(r8) (:,:) ] fraction of ice relative to the tot water - h2osoi_ice => col_ws%h2osoi_ice , & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) - h2osoi_liq => col_ws%h2osoi_liq , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) - swe_old => col_ws%swe_old , & ! Output: [real(r8) (:,:) ] snow water before update + h2osoi_ice => col_ws%h2osoi_ice , & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) + h2osoi_liq => col_ws%h2osoi_liq , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) + swe_old => col_ws%swe_old , & ! Output: [real(r8) (:,:) ] snow water before update irrig_rate => veg_wf%irrig_rate , & ! Input: [real(r8) (:) ] current irrigation rate (applied if n_irrig_steps_left > 0) [mm/s] n_irrig_steps_left => veg_wf%n_irrig_steps_left , & ! Output: [integer (:) ] number of time steps for which we still need to irrigate today - qflx_floodc => col_wf%qflx_floodc , & ! Output: [real(r8) (:) ] column flux of flood water from RTM - qflx_snow_melt => col_wf%qflx_snow_melt , & ! Output: [real(r8) (:) ] snow melt from previous time step - qflx_snow_h2osfc => col_wf%qflx_snow_h2osfc , & ! Output: [real(r8) (:) ] snow falling on surface water (mm/s) + qflx_floodc => col_wf%qflx_floodc , & ! Output: [real(r8) (:) ] column flux of flood water from RTM + qflx_snow_melt => col_wf%qflx_snow_melt , & ! Output: [real(r8) (:) ] snow melt from previous time step + qflx_snow_h2osfc => col_wf%qflx_snow_h2osfc , & ! Output: [real(r8) (:) ] snow falling on surface water (mm/s) qflx_snwcp_liq => veg_wf%qflx_snwcp_liq , & ! Output: [real(r8) (:) ] excess rainfall due to snow capping (mm H2O /s) [+] qflx_snwcp_ice => veg_wf%qflx_snwcp_ice , & ! Output: [real(r8) (:) ] excess snowfall due to snow capping (mm H2O /s) [+] qflx_snow_grnd_col => col_wf%qflx_snow_grnd , & ! Output: [real(r8) (:) ] snow on ground after interception (mm H2O/s) [+] qflx_snow_grnd_patch => veg_wf%qflx_snow_grnd , & ! Output: [real(r8) (:) ] snow on ground after interception (mm H2O/s) [+] - qflx_prec_intr => veg_wf%qflx_prec_intr , & ! Output: [real(r8) (:) ] interception of precipitation [mm/s] + qflx_prec_intr => veg_wf%qflx_prec_intr , & ! Output: [real(r8) (:) ] interception of precipitation [mm/s] qflx_prec_grnd => veg_wf%qflx_prec_grnd , & ! Output: [real(r8) (:) ] water onto ground including canopy runoff [kg/(m2 s)] qflx_rain_grnd => veg_wf%qflx_rain_grnd , & ! Output: [real(r8) (:) ] rain on ground after interception (mm H2O/s) [+] - qflx_dirct_rain => veg_wf%qflx_dirct_rain , & ! Output: [real(r8) (:) ] direct rain throughfall on ground (mm H2O/s) + qflx_dirct_rain => veg_wf%qflx_dirct_rain , & ! Output: [real(r8) (:) ] direct rain throughfall on ground (mm H2O/s) qflx_leafdrip => veg_wf%qflx_leafdrip , & ! Output: [real(r8) (:) ] leap rain drip on ground (mm H2O/s) - qflx_irrig => veg_wf%qflx_irrig_patch , & ! Output: [real(r8) (:) ] total water demand or irrigation amount if one-way (mm/s) - qflx_real_irrig => veg_wf%qflx_real_irrig_patch , & ! Output: [real(r8) (:) ] actual irrigation amount (mm/s) - qflx_surf_irrig_col => col_wf%qflx_surf_irrig , & ! Output: [real(r8) (:) ] col real surface water irrigation flux (mm H2O /s) - qflx_grnd_irrig_col => col_wf%qflx_grnd_irrig , & ! Output: [real(r8) (:) ] col real groundwater irrigation flux (mm H2O /s) - qflx_over_supply_col => col_wf%qflx_over_supply , & ! Output: [real(r8) (:) ] over supply irrigation flux (mm H2O /s) - qflx_supply => veg_wf%qflx_supply_patch , & ! Output: [real(r8) (:) ] irrigation supply (mm/s) - qflx_surf_irrig => veg_wf%qflx_surf_irrig_patch , & ! Output: [real(r8) (:) ] actual surface water irrigation (mm/s) - qflx_grnd_irrig => veg_wf%qflx_grnd_irrig_patch , & ! Output: [real(r8) (:) ] actual groundwater irrigation (mm/s) - qflx_over_supply => veg_wf%qflx_over_supply_patch & ! Output: [real(r8) (:) ] the portion of supply that exceed total demand (mm/s) + qflx_irrig => veg_wf%qflx_irrig_patch , & ! Output: [real(r8) (:) ] total water demand or irrigation amount if one-way (mm/s) + qflx_real_irrig => veg_wf%qflx_real_irrig_patch , & ! Output: [real(r8) (:) ] actual irrigation amount (mm/s) + qflx_surf_irrig_col => col_wf%qflx_surf_irrig , & ! Output: [real(r8) (:) ] col real surface water irrigation flux (mm H2O /s) + qflx_grnd_irrig_col => col_wf%qflx_grnd_irrig , & ! Output: [real(r8) (:) ] col real groundwater irrigation flux (mm H2O /s) + qflx_over_supply_col => col_wf%qflx_over_supply , & ! Output: [real(r8) (:) ] over supply irrigation flux (mm H2O /s) + qflx_supply => veg_wf%qflx_supply_patch , & ! Output: [real(r8) (:) ] irrigation supply (mm/s) + qflx_surf_irrig => veg_wf%qflx_surf_irrig_patch , & ! Output: [real(r8) (:) ] actual surface water irrigation (mm/s) + qflx_grnd_irrig => veg_wf%qflx_grnd_irrig_patch , & ! Output: [real(r8) (:) ] actual groundwater irrigation (mm/s) + qflx_over_supply => veg_wf%qflx_over_supply_patch & ! Output: [real(r8) (:) ] the portion of supply that exceed total demand (mm/s) ) ! Compute time step - + dtime = get_step_size() do gg = bounds%begg,bounds%endg irrigated_ppg(gg) = 0 end do - + do f = 1, num_nolakep p = filter_nolakep(f) g = pgridcell(p) @@ -380,26 +380,26 @@ subroutine CanopyHydrology(bounds, & ! Add irrigation water directly onto ground (bypassing canopy interception) ! Note that it's still possible that (some of) this irrigation water will runoff (as runoff is computed later) - if(irrigate) then - if (tw_irr) then ! else one way + if(irrigate) then + if (tw_irr) then ! else one way if (pgwgt(p) > 0._r8) then qflx_supply(p) = atm2lnd_vars%supply_grc(g)*(1-wt_lunit(g,tpu_ind,istdlak))/pgwgt(p) ! original supply at grid level (mm/s) concentrate to pft level. Take lake fraction into consideration else qflx_supply(p) = 0._r8 end if - qflx_real_irrig(p) = 0._r8 + qflx_real_irrig(p) = 0._r8 qflx_surf_irrig(p) = 0._r8 qflx_grnd_irrig(p) = 0._r8 qflx_over_supply(p) = 0._r8 - - if (qflx_irrig(p) > 0._r8 .or. qflx_supply(p) > 0._r8) then !this pft needs water or have supply + + if (qflx_irrig(p) > 0._r8 .or. qflx_supply(p) > 0._r8) then !this pft needs water or have supply if (irrigated(veg_pp%itype(p)) == 1._r8) then ! this pft is irrigated - qflx_surf_irrig(p) = qflx_supply(p)/irrigated_ppg(g) ! surface water irrgation from MOSART, with time step shift + qflx_surf_irrig(p) = qflx_supply(p)/irrigated_ppg(g) ! surface water irrgation from MOSART, with time step shift qflx_grnd_irrig(p) = f_grd(g,tpu_ind)*qflx_irrig(p) ! groundwater irrigation based on demand in ELM, same time step - + if (extra_gw_irr) then ! if always met by additional extra gw pumping - if (qflx_supply(p) > 0._r8) then - if (pgwgt(p) > 0._r8) then + if (qflx_supply(p) > 0._r8) then + if (pgwgt(p) > 0._r8) then qflx_grnd_irrig(p) = atm2lnd_vars%deficit_grc(g)/pgwgt(p)/irrigated_ppg(g) + f_grd(g,tpu_ind)*qflx_irrig(p) !groundwater irrigation based on deficit from MOSART (with time step shift), not demand from ELM end if @@ -407,24 +407,24 @@ subroutine CanopyHydrology(bounds, & qflx_grnd_irrig(p) = f_grd(g,tpu_ind)*qflx_irrig(p) else qflx_grnd_irrig(p) = 0._r8 - endif - endif + endif + endif qflx_real_irrig(p) = qflx_surf_irrig(p) + qflx_grnd_irrig(p) ! actual irrigation, including groundwater irrigation - qflx_prec_grnd_rain(p) = qflx_prec_grnd_rain(p) + qflx_real_irrig(p) - end if - end if + qflx_prec_grnd_rain(p) = qflx_prec_grnd_rain(p) + qflx_real_irrig(p) + end if + end if else ! one way coupling qflx_surf_irrig(p) = f_surf(g,tpu_ind)*qflx_irrig(p) qflx_grnd_irrig(p) = f_grd(g,tpu_ind)*qflx_irrig(p) qflx_real_irrig(p) = qflx_surf_irrig(p) + qflx_grnd_irrig(p) - qflx_prec_grnd_rain(p) = qflx_prec_grnd_rain(p) + qflx_real_irrig(p) + qflx_prec_grnd_rain(p) = qflx_prec_grnd_rain(p) + qflx_real_irrig(p) qflx_over_supply(p) = 0._r8 - qflx_supply(p) = 0._r8 !no water supplied by MOSART + qflx_supply(p) = 0._r8 !no water supplied by MOSART end if else qflx_surf_irrig(p) = 0._r8 qflx_grnd_irrig(p) = 0._r8 - qflx_real_irrig(p) = qflx_surf_irrig(p) + qflx_grnd_irrig(p) + qflx_real_irrig(p) = qflx_surf_irrig(p) + qflx_grnd_irrig(p) qflx_over_supply(p) = 0._r8 qflx_supply(p) = 0._r8 !no water supplied by MOSART end if @@ -448,7 +448,7 @@ subroutine CanopyHydrology(bounds, & end if else qflx_snow_grnd_patch(p) = qflx_prec_grnd_snow(p) ! ice onto ground (mm/s) - qflx_rain_grnd(p) = qflx_prec_grnd_rain(p) ! liquid water onto ground (mm/s) + qflx_rain_grnd(p) = qflx_prec_grnd_rain(p) ! liquid water onto ground (mm/s) endif end do ! (end pft loop) @@ -464,17 +464,17 @@ subroutine CanopyHydrology(bounds, & call p2c(bounds, num_nolakec, filter_nolakec, & qflx_snow_grnd_patch(bounds%begp:bounds%endp), & qflx_snow_grnd_col(bounds%begc:bounds%endc)) - - ! Update column level irrigation supple for balance check - call p2c(bounds, num_nolakec, filter_nolakec, & + + ! Update column level irrigation supple for balance check + call p2c(bounds, num_nolakec, filter_nolakec, & qflx_over_supply(bounds%begp:bounds%endp), & qflx_over_supply_col(bounds%begc:bounds%endc)) - call p2c(bounds, num_nolakec, filter_nolakec, & + call p2c(bounds, num_nolakec, filter_nolakec, & qflx_surf_irrig(bounds%begp:bounds%endp), & qflx_surf_irrig_col(bounds%begc:bounds%endc)) - - call p2c(bounds, num_nolakec, filter_nolakec, & + + call p2c(bounds, num_nolakec, filter_nolakec, & qflx_grnd_irrig(bounds%begp:bounds%endp), & qflx_grnd_irrig_col(bounds%begc:bounds%endc)) @@ -482,7 +482,7 @@ subroutine CanopyHydrology(bounds, & do f = 1, num_nolakec c = filter_nolakec(f) g = cgridcell(c) - if (ctype(c) /= icol_sunwall .and. ctype(c) /= icol_shadewall) then + if (ctype(c) /= icol_sunwall .and. ctype(c) /= icol_shadewall) then qflx_floodc(c) = qflx_floodg(g) else qflx_floodc(c) = 0._r8 @@ -490,7 +490,7 @@ subroutine CanopyHydrology(bounds, & enddo ! Determine snow height and snow water - + if (use_extrasnowlayers) then call NewSnowBulkDensity(bounds, num_nolakec, filter_nolakec, & top_as, bifall(bounds%begc:bounds%endc)) @@ -584,7 +584,7 @@ subroutine CanopyHydrology(bounds, & endif ! use original fsca formulation (n&y 07) - if (oldfflag == 1) then + if (oldfflag == 1) then ! snow cover fraction in Niu et al. 2007 if(snow_depth(c) > 0.0_r8) then frac_sno(c) = tanh(snow_depth(c)/(2.5_r8*zlnd* & @@ -597,7 +597,7 @@ subroutine CanopyHydrology(bounds, & else !h2osno == 0 ! initialize frac_sno and snow_depth when no snow present initially - if (newsnow(c) > 0._r8) then + if (newsnow(c) > 0._r8) then z_avg = newsnow(c)/bifall(c) fmelt=newsnow(c) frac_sno(c) = tanh(accum_factor*newsnow(c)) @@ -615,7 +615,7 @@ subroutine CanopyHydrology(bounds, & snow_depth(c)=newsnow(c)/bifall(c) endif ! use n&y07 formulation - if (oldfflag == 1) then + if (oldfflag == 1) then ! snow cover fraction in Niu et al. 2007 if(snow_depth(c) > 0.0_r8) then frac_sno(c) = tanh(snow_depth(c)/(2.5_r8*zlnd* & @@ -633,7 +633,7 @@ subroutine CanopyHydrology(bounds, & qflx_snow_h2osfc(c) = 0._r8 ! update h2osno for new snow - h2osno(c) = h2osno(c) + newsnow(c) + h2osno(c) = h2osno(c) + newsnow(c) int_snow(c) = int_snow(c) + newsnow(c) ! update change in snow depth @@ -643,7 +643,7 @@ subroutine CanopyHydrology(bounds, & ! set frac_sno_eff variable if (ltype(l) == istsoil .or. ltype(l) == istcrop) then - if (subgridflag ==1) then + if (subgridflag ==1) then frac_sno_eff(c) = frac_sno(c) else frac_sno_eff(c) = 1._r8 @@ -695,7 +695,7 @@ subroutine CanopyHydrology(bounds, & call aerosol_vars%Reset(column=c) ! call waterstate_vars%Reset(column=c) col_ws%snw_rds(c,0) = snw_rds_min - end if + end if end if ! The change of ice partial density of surface node due to precipitation. @@ -713,7 +713,7 @@ subroutine CanopyHydrology(bounds, & call FracH2oSfc(bounds, num_nolakec, filter_nolakec, & col_wf%qflx_h2osfc2topsoi, dtime) - end associate + end associate end subroutine CanopyHydrology @@ -739,15 +739,15 @@ subroutine FracWet(numf, filter, canopystate_vars ) real(r8) :: dewmxi ! inverse of maximum allowed dew [1/mm] !----------------------------------------------------------------------- - associate( & + associate( & frac_veg_nosno => canopystate_vars%frac_veg_nosno_patch , & ! Input: [integer (:)] fraction of veg not covered by snow (0/1 now) [-] - dewmx => canopystate_vars%dewmx_patch , & ! Input: [real(r8) (:) ] Maximum allowed dew [mm] + dewmx => canopystate_vars%dewmx_patch , & ! Input: [real(r8) (:) ] Maximum allowed dew [mm] elai => canopystate_vars%elai_patch , & ! Input: [real(r8) (:) ] one-sided leaf area index with burying by snow esai => canopystate_vars%esai_patch , & ! Input: [real(r8) (:) ] one-sided stem area index with burying by snow - h2ocan => veg_ws%h2ocan , & ! Input: [real(r8) (:) ] total canopy water (mm H2O) - - fwet => veg_ws%fwet , & ! Output: [real(r8) (:) ] fraction of canopy that is wet (0 to 1) + h2ocan => veg_ws%h2ocan , & ! Input: [real(r8) (:) ] total canopy water (mm H2O) + + fwet => veg_ws%fwet , & ! Output: [real(r8) (:) ] fraction of canopy that is wet (0 to 1) fdry => veg_ws%fdry & ! Output: [real(r8) (:) ] fraction of foliage that is green and dry [-] (new) ) @@ -777,7 +777,7 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, & qflx_h2osfc2topsoi,dtime, no_update) ! ! !DESCRIPTION: - ! Determine fraction of land surfaces which are submerged + ! Determine fraction of land surfaces which are submerged ! based on surface microtopography and surface water storage. ! ! !USES: @@ -786,11 +786,11 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, & use landunit_varcon , only : istsoil, istcrop ! ! !ARGUMENTS: - type(bounds_type) , intent(in) :: bounds + type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_h2osfc ! number of column points in column filter - integer , intent(in) :: filter_h2osfc(:) ! column filter - real(r8) , intent(inout) :: qflx_h2osfc2topsoi(bounds%begc:bounds%endc) - real(r8) , intent(in) :: dtime + integer , intent(in) :: filter_h2osfc(:) ! column filter + real(r8) , intent(inout) :: qflx_h2osfc2topsoi(bounds%begc:bounds%endc) + real(r8) , intent(in) :: dtime integer , intent(in), optional :: no_update ! flag to make calculation w/o updating variables ! ! !LOCAL VARIABLES: @@ -800,17 +800,18 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, & real(r8):: min_h2osfc !----------------------------------------------------------------------- - associate( & - micro_sigma => col_pp%micro_sigma , & ! Input: [real(r8) (:) ] microtopography pdf sigma (m) - - h2osno => col_ws%h2osno , & ! Input: [real(r8) (:) ] snow water (mm H2O) - - h2osoi_liq => col_ws%h2osoi_liq , & ! Output: [real(r8) (:,:) ] liquid water (col,lyr) [kg/m2] - h2osfc => col_ws%h2osfc , & ! Output: [real(r8) (:) ] surface water (mm) - frac_sno => col_ws%frac_sno , & ! Output: [real(r8) (:) ] fraction of ground covered by snow (0 to 1) - frac_sno_eff => col_ws%frac_sno_eff , & ! Output: [real(r8) (:) ] eff. fraction of ground covered by snow (0 to 1) - frac_h2osfc => col_ws%frac_h2osfc , & ! Output: [real(r8) (:) ] col fractional area with surface water greater than zero - frac_h2osfc_act => col_ws%frac_h2osfc_act & ! Output: [real(r8) (:) ] col fractional area with surface water greater than zero + associate( & + micro_sigma => col_pp%micro_sigma , & ! Input: [real(r8) (:) ] microtopography pdf sigma (m) + + h2osno => col_ws%h2osno , & ! Input: [real(r8) (:) ] snow water (mm H2O) + excess_ice => col_ws%excess_ice , & ! Input: [real(r8) (:) ] excess ground ice [kg/m2] + + h2osoi_liq => col_ws%h2osoi_liq , & ! Output: [real(r8) (:,:) ] liquid water (col,lyr) [kg/m2] + h2osfc => col_ws%h2osfc , & ! Output: [real(r8) (:) ] surface water (mm) + frac_sno => col_ws%frac_sno , & ! Output: [real(r8) (:) ] fraction of ground covered by snow (0 to 1) + frac_sno_eff => col_ws%frac_sno_eff , & ! Output: [real(r8) (:) ] eff. fraction of ground covered by snow (0 to 1) + frac_h2osfc => col_ws%frac_h2osfc , & ! Output: [real(r8) (:) ] col fractional area with surface water greater than zero + frac_h2osfc_act => col_ws%frac_h2osfc_act & ! Output: [real(r8) (:) ] col fractional area with surface water greater than zero ) ! arbitrary lower limit on h2osfc for safer numerics... @@ -824,7 +825,7 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, & if (lun_pp%itype(l) == istsoil .or. lun_pp%itype(l) == istcrop) then ! Use newton-raphson method to iteratively determine frac_h20sfc - ! based on amount of surface water storage (h2osfc) and + ! based on amount of surface water storage (h2osfc) and ! microtopography variability (micro_sigma) if (h2osfc(c) > min_h2osfc) then @@ -846,10 +847,10 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, & else frac_h2osfc(c) = 0._r8 h2osoi_liq(c,1) = h2osoi_liq(c,1) + h2osfc(c) - qflx_h2osfc2topsoi(c) = h2osfc(c)/dtime + qflx_h2osfc2topsoi(c) = h2osfc(c)/dtime h2osfc(c)=0._r8 endif - + frac_h2osfc_act(c) = frac_h2osfc(c) if (.not. present(no_update)) then @@ -857,7 +858,7 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, & ! adjust fh2o, fsno when sum is greater than zero if (frac_sno(c) > (1._r8 - frac_h2osfc(c)) .and. h2osno(c) > 0) then - if (frac_h2osfc(c) > 0.01_r8) then + if (frac_h2osfc(c) > 0.01_r8) then frac_h2osfc(c) = max(1.0_r8 - frac_sno(c),0.01_r8) frac_sno(c) = 1.0_r8 - frac_h2osfc(c) else @@ -873,7 +874,7 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, & frac_h2osfc(c) = 0._r8 frac_h2osfc_act(c) = 0._r8 - + endif end do diff --git a/components/elm/src/biogeophys/WaterStateType.F90 b/components/elm/src/biogeophys/WaterStateType.F90 index c8241db7ab11..1233b24d863e 100644 --- a/components/elm/src/biogeophys/WaterStateType.F90 +++ b/components/elm/src/biogeophys/WaterStateType.F90 @@ -12,10 +12,10 @@ module WaterstateType use decompMod , only : bounds_type use elm_varctl , only : use_vancouver, use_mexicocity, use_cn, iulog, use_fates_planthydro, & use_hydrstress, use_fan - use elm_varpar , only : nlevgrnd, nlevurb, nlevsno + use elm_varpar , only : nlevgrnd, nlevurb, nlevsno use elm_varcon , only : spval - use LandunitType , only : lun_pp - use ColumnType , only : col_pp + use LandunitType , only : lun_pp + use ColumnType , only : col_pp ! implicit none save @@ -32,24 +32,25 @@ module WaterstateType real(r8), pointer :: snowliq_col (:) ! col average snow liquid water real(r8), pointer :: int_snow_col (:) ! col integrated snowfall (mm H2O) real(r8), pointer :: snow_layer_unity_col (:,:) ! value 1 for each snow layer, used for history diagnostics - real(r8), pointer :: bw_col (:,:) ! col partial density of water in the snow pack (ice + liquid) [kg/m3] + real(r8), pointer :: bw_col (:,:) ! col partial density of water in the snow pack (ice + liquid) [kg/m3] real(r8), pointer :: finundated_col (:) ! fraction of column that is inundated, this is for bgc caclulation in betr real(r8), pointer :: h2osoi_tend_tsl_col (:) ! col moisture tendency due to vertical movement at topmost layer (m3/m3/s) - + real(r8), pointer :: excess_ice (:) ! excess ground ice in polygonal tundra areas [kg/m2] + real(r8), pointer :: rhvap_soi_col (:,:) real(r8), pointer :: rho_vap_col (:,:) real(r8), pointer :: smp_l_col (:,:) ! col liquid phase soil matric potential, mm real(r8), pointer :: h2osno_col (:) ! col snow water (mm H2O) real(r8), pointer :: h2osno_old_col (:) ! col snow mass for previous time step (kg/m2) (new) - real(r8), pointer :: h2osoi_liq_col (:,:) ! col liquid water (kg/m2) (new) (-nlevsno+1:nlevgrnd) - real(r8), pointer :: h2osoi_ice_col (:,:) ! col ice lens (kg/m2) (new) (-nlevsno+1:nlevgrnd) + real(r8), pointer :: h2osoi_liq_col (:,:) ! col liquid water (kg/m2) (new) (-nlevsno+1:nlevgrnd) + real(r8), pointer :: h2osoi_ice_col (:,:) ! col ice lens (kg/m2) (new) (-nlevsno+1:nlevgrnd) real(r8), pointer :: h2osoi_liq_old_col (:,:) real(r8), pointer :: h2osoi_ice_old_col (:,:) real(r8), pointer :: h2osoi_liqice_10cm_col (:) ! col liquid water + ice lens in top 10cm of soil (kg/m2) real(r8), pointer :: h2osoi_vol_col (:,:) ! col volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3] (nlevgrnd) - real(r8), pointer :: air_vol_col (:,:) ! col air filled porosity + real(r8), pointer :: air_vol_col (:,:) ! col air filled porosity real(r8), pointer :: h2osoi_liqvol_col (:,:) ! col volumetric liquid water content (v/v) - real(r8), pointer :: h2osoi_icevol_col (:,:) ! col volumetric ice content (v/v) + real(r8), pointer :: h2osoi_icevol_col (:,:) ! col volumetric ice content (v/v) real(r8), pointer :: h2ocan_patch (:) ! patch canopy water (mm H2O) real(r8), pointer :: h2ocan_col (:) ! col canopy water (mm H2O) real(r8), pointer :: h2osfc_col (:) ! col surface water (mm H2O) @@ -88,10 +89,10 @@ module WaterstateType ! Fractions real(r8), pointer :: frac_sno_col (:) ! col fraction of ground covered by snow (0 to 1) real(r8), pointer :: frac_sno_eff_col (:) ! col fraction of ground covered by snow (0 to 1) - real(r8), pointer :: frac_iceold_col (:,:) ! col fraction of ice relative to the tot water (new) (-nlevsno+1:nlevgrnd) + real(r8), pointer :: frac_iceold_col (:,:) ! col fraction of ice relative to the tot water (new) (-nlevsno+1:nlevgrnd) real(r8), pointer :: frac_h2osfc_col (:) ! col fractional area with surface water greater than zero - real(r8), pointer :: wf_col (:) ! col soil water as frac. of whc for top 0.05 m (0-1) - real(r8), pointer :: wf2_col (:) ! col soil water as frac. of whc for top 0.17 m (0-1) + real(r8), pointer :: wf_col (:) ! col soil water as frac. of whc for top 0.05 m (0-1) + real(r8), pointer :: wf2_col (:) ! col soil water as frac. of whc for top 0.17 m (0-1) real(r8), pointer :: fwet_patch (:) ! patch canopy fraction that is wet (0 to 1) real(r8), pointer :: fdry_patch (:) ! patch canopy fraction of foliage that is green and dry [-] (new) @@ -133,17 +134,17 @@ module WaterstateType contains - procedure :: Init - procedure :: Restart - procedure, public :: Reset - procedure, private :: InitAllocate - procedure, private :: InitHistory - procedure, private :: InitCold + procedure :: Init + procedure :: Restart + procedure, public :: Reset + procedure, private :: InitAllocate + procedure, private :: InitHistory + procedure, private :: InitCold procedure, public :: save_h2osoi_old end type waterstate_type ! minimum allowed snow effective radius (also "fresh snow" value) [microns] - real(r8), public, parameter :: snw_rds_min = 54.526_r8 - + real(r8), public, parameter :: snw_rds_min = 54.526_r8 + type(waterstate_type) ,public :: waterstate_vars !------------------------------------------------------------------------ @@ -154,7 +155,7 @@ subroutine Init(this, bounds, & h2osno_input_col, snow_depth_input_col, watsat_col, t_soisno_col) class(waterstate_type) :: this - type(bounds_type) , intent(in) :: bounds + type(bounds_type) , intent(in) :: bounds real(r8) , intent(inout) :: h2osno_input_col(bounds%begc:) real(r8) , intent(inout) :: snow_depth_input_col(bounds%begc:) !real(r8) , intent(inout) :: watsat_col(bounds%begc:, 1:) ! volumetric soil water at saturation (porosity) @@ -162,7 +163,7 @@ subroutine Init(this, bounds, & real(r8) , intent(inout) :: watsat_col(bounds%begc:bounds%endc, 1:nlevgrnd) ! volumetric soil water at saturation (porosity) real(r8) , intent(inout) :: t_soisno_col(bounds%begc:bounds%endc, -nlevsno+1:nlevgrnd) ! col soil temperature (Kelvin) - call this%InitAllocate(bounds) + call this%InitAllocate(bounds) call this%InitHistory(bounds) @@ -182,7 +183,7 @@ subroutine InitAllocate(this, bounds) ! ! !ARGUMENTS: class(waterstate_type) :: this - type(bounds_type), intent(in) :: bounds + type(bounds_type), intent(in) :: bounds ! ! !LOCAL VARIABLES: integer :: begp, endp @@ -197,39 +198,40 @@ subroutine InitAllocate(this, bounds) begl = bounds%begl; endl= bounds%endl begg = bounds%begg; endg= bounds%endg - allocate(this%do_capsnow_col (begc:endc)) + allocate(this%do_capsnow_col (begc:endc)) allocate(this%snow_depth_col (begc:endc)) ; this%snow_depth_col (:) = nan allocate(this%snow_persistence_col (begc:endc)) ; this%snow_persistence_col (:) = nan allocate(this%snowdp_col (begc:endc)) ; this%snowdp_col (:) = nan - allocate(this%snowice_col (begc:endc)) ; this%snowice_col (:) = nan - allocate(this%snowliq_col (begc:endc)) ; this%snowliq_col (:) = nan - allocate(this%int_snow_col (begc:endc)) ; this%int_snow_col (:) = nan + allocate(this%snowice_col (begc:endc)) ; this%snowice_col (:) = nan + allocate(this%snowliq_col (begc:endc)) ; this%snowliq_col (:) = nan + allocate(this%int_snow_col (begc:endc)) ; this%int_snow_col (:) = nan allocate(this%snow_layer_unity_col (begc:endc,-nlevsno+1:0)) ; this%snow_layer_unity_col (:,:) = nan - allocate(this%bw_col (begc:endc,-nlevsno+1:0)) ; this%bw_col (:,:) = nan + allocate(this%bw_col (begc:endc,-nlevsno+1:0)) ; this%bw_col (:,:) = nan allocate(this%smp_l_col (begc:endc,-nlevsno+1:nlevgrnd)) ; this%smp_l_col (:,:) = nan allocate(this%finundated_col (begc:endc)) ; this%finundated_col (:) = nan if (use_fan) then allocate(this%h2osoi_tend_tsl_col(begc:endc)) ; this%h2osoi_tend_tsl_col (:) = nan end if + allocate(this%excess_ice (begc:endc)) ; this%excess_ice (:) = nan - allocate(this%h2osno_col (begc:endc)) ; this%h2osno_col (:) = nan - allocate(this%h2osno_old_col (begc:endc)) ; this%h2osno_old_col (:) = nan + allocate(this%h2osno_col (begc:endc)) ; this%h2osno_col (:) = nan + allocate(this%h2osno_old_col (begc:endc)) ; this%h2osno_old_col (:) = nan allocate(this%h2osoi_liqice_10cm_col (begc:endc)) ; this%h2osoi_liqice_10cm_col (:) = nan - + allocate(this%h2osoi_vol_col (begc:endc, 1:nlevgrnd)) ; this%h2osoi_vol_col (:,:) = nan allocate(this%air_vol_col (begc:endc, 1:nlevgrnd)) ; this%air_vol_col (:,:) = nan allocate(this%h2osoi_liqvol_col (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_liqvol_col (:,:) = nan - allocate(this%h2osoi_icevol_col (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_icevol_col (:,:) = nan + allocate(this%h2osoi_icevol_col (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_icevol_col (:,:) = nan allocate(this%rho_vap_col (begc:endc,-nlevsno+1:nlevgrnd)) ; this%rho_vap_col (:,:) = nan allocate(this%rhvap_soi_col (begc:endc,-nlevsno+1:nlevgrnd)) ; this%rhvap_soi_col (:,:) = nan allocate(this%h2osoi_liq_old_col (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_liq_old_col (:,:) = nan - allocate(this%h2osoi_ice_old_col (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_ice_old_col (:,:) = nan + allocate(this%h2osoi_ice_old_col (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_ice_old_col (:,:) = nan allocate(this%h2osoi_ice_col (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_ice_col (:,:) = nan allocate(this%h2osoi_liq_col (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_liq_col (:,:) = nan - allocate(this%h2ocan_patch (begp:endp)) ; this%h2ocan_patch (:) = nan - allocate(this%h2ocan_col (begc:endc)) ; this%h2ocan_col (:) = nan - allocate(this%h2osfc_col (begc:endc)) ; this%h2osfc_col (:) = nan - allocate(this%swe_old_col (begc:endc,-nlevsno+1:0)) ; this%swe_old_col (:,:) = nan + allocate(this%h2ocan_patch (begp:endp)) ; this%h2ocan_patch (:) = nan + allocate(this%h2ocan_col (begc:endc)) ; this%h2ocan_col (:) = nan + allocate(this%h2osfc_col (begc:endc)) ; this%h2osfc_col (:) = nan + allocate(this%swe_old_col (begc:endc,-nlevsno+1:0)) ; this%swe_old_col (:,:) = nan allocate(this%liq1_grc (begg:endg)) ; this%liq1_grc (:) = nan allocate(this%liq2_grc (begg:endg)) ; this%liq2_grc (:) = nan allocate(this%ice1_grc (begg:endg)) ; this%ice1_grc (:) = nan @@ -245,11 +247,11 @@ subroutine InitAllocate(this, bounds) allocate(this%h2osno_top_col (begc:endc)) ; this%h2osno_top_col (:) = nan allocate(this%sno_liq_top_col (begc:endc)) ; this%sno_liq_top_col (:) = nan - allocate(this%qg_snow_col (begc:endc)) ; this%qg_snow_col (:) = nan - allocate(this%qg_soil_col (begc:endc)) ; this%qg_soil_col (:) = nan - allocate(this%qg_h2osfc_col (begc:endc)) ; this%qg_h2osfc_col (:) = nan - allocate(this%qg_col (begc:endc)) ; this%qg_col (:) = nan - allocate(this%dqgdT_col (begc:endc)) ; this%dqgdT_col (:) = nan + allocate(this%qg_snow_col (begc:endc)) ; this%qg_snow_col (:) = nan + allocate(this%qg_soil_col (begc:endc)) ; this%qg_soil_col (:) = nan + allocate(this%qg_h2osfc_col (begc:endc)) ; this%qg_h2osfc_col (:) = nan + allocate(this%qg_col (begc:endc)) ; this%qg_col (:) = nan + allocate(this%dqgdT_col (begc:endc)) ; this%dqgdT_col (:) = nan allocate(this%qaf_lun (begl:endl)) ; this%qaf_lun (:) = nan allocate(this%q_ref2m_patch (begp:endp)) ; this%q_ref2m_patch (:) = nan allocate(this%rh_ref2m_patch (begp:endp)) ; this%rh_ref2m_patch (:) = nan @@ -260,9 +262,9 @@ subroutine InitAllocate(this, bounds) allocate(this%frac_sno_col (begc:endc)) ; this%frac_sno_col (:) = nan allocate(this%frac_sno_eff_col (begc:endc)) ; this%frac_sno_eff_col (:) = nan allocate(this%frac_iceold_col (begc:endc,-nlevsno+1:nlevgrnd)) ; this%frac_iceold_col (:,:) = nan - allocate(this%frac_h2osfc_col (begc:endc)) ; this%frac_h2osfc_col (:) = nan + allocate(this%frac_h2osfc_col (begc:endc)) ; this%frac_h2osfc_col (:) = nan allocate(this%wf_col (begc:endc)) ; this%wf_col (:) = nan - allocate(this%wf2_col (begc:endc)) ; + allocate(this%wf2_col (begc:endc)) ; allocate(this%fwet_patch (begp:endp)) ; this%fwet_patch (:) = nan allocate(this%fdry_patch (begp:endp)) ; this%fdry_patch (:) = nan @@ -312,12 +314,12 @@ subroutine InitHistory(this, bounds) use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) use elm_varctl , only : create_glacier_mec_landunit, use_cn, use_lch4 use elm_varctl , only : hist_wrtch4diag, use_lake_wat_storage - use elm_varpar , only : nlevsno, crop_prog + use elm_varpar , only : nlevsno, crop_prog use histFileMod , only : hist_addfld1d, hist_addfld2d, no_snow_normal, no_snow_zero ! ! !ARGUMENTS: class(waterstate_type) :: this - type(bounds_type), intent(in) :: bounds + type(bounds_type), intent(in) :: bounds ! ! !LOCAL VARIABLES: integer :: begp, endp @@ -332,8 +334,8 @@ subroutine InitHistory(this, bounds) begg = bounds%begg; endg= bounds%endg - ! h2osno also includes snow that is part of the soil column (an - ! initial snow layer is only created if h2osno > 10mm). + ! h2osno also includes snow that is part of the soil column (an + ! initial snow layer is only created if h2osno > 10mm). ! Snow properties - these will be vertically averaged over the snow profile @@ -353,7 +355,7 @@ subroutine InitCold(this, bounds, & h2osno_input_col, snow_depth_input_col, watsat_col, t_soisno_col) ! ! !DESCRIPTION: - ! Initialize time constant variables and cold start conditions + ! Initialize time constant variables and cold start conditions ! ! !USES: use shr_const_mod , only : shr_const_pi @@ -362,10 +364,10 @@ subroutine InitCold(this, bounds, & use shr_kind_mod , only : r8 => shr_kind_r8 use shr_const_mod , only : SHR_CONST_TKFRZ use elm_varpar , only : nlevsoi, nlevgrnd, nlevsno, nlevlak, nlevurb - use landunit_varcon , only : istice, istwet, istsoil, istdlak, istcrop, istice_mec + use landunit_varcon , only : istice, istwet, istsoil, istdlak, istcrop, istice_mec use column_varcon , only : icol_shadewall, icol_road_perv use column_varcon , only : icol_road_imperv, icol_roof, icol_sunwall - use elm_varcon , only : denice, denh2o, spval, sb, bdsno + use elm_varcon , only : denice, denh2o, spval, sb, bdsno use elm_varcon , only : h2osno_max, zlnd, tfrz, spval use elm_varctl , only : fsurdat, iulog use spmdMod , only : masterproc @@ -382,13 +384,13 @@ subroutine InitCold(this, bounds, & real(r8) , intent(in) :: t_soisno_col(bounds%begc:, -nlevsno+1:) ! col soil temperature (Kelvin) ! ! !LOCAL VARIABLES: - integer :: p,c,j,l,g,lev,nlevs,nlevbed + integer :: p,c,j,l,g,lev,nlevs,nlevbed real(r8) :: maxslope, slopemax, minslope real(r8) :: d, fd, dfdd, slope0,slopebeta - real(r8) ,pointer :: std (:) - logical :: readvar - type(file_desc_t) :: ncid - character(len=256) :: locfn + real(r8) ,pointer :: std (:) + logical :: readvar + type(file_desc_t) :: ncid + character(len=256) :: locfn real(r8) :: snowbd ! temporary calculation of snow bulk density (kg/m3) real(r8) :: fmelt ! snowbd/100 !----------------------------------------------------------------------- @@ -400,8 +402,8 @@ subroutine InitCold(this, bounds, & ! The first three arrays are initialized from the input argument do c = bounds%begc,bounds%endc - this%h2osno_col(c) = h2osno_input_col(c) - this%int_snow_col(c) = h2osno_input_col(c) + this%h2osno_col(c) = h2osno_input_col(c) + this%int_snow_col(c) = h2osno_input_col(c) this%snow_depth_col(c) = snow_depth_input_col(c) this%snow_persistence_col(c) = 0._r8 this%snow_layer_unity_col(c,:) = 1._r8 @@ -412,7 +414,7 @@ subroutine InitCold(this, bounds, & this%wf2_col(c) = spval end do - do l = bounds%begl, bounds%endl + do l = bounds%begl, bounds%endl if (lun_pp%urbpoi(l)) then if (use_vancouver) then this%qaf_lun(l) = 0.0111_r8 @@ -427,8 +429,8 @@ subroutine InitCold(this, bounds, & ! Water Stored in plants is almost always a static entity, with the exception ! of when FATES-hydraulics is used. As such, this is trivially set to 0.0 (rgk 03-2017) this%total_plant_stored_h2o_col(bounds%begc:bounds%endc) = 0.0_r8 - - associate(snl => col_pp%snl, nlev2bed => col_pp%nlevbed) + + associate(snl => col_pp%snl, nlev2bed => col_pp%nlevbed) this%h2osfc_col(bounds%begc:bounds%endc) = 0._r8 this%h2ocan_patch(bounds%begp:bounds%endp) = 0._r8 @@ -491,14 +493,14 @@ end subroutine InitCold !------------------------------------------------------------------------ subroutine Restart(this, bounds, ncid, flag, & watsat_col) - ! + ! ! !DESCRIPTION: ! Read/Write module information to/from restart file. ! ! !USES: use spmdMod , only : masterproc - use elm_varcon , only : denice, denh2o, pondmx, watmin, spval - use landunit_varcon , only : istcrop, istdlak, istsoil + use elm_varcon , only : denice, denh2o, pondmx, watmin, spval + use landunit_varcon , only : istcrop, istdlak, istsoil use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall use elm_time_manager , only : is_first_step use elm_varctl , only : bound_h2osoi, use_lake_wat_storage @@ -508,7 +510,7 @@ subroutine Restart(this, bounds, ncid, flag, & ! ! !ARGUMENTS: class(waterstate_type) :: this - type(bounds_type), intent(in) :: bounds + type(bounds_type), intent(in) :: bounds type(file_desc_t), intent(inout) :: ncid ! netcdf id character(len=*) , intent(in) :: flag ! 'read' or 'write' real(r8) , intent(in) :: watsat_col (bounds%begc:, 1:) ! volumetric soil water at saturation (porosity) @@ -516,7 +518,7 @@ subroutine Restart(this, bounds, ncid, flag, & ! !LOCAL VARIABLES: integer :: c,l,j,nlevs logical :: readvar - real(r8) :: maxwatsat ! maximum porosity + real(r8) :: maxwatsat ! maximum porosity real(r8) :: excess ! excess volumetric soil water real(r8) :: totwat ! total soil water (mm) !------------------------------------------------------------------------ @@ -552,7 +554,7 @@ subroutine save_h2osoi_old(this,bounds) ! ! !ARGUMENTS: class(waterstate_type) :: this - type(bounds_type) , intent(in) :: bounds + type(bounds_type) , intent(in) :: bounds integer :: begc, endc @@ -561,6 +563,6 @@ subroutine save_h2osoi_old(this,bounds) this%h2osoi_liq_old_col(begc:endc,:) = this%h2osoi_liq_col(begc:endc,:) this%h2osoi_ice_old_col(begc:endc,:) = this%h2osoi_ice_col(begc:endc,:) - + end subroutine save_h2osoi_old end module WaterstateType diff --git a/components/elm/src/data_types/ColumnDataType.F90 b/components/elm/src/data_types/ColumnDataType.F90 index 66f167fbb4c0..0d968c0cb899 100644 --- a/components/elm/src/data_types/ColumnDataType.F90 +++ b/components/elm/src/data_types/ColumnDataType.F90 @@ -1568,7 +1568,7 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ this%excess_ice(begc:endc) = spval call hist_addfld1d (fname='EXCESS_ICE', units = 'kg/m2', & avgflag='A', long_name='Excess ground ice', & - ptr_col=this$excess_ice, l2g_scale_type='veg') ! <- RPF: should this be natveg? + ptr_col=this%excess_ice, l2g_scale_type='veg') ! <- RPF: should this be natveg? this%frac_sno(begc:endc) = spval call hist_addfld1d (fname='FSNO', units='1', & @@ -1641,6 +1641,7 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ this%frac_h2osfc_act(c) = 0._r8 this%h2orof(c) = 0._r8 this%frac_h2orof(c) = 0._r8 + this%excess_ice(c) = 500._r8 if (lun_pp%urbpoi(l)) then ! From Bonan 1996 (LSM technical note) From 7a7ae66f33a090feaad088e3f2f57954af57f29e Mon Sep 17 00:00:00 2001 From: Rich Fiorella Date: Tue, 16 Apr 2024 06:56:30 -0600 Subject: [PATCH 6/7] Change excess ground ice units to unitless Excess ground ice had been templated in these changes as mass/area - changing to a volume/volume basis (and therefore, now unitless). --- components/elm/src/data_types/ColumnDataType.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/components/elm/src/data_types/ColumnDataType.F90 b/components/elm/src/data_types/ColumnDataType.F90 index 0d968c0cb899..d83949b2016b 100644 --- a/components/elm/src/data_types/ColumnDataType.F90 +++ b/components/elm/src/data_types/ColumnDataType.F90 @@ -1566,8 +1566,8 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ end if this%excess_ice(begc:endc) = spval - call hist_addfld1d (fname='EXCESS_ICE', units = 'kg/m2', & - avgflag='A', long_name='Excess ground ice', & + call hist_addfld1d (fname='EXCESS_ICE', units = '1', & + avgflag='A', long_name='Excess ground ice (0 to 1)', & ptr_col=this%excess_ice, l2g_scale_type='veg') ! <- RPF: should this be natveg? this%frac_sno(begc:endc) = spval @@ -1641,7 +1641,7 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ this%frac_h2osfc_act(c) = 0._r8 this%h2orof(c) = 0._r8 this%frac_h2orof(c) = 0._r8 - this%excess_ice(c) = 500._r8 + this%excess_ice(c) = 0.5_r8 if (lun_pp%urbpoi(l)) then ! From Bonan 1996 (LSM technical note) @@ -1909,7 +1909,7 @@ subroutine col_ws_restart(this, bounds, ncid, flag, watsat_input) call restartvar(ncid=ncid, flag=flag, varname='EXCESS_ICE', xtype=ncd_double, & dim1name='column', & - long_name='excess ground ice', units='kg/m2', & + long_name='excess ground ice (0 to 1)', units='1', & interpinic_flag='interp', readvar=readvar, data=this%excess_ice) call restartvar(ncid=ncid, flag=flag, varname='frac_sno', xtype=ncd_double, & From f6c4c8916fce3a0a772793e9b1e3ff08cd6d6ed4 Mon Sep 17 00:00:00 2001 From: Rich Fiorella Date: Wed, 1 May 2024 10:29:45 -0600 Subject: [PATCH 7/7] Update excess_ice from single column value to depth varying Changes dimension on excess_ice parameter in water state types from being one value per column to having dimensions (columns, depth in soil) --- .../elm/src/biogeophys/CanopyHydrologyMod.F90 | 2 +- .../elm/src/biogeophys/WaterStateType.F90 | 4 +- .../elm/src/data_types/ColumnDataType.F90 | 42 ++++++++++++------- 3 files changed, 29 insertions(+), 19 deletions(-) diff --git a/components/elm/src/biogeophys/CanopyHydrologyMod.F90 b/components/elm/src/biogeophys/CanopyHydrologyMod.F90 index 13c83ed010e5..1f391bd67b1f 100755 --- a/components/elm/src/biogeophys/CanopyHydrologyMod.F90 +++ b/components/elm/src/biogeophys/CanopyHydrologyMod.F90 @@ -804,7 +804,7 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, & micro_sigma => col_pp%micro_sigma , & ! Input: [real(r8) (:) ] microtopography pdf sigma (m) h2osno => col_ws%h2osno , & ! Input: [real(r8) (:) ] snow water (mm H2O) - excess_ice => col_ws%excess_ice , & ! Input: [real(r8) (:) ] excess ground ice [kg/m2] + excess_ice => col_ws%excess_ice , & ! Input: [real(r8) (:,:) ] excess ground ice [kg/m2] h2osoi_liq => col_ws%h2osoi_liq , & ! Output: [real(r8) (:,:) ] liquid water (col,lyr) [kg/m2] h2osfc => col_ws%h2osfc , & ! Output: [real(r8) (:) ] surface water (mm) diff --git a/components/elm/src/biogeophys/WaterStateType.F90 b/components/elm/src/biogeophys/WaterStateType.F90 index 1233b24d863e..e223a0ec53ee 100644 --- a/components/elm/src/biogeophys/WaterStateType.F90 +++ b/components/elm/src/biogeophys/WaterStateType.F90 @@ -35,7 +35,7 @@ module WaterstateType real(r8), pointer :: bw_col (:,:) ! col partial density of water in the snow pack (ice + liquid) [kg/m3] real(r8), pointer :: finundated_col (:) ! fraction of column that is inundated, this is for bgc caclulation in betr real(r8), pointer :: h2osoi_tend_tsl_col (:) ! col moisture tendency due to vertical movement at topmost layer (m3/m3/s) - real(r8), pointer :: excess_ice (:) ! excess ground ice in polygonal tundra areas [kg/m2] + real(r8), pointer :: excess_ice (:,:) ! excess ground ice in polygonal tundra areas [kg/m2] real(r8), pointer :: rhvap_soi_col (:,:) real(r8), pointer :: rho_vap_col (:,:) @@ -212,7 +212,7 @@ subroutine InitAllocate(this, bounds) if (use_fan) then allocate(this%h2osoi_tend_tsl_col(begc:endc)) ; this%h2osoi_tend_tsl_col (:) = nan end if - allocate(this%excess_ice (begc:endc)) ; this%excess_ice (:) = nan + allocate(this%excess_ice (begc:endc, 1:nlevgrnd)) ; this%excess_ice (:,:) = nan allocate(this%h2osno_col (begc:endc)) ; this%h2osno_col (:) = nan allocate(this%h2osno_old_col (begc:endc)) ; this%h2osno_old_col (:) = nan diff --git a/components/elm/src/data_types/ColumnDataType.F90 b/components/elm/src/data_types/ColumnDataType.F90 index d83949b2016b..2523230868b4 100644 --- a/components/elm/src/data_types/ColumnDataType.F90 +++ b/components/elm/src/data_types/ColumnDataType.F90 @@ -108,6 +108,7 @@ module ColumnDataType real(r8), pointer :: h2osoi_liq (:,:) => null() ! liquid water (-nlevsno+1:nlevgrnd) (kg/m2) real(r8), pointer :: h2osoi_ice (:,:) => null() ! ice lens (-nlevsno+1:nlevgrnd) (kg/m2) real(r8), pointer :: h2osoi_vol (:,:) => null() ! volumetric soil water (0<=h2osoi_vol<=watsat) (1:nlevgrnd) (m3/m3) + real(r8), pointer :: excess_ice (:,:) => null() ! NGEE Arctic: excess ground ice in column (1:nlevgrnd) (0 to 1) real(r8), pointer :: h2osfc (:) => null() ! surface water (kg/m2) real(r8), pointer :: h2ocan (:) => null() ! canopy water integrated to column (kg/m2) real(r8), pointer :: total_plant_stored_h2o(:)=> null() ! total water in plants (kg/m2) @@ -143,7 +144,6 @@ module ColumnDataType real(r8), pointer :: snw_rds_top (:) => null() ! snow grain radius (top layer) (m^-6, microns) logical , pointer :: do_capsnow (:) => null() ! true => do snow capping real(r8), pointer :: h2osoi_tend_tsl_col(:) => null() ! col moisture tendency due to vertical movement at topmost layer (m3/m3/s) - real(r8), pointer :: excess_ice (:) => null() ! NGEE-Arctic: tracking excess ground ice ! Area fractions real(r8), pointer :: frac_sno (:) => null() ! fraction of ground covered by snow (0 to 1) real(r8), pointer :: frac_sno_eff (:) => null() ! fraction of ground covered by snow (0 to 1) @@ -1384,6 +1384,7 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ allocate(this%h2osoi_liq (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_liq (:,:) = spval allocate(this%h2osoi_ice (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_ice (:,:) = spval allocate(this%h2osoi_vol (begc:endc, 1:nlevgrnd)) ; this%h2osoi_vol (:,:) = spval + allocate(this%excess_ice (begc:endc, 1:nlevgrnd)) ; this%excess_ice (:,:) = spval allocate(this%h2osfc (begc:endc)) ; this%h2osfc (:) = spval allocate(this%h2ocan (begc:endc)) ; this%h2ocan (:) = spval allocate(this%wslake_col (begc:endc)) ; this%wslake_col (:) = spval @@ -1417,7 +1418,6 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ if (use_fan) then allocate(this%h2osoi_tend_tsl_col(begc:endc)) ; this%h2osoi_tend_tsl_col(:) = spval end if - allocate(this%excess_ice (begc:endc)) ; this%excess_ice (:) = spval allocate(this%snw_rds_top (begc:endc)) ; this%snw_rds_top (:) = spval allocate(this%do_capsnow (begc:endc)) allocate(this%frac_sno (begc:endc)) ; this%frac_sno (:) = spval @@ -1473,9 +1473,14 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ ptr_col=this%h2osoi_ice, l2g_scale_type='veg') this%h2osoi_ice(begc:endc,:) = spval - call hist_addfld2d (fname='SOILICE_ICE', units='kg/m2', type2d='levgrnd', & - avgflag='A', long_name='soil ice (ice landunits only)', & - ptr_col=this%h2osoi_ice, l2g_scale_type='ice') + call hist_addfld2d (fname='SOILICE_ICE', units='kg/m2', type2d='levgrnd', & + avgflag='A', long_name='soil ice (ice landunits only)', & + ptr_col=this%h2osoi_ice, l2g_scale_type='ice') + + this%excess_ice(begc:endc, :) = spval + call hist_addfld2d (fname='EXCESS_ICE', units = '1', type2d='levgrnd', & + avgflag='A', long_name='Excess ground ice (0 to 1)', & + ptr_col=this%excess_ice, l2g_scale_type='veg') ! <- RPF: should this be natveg? this%h2osfc(begc:endc) = spval call hist_addfld1d (fname='H2OSFC', units='mm', & @@ -1565,11 +1570,6 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ default='inactive') end if - this%excess_ice(begc:endc) = spval - call hist_addfld1d (fname='EXCESS_ICE', units = '1', & - avgflag='A', long_name='Excess ground ice (0 to 1)', & - ptr_col=this%excess_ice, l2g_scale_type='veg') ! <- RPF: should this be natveg? - this%frac_sno(begc:endc) = spval call hist_addfld1d (fname='FSNO', units='1', & avgflag='A', long_name='fraction of ground covered by snow', & @@ -1641,7 +1641,6 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ this%frac_h2osfc_act(c) = 0._r8 this%h2orof(c) = 0._r8 this%frac_h2orof(c) = 0._r8 - this%excess_ice(c) = 0.5_r8 if (lun_pp%urbpoi(l)) then ! From Bonan 1996 (LSM technical note) @@ -1787,6 +1786,17 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ this%h2osoi_ice(c,j) = 0._r8 this%h2osoi_liq(c,j) = col_pp%dz(c,j)*denh2o*this%h2osoi_vol(c,j) endif + ! RPF - to do: a) apply to only polygonal ground + ! b) pull in ALT information to get starting point right. + ! for now, assuming no excess ice in top meter, 0.5 for 1-4 m, + ! 0.2 below 4 m. Also: need to check signs! + if (col_pp%z(c,j) > -1._r8) then + this%excess_ice(c,j) = 0._r8 + else if (col_pp%z(c,j) > -4._r8 .and. col_pp%z(c,j) < -1._r8) then + this%excess_ice(c,j) = 0.5_r8 + else + this%excess_ice(c,j) = 0.2_r8 + end if end do this%h2osoi_liq_old(c,:) = this%h2osoi_liq(c,:) @@ -1844,6 +1854,11 @@ subroutine col_ws_restart(this, bounds, ncid, flag, watsat_input) long_name='ice lens', units='kg/m2', & interpinic_flag='interp', readvar=readvar, data=this%h2osoi_ice) + call restartvar(ncid=ncid, flag=flag, varname='EXCESS_ICE', xtype=ncd_double, & + dim1name='column', dim2name='levgrnd', switchdim=.true., & + long_name='excess ground ice (0 to 1)', units='1', & + interpinic_flag='interp', readvar=readvar, data=this%excess_ice) + call restartvar(ncid=ncid, flag=flag, varname='SOILP', xtype=ncd_double, & dim1name='column', dim2name='levgrnd', switchdim=.true., & long_name='soil pressure ', units='Pa', & @@ -1907,11 +1922,6 @@ subroutine col_ws_restart(this, bounds, ncid, flag, watsat_input) this%snow_persistence(:) = 0.0_r8 end if - call restartvar(ncid=ncid, flag=flag, varname='EXCESS_ICE', xtype=ncd_double, & - dim1name='column', & - long_name='excess ground ice (0 to 1)', units='1', & - interpinic_flag='interp', readvar=readvar, data=this%excess_ice) - call restartvar(ncid=ncid, flag=flag, varname='frac_sno', xtype=ncd_double, & dim1name='column', & long_name='fraction of ground covered by snow (0 to 1)',units='1',&