diff --git a/components/elm/src/biogeophys/BalanceCheckMod.F90 b/components/elm/src/biogeophys/BalanceCheckMod.F90 index f46091972852..43c0fe58aa46 100644 --- a/components/elm/src/biogeophys/BalanceCheckMod.F90 +++ b/components/elm/src/biogeophys/BalanceCheckMod.F90 @@ -469,6 +469,27 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, & snow_sinks(c) = qflx_sub_snow(c) + qflx_evap_grnd(c) + qflx_snow_melt(c) & + qflx_snwcp_ice(c) + qflx_snwcp_liq(c) + qflx_sl_top_soil(c) + if (lun_pp%itype(l) == istice .or. lun_pp%itype(l) == istice_mec) then + ! With fractional frac_sno_eff on glacier columns, only the + ! snow-covered fraction receives dew/rain and loses sub/evap. + if (do_capsnow(c) .and. .not. use_firn_percolation_and_compaction) then + snow_sources(c) = qflx_prec_grnd(c) & + + frac_sno_eff(c) * (qflx_dew_snow(c) + qflx_dew_grnd(c)) + + snow_sinks(c) = frac_sno_eff(c) * (qflx_sub_snow(c) + qflx_evap_grnd(c)) & + + qflx_snwcp_ice(c) + qflx_snwcp_liq(c) & + + qflx_snow_melt(c) + qflx_sl_top_soil(c) + else + snow_sources(c) = qflx_snow_grnd_col(c) & + + frac_sno_eff(c) * (qflx_rain_grnd_col(c) & + + qflx_dew_snow(c) + qflx_dew_grnd(c)) + + snow_sinks(c) = frac_sno_eff(c) * (qflx_sub_snow(c) + qflx_evap_grnd(c)) & + + qflx_snwcp_ice(c) + qflx_snwcp_liq(c) & + + qflx_snow_melt(c) + qflx_sl_top_soil(c) + endif + endif + if (col_pp%is_lake(c)) then if (.not. use_firn_percolation_and_compaction) then if ( do_capsnow(c)) then diff --git a/components/elm/src/biogeophys/CanopyHydrologyMod.F90 b/components/elm/src/biogeophys/CanopyHydrologyMod.F90 index 3d2751f88cd9..6a95def159b2 100644 --- a/components/elm/src/biogeophys/CanopyHydrologyMod.F90 +++ b/components/elm/src/biogeophys/CanopyHydrologyMod.F90 @@ -642,14 +642,15 @@ subroutine CanopyHydrology(bounds, & end if !end of do_capsnow construct ! set frac_sno_eff variable - if (col_pp%is_soil(c) .or. col_pp%is_crop(c)) then - if (subgridflag ==1) then - frac_sno_eff(c) = frac_sno(c) - else - frac_sno_eff(c) = 1._r8 - endif - else + ! Port of CLM fractional-snow energy fix: when subgrid fluxes are on, + ! allow fractional frac_sno_eff for all non-urban, non-lake columns + ! (including glacier and wetland), matching soil/crop treatment. + if (subgridflag == 1 .and. .not. urbpoi(l) .and. ltype(l) /= istdlak) then + frac_sno_eff(c) = frac_sno(c) + else if (frac_sno(c) > 0._r8) then frac_sno_eff(c) = 1._r8 + else + frac_sno_eff(c) = 0._r8 endif if (ltype(l) == istwet .and. t_grnd(c)>tfrz) then diff --git a/components/elm/src/biogeophys/SnowHydrologyMod.F90 b/components/elm/src/biogeophys/SnowHydrologyMod.F90 index da648f88ca2c..44636b0ab8da 100644 --- a/components/elm/src/biogeophys/SnowHydrologyMod.F90 +++ b/components/elm/src/biogeophys/SnowHydrologyMod.F90 @@ -568,6 +568,7 @@ subroutine SnowCompaction(bounds, num_snowc, filter_snowc, & ! ! !LOCAL VARIABLES: integer :: j, l, c, fc, t ! indices + ! (urbpoi / lakpoi used for CLM-consistent fractional melt compaction gate) ! parameters real(r8), parameter :: c2 = 23.e-3_r8 ! [m3/kg] real(r8), parameter :: c3 = 2.777e-6_r8 ! [1/s] @@ -605,6 +606,8 @@ subroutine SnowCompaction(bounds, num_snowc, filter_snowc, & snl => col_pp%snl , & ! Input: [integer (:) ] number of snow layers n_melt => col_pp%n_melt , & ! Input: [real(r8) (:) ] SCA shape parameter ltype => lun_pp%itype , & ! Input: [integer (:) ] landunit type + urbpoi => lun_pp%urbpoi , & ! Input: [logical (:) ] true => landunit is an urban point + lakpoi => lun_pp%lakpoi , & ! Input: [logical (:) ] true => landunit is a lake point forc_wind => top_as_inst%windbot , & ! Input: [real(r8) (:) ] atmospheric wind speed (m/s) t_soisno => col_es%t_soisno , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) @@ -703,9 +706,13 @@ subroutine SnowCompaction(bounds, num_snowc, filter_snowc, & endif ! Compaction occurring during melt + ! Port of CLM fractional-snow melt compaction fix: use the + ! fractional ddz3 form for all non-lake, non-urban columns when + ! subgridflag is on (includes glacier and wetland). if (imelt(c,j) == 1) then - if(subgridflag==1 .and. (col_pp%is_soil(c) .or. col_pp%is_crop(c))) then + l = col_pp%landunit(c) + if (subgridflag == 1 .and. (.not. lakpoi(l) .and. .not. urbpoi(l))) then ! first term is delta mass over mass ddz3 = max(0._r8,min(1._r8,(swe_old(c,j) - wx)/wx)) diff --git a/components/elm/src/biogeophys/SoilHydrologyMod.F90 b/components/elm/src/biogeophys/SoilHydrologyMod.F90 index 51f1dd28711f..3283dbc47f1d 100644 --- a/components/elm/src/biogeophys/SoilHydrologyMod.F90 +++ b/components/elm/src/biogeophys/SoilHydrologyMod.F90 @@ -554,23 +554,23 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f !5. surface runoff from h2osfc if (h2osfcflag==1) then ! calculate runoff from h2osfc ------------------------------------- - if (use_modified_infil) then + !if (use_modified_infil) then if (frac_h2osfc_act(c) <= pc .and. frac_h2osfc(c) <= pc) then frac_infclust=0.0_r8 else - if (frac_h2osfc(c) <= pc) then + if (frac_h2osfc_act(c) <= pc) then frac_infclust=(frac_h2osfc_act(c)-pc)**mu else - frac_infclust=(frac_h2osfc(c)-pc)**mu + frac_infclust=(frac_h2osfc_act(c)-pc)**mu endif endif - else - if (frac_h2osfc(c) <= pc) then - frac_infclust=0.0_r8 - else - frac_infclust=(frac_h2osfc(c)-pc)**mu - endif - endif + !else + ! if (frac_h2osfc(c) <= pc) then + ! frac_infclust=0.0_r8 + ! else + ! frac_infclust=(frac_h2osfc(c)-pc)**mu + ! endif + !endif endif if (lun_pp%ispolygon(col_pp%landunit(c))) then @@ -596,7 +596,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f ! limit runoff to value of storage above S(pc) if(h2osfc(c) >= h2osfc_thresh(c) .and. h2osfcflag/=0) then ! spatially variable k_wet - k_wet=1.0_r8 * sin((rpi/180._r8) * col_pp%topo_slope(c)) + k_wet=1.0e-4_r8 * sin((rpi/180._r8) * max(col_pp%topo_slope(c), 1.0e-3_r8)) qflx_h2osfc_surf(c) = k_wet * frac_infclust * (h2osfc(c) - h2osfc_thresh(c)) qflx_h2osfc_surf(c)=min(qflx_h2osfc_surf(c),(h2osfc(c) - h2osfc_thresh(c))/dtime) diff --git a/components/elm/src/biogeophys/SoilTemperatureMod.F90 b/components/elm/src/biogeophys/SoilTemperatureMod.F90 index bea06e4e89eb..555c2c4966ab 100644 --- a/components/elm/src/biogeophys/SoilTemperatureMod.F90 +++ b/components/elm/src/biogeophys/SoilTemperatureMod.F90 @@ -1044,6 +1044,10 @@ subroutine SoilThermProp (bounds, num_nolakec, filter_nolakec, & .AND. col_pp%itype(c) /= icol_sunwall .AND. col_pp%itype(c) /= icol_shadewall .AND. & col_pp%itype(c) /= icol_roof) then cv(c,j) = csol(c,j)*(1._r8-watsat(c,j))*dz(c,j) + (h2osoi_ice(c,j)*cpice + h2osoi_liq(c,j)*cpliq) + ! Port of CLM bedrock heat-capacity fix: do not reduce solid-rock + ! heat capacity by soil porosity below the soil-bedrock boundary. + ! csol is already set to bedrock mineral heat capacity for these layers. + if (j > nlevbed) cv(c,j) = csol(c,j)*dz(c,j) else if (lun_pp%itype(l) == istwet) then cv(c,j) = (h2osoi_ice(c,j)*cpice + h2osoi_liq(c,j)*cpliq) if (j > nlevbed) cv(c,j) = csol(c,j)*dz(c,j) diff --git a/components/elm/src/main/accumulMod.F90 b/components/elm/src/main/accumulMod.F90 index a3e2185b1ed7..51ef4f108d44 100644 --- a/components/elm/src/main/accumulMod.F90 +++ b/components/elm/src/main/accumulMod.F90 @@ -18,10 +18,17 @@ module accumulMod ! continuously accumulated. The trigger value "-99999." resets ! the accumulation to zero. ! + ! Inactive / zero-weight subgrid elements are skipped during update, + ! following the CLM (ctsm5.4.042) active-element accumulator fix. + ! ! !USES: use shr_kind_mod, only: r8 => shr_kind_r8 use shr_sys_mod , only: shr_sys_abort use elm_varctl , only: iulog + use TopounitType , only : top_pp + use LandunitType , only : lun_pp + use ColumnType , only : col_pp + use VegetationType , only : veg_pp ! ! !PUBLIC TYPES: implicit none @@ -56,8 +63,11 @@ module accumulMod integer :: end1d !subgrid type ending index integer :: num1d !total subgrid points integer :: numlev !number of vertical levels in field + logical, pointer :: active(:) => null() ! whether each point is active (null => all active) + logical, pointer :: active_storage(:) => null() ! owned storage when no type-level active flag exists real(r8) :: initval !initial value of accumulated field - real(r8), pointer :: val(:,:) !accumulated field + real(r8), pointer :: val(:,:) => null() !accumulated field + integer , pointer :: nsteps(:,:) => null() ! active accumulation step counts since last reset integer :: period !field accumulation period (in model time steps) end type accum_field @@ -154,22 +164,30 @@ subroutine init_accum_field (name, units, desc, & beg1d = begg end1d = endg num1d = numg + ! Gridcells have no type-level active flag; treat all as active. + allocate(accum(nf)%active_storage(beg1d:end1d)) + accum(nf)%active_storage(:) = .true. + accum(nf)%active => accum(nf)%active_storage case ('topounit') beg1d = begt end1d = endt num1d = numt + accum(nf)%active => top_pp%active case ('landunit') beg1d = begl end1d = endl num1d = numl + accum(nf)%active => lun_pp%active case ('column') beg1d = begc end1d = endc num1d = numc + accum(nf)%active => col_pp%active case ('pft') beg1d = begp end1d = endp num1d = nump + accum(nf)%active => veg_pp%active case default write(iulog,*)'ACCUMULINIT: unknown subgrid type ',subgrid_type call shr_sys_abort () @@ -190,7 +208,10 @@ subroutine init_accum_field (name, units, desc, & ! Allocate and initialize accumulation field allocate(accum(nf)%val(beg1d:end1d,numlev)) - accum(nf)%val(beg1d:end1d,numlev) = init_value + accum(nf)%val(beg1d:end1d,1:numlev) = init_value + + allocate(accum(nf)%nsteps(beg1d:end1d,numlev)) + accum(nf)%nsteps(beg1d:end1d,1:numlev) = 0 end subroutine init_accum_field @@ -378,6 +399,7 @@ subroutine update_accum_field_sl (name, field, nstep) ! !DESCRIPTION: ! Accumulate single level field over specified time interval. ! The appropriate field is accumulated in the array [accval]. + ! Values of 'field' are ignored at inactive points. ! ! !ARGUMENTS: implicit none @@ -389,6 +411,7 @@ subroutine update_accum_field_sl (name, field, nstep) integer :: i,k,nf !indices integer :: accper !temporary accumulation period integer :: beg,end !subgrid beginning,ending indices + logical :: time_to_reset !------------------------------------------------------------------------ ! find field index. return if "name" is not on list @@ -433,33 +456,54 @@ subroutine update_accum_field_sl (name, field, nstep) if (accum(nf)%acctype == 'timeavg') then !time average field reset every accumulation period - !normalize at end of accumulation period + !normalize at end of accumulation period using active-point nsteps - if ((mod(nstep,accum(nf)%period) == 1 .or. accum(nf)%period == 1) .and. (nstep /= 0))then - accum(nf)%val(beg:end,1) = 0._r8 - end if - accum(nf)%val(beg:end,1) = accum(nf)%val(beg:end,1) + field(beg:end) - if (mod(nstep,accum(nf)%period) == 0) then - accum(nf)%val(beg:end,1) = accum(nf)%val(beg:end,1) / accum(nf)%period - endif + time_to_reset = (mod(nstep,accum(nf)%period) == 1 .or. accum(nf)%period == 1) .and. (nstep /= 0) + do k = beg,end + if (accum(nf)%active(k)) then + if (time_to_reset) then + accum(nf)%val(k,1) = 0._r8 + accum(nf)%nsteps(k,1) = 0 + end if + accum(nf)%val(k,1) = accum(nf)%val(k,1) + field(k) + accum(nf)%nsteps(k,1) = accum(nf)%nsteps(k,1) + 1 + if (mod(nstep,accum(nf)%period) == 0) then + if (accum(nf)%nsteps(k,1) > 0) then + accum(nf)%val(k,1) = accum(nf)%val(k,1) / real(accum(nf)%nsteps(k,1),r8) + else + accum(nf)%val(k,1) = accum(nf)%initval + end if + end if + end if + end do else if (accum(nf)%acctype == 'runmean') then - !running mean - reset accumulation period until greater than nstep + !running mean - only update active points; cap nsteps at period - accper = min (nstep,accum(nf)%period) - accum(nf)%val(beg:end,1) = ((accper-1)*accum(nf)%val(beg:end,1) + field(beg:end)) / accper + do k = beg,end + if (accum(nf)%active(k)) then + accum(nf)%nsteps(k,1) = min(accum(nf)%nsteps(k,1) + 1, accum(nf)%period) + accper = accum(nf)%nsteps(k,1) + accum(nf)%val(k,1) = ((accper-1)*accum(nf)%val(k,1) + field(k)) / real(accper,r8) + end if + end do else if (accum(nf)%acctype == 'runaccum') then !running accumulation field reset at trigger -99999 do k = beg,end - if (nint(field(k)) == -99999) then - accum(nf)%val(k,1) = 0._r8 + if (accum(nf)%active(k)) then + if (nint(field(k)) == -99999) then + accum(nf)%val(k,1) = 0._r8 + accum(nf)%nsteps(k,1) = 0 + else + accum(nf)%val(k,1) = min(max(accum(nf)%val(k,1) + field(k), 0._r8), 99999._r8) + accum(nf)%nsteps(k,1) = accum(nf)%nsteps(k,1) + 1 + end if end if end do - accum(nf)%val(beg:end,1) = min(max(accum(nf)%val(beg:end,1) + field(beg:end), 0._r8), 99999._r8) end if @@ -470,6 +514,7 @@ subroutine update_accum_field_ml (name, field, nstep) ! ! !DESCRIPTION: ! Accumulate multi level field over specified time interval. + ! Values of 'field' are ignored at inactive points. ! ! !ARGUMENTS: implicit none @@ -482,6 +527,7 @@ subroutine update_accum_field_ml (name, field, nstep) integer :: accper !temporary accumulation period integer :: beg,end !subgrid beginning,ending indices integer :: numlev !number of vertical levels + logical :: time_to_reset !------------------------------------------------------------------------ ! find field index. return if "name" is not on list @@ -532,38 +578,54 @@ subroutine update_accum_field_ml (name, field, nstep) if (accum(nf)%acctype == 'timeavg') then - !time average field reset every accumulation period - !normalize at end of accumulation period - - if ((mod(nstep,accum(nf)%period) == 1 .or. accum(nf)%period == 1) .and. (nstep /= 0))then - accum(nf)%val(beg:end,1:numlev) = 0._r8 - endif - accum(nf)%val(beg:end,1:numlev) = accum(nf)%val(beg:end,1:numlev) + field(beg:end,1:numlev) - if (mod(nstep,accum(nf)%period) == 0) then - accum(nf)%val(beg:end,1:numlev) = accum(nf)%val(beg:end,1:numlev) / accum(nf)%period - endif + time_to_reset = (mod(nstep,accum(nf)%period) == 1 .or. accum(nf)%period == 1) .and. (nstep /= 0) + do j = 1,numlev + do k = beg,end + if (accum(nf)%active(k)) then + if (time_to_reset) then + accum(nf)%val(k,j) = 0._r8 + accum(nf)%nsteps(k,j) = 0 + end if + accum(nf)%val(k,j) = accum(nf)%val(k,j) + field(k,j) + accum(nf)%nsteps(k,j) = accum(nf)%nsteps(k,j) + 1 + if (mod(nstep,accum(nf)%period) == 0) then + if (accum(nf)%nsteps(k,j) > 0) then + accum(nf)%val(k,j) = accum(nf)%val(k,j) / real(accum(nf)%nsteps(k,j),r8) + else + accum(nf)%val(k,j) = accum(nf)%initval + end if + end if + end if + end do + end do else if (accum(nf)%acctype == 'runmean') then - !running mean - reset accumulation period until greater than nstep - - accper = min (nstep,accum(nf)%period) - accum(nf)%val(beg:end,1:numlev) = & - ((accper-1)*accum(nf)%val(beg:end,1:numlev) + field(beg:end,1:numlev)) / accper + do j = 1,numlev + do k = beg,end + if (accum(nf)%active(k)) then + accum(nf)%nsteps(k,j) = min(accum(nf)%nsteps(k,j) + 1, accum(nf)%period) + accper = accum(nf)%nsteps(k,j) + accum(nf)%val(k,j) = ((accper-1)*accum(nf)%val(k,j) + field(k,j)) / real(accper,r8) + end if + end do + end do else if (accum(nf)%acctype == 'runaccum') then - !running accumulation field reset at trigger -99999 - do j = 1,numlev do k = beg,end - if (nint(field(k,j)) == -99999) then - accum(nf)%val(k,j) = 0._r8 + if (accum(nf)%active(k)) then + if (nint(field(k,j)) == -99999) then + accum(nf)%val(k,j) = 0._r8 + accum(nf)%nsteps(k,j) = 0 + else + accum(nf)%val(k,j) = min(max(accum(nf)%val(k,j) + field(k,j), 0._r8), 99999._r8) + accum(nf)%nsteps(k,j) = accum(nf)%nsteps(k,j) + 1 + end if end if end do end do - accum(nf)%val(beg:end,1:numlev) = & - min(max(accum(nf)%val(beg:end,1:numlev) + field(beg:end,1:numlev), 0._r8), 99999._r8) end if @@ -612,6 +674,26 @@ subroutine accumulRest( ncid, flag ) data=accum(nf)%val, readvar=readvar) end if + varname = trim(accum(nf)%name) // '_NSTEPS' + if (accum(nf)%numlev == 1) then + call restartvar(ncid=ncid, flag=flag, varname=varname, xtype=ncd_int, & + dim1name=accum(nf)%type1d, & + long_name='number of accumulated steps for '//trim(accum(nf)%name), & + units='-', & + interpinic_flag='interp', & + data=accum(nf)%nsteps, readvar=readvar) + else + call restartvar(ncid=ncid, flag=flag, varname=varname, xtype=ncd_int, & + dim1name=accum(nf)%type1d, dim2name=accum(nf)%type2d, & + long_name='number of accumulated steps for '//trim(accum(nf)%name), & + units='-', & + interpinic_flag='interp', & + data=accum(nf)%nsteps, readvar=readvar) + end if + if (flag == 'read' .and. .not. readvar) then + accum(nf)%nsteps(accum(nf)%beg1d:accum(nf)%end1d,1:accum(nf)%numlev) = 0 + end if + varname = trim(accum(nf)%name) // '_PERIOD' call restartvar(ncid=ncid, flag=flag, varname=varname, xtype=ncd_int, & long_name=trim(accum(nf)%desc), units='time steps', &