Skip to content
33 changes: 15 additions & 18 deletions src/biogeophys/SoilTemperatureMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ module SoilTemperatureMod
! !PRIVATE MEMBER FUNCTIONS:
private :: SoilThermProp ! Set therm conduct. and heat cap of snow/soil layers
private :: PhaseChangeH2osfc ! When surface water freezes move ice to bottom snow layer
private :: PhaseChange_beta ! Calculation of the phase change within snow and soil layers
private :: PhaseChange ! Calculation of the phase change within snow and soil layers
private :: BuildingHAC ! Building Heating and Cooling for simpler method (introduced in CLM4.5)

real(r8), private, parameter :: thin_sfclayer = 1.0e-6_r8 ! Threshold for thin surface layer
Expand Down Expand Up @@ -517,7 +517,7 @@ subroutine SoilTemperature(bounds, num_urbanl, filter_urbanl, num_urbanc, filter
dhsdT(bounds%begc:bounds%endc), &
waterstatebulk_inst, waterdiagnosticbulk_inst, waterfluxbulk_inst, temperature_inst,energyflux_inst)

call Phasechange_beta (bounds, num_nolakec, filter_nolakec, &
call Phasechange (bounds, num_nolakec, filter_nolakec, &
dhsdT(bounds%begc:bounds%endc), &
soilstate_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, waterfluxbulk_inst, energyflux_inst, temperature_inst)

Expand Down Expand Up @@ -1130,7 +1130,7 @@ subroutine PhaseChangeH2osfc (bounds, num_nolakec, filter_nolakec, &
end subroutine PhaseChangeH2osfc

!-----------------------------------------------------------------------
subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &
subroutine Phasechange (bounds, num_nolakec, filter_nolakec, dhsdT, &
soilstate_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, waterfluxbulk_inst, energyflux_inst, temperature_inst)
!
! !DESCRIPTION:
Expand Down Expand Up @@ -1186,7 +1186,7 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &

!-----------------------------------------------------------------------

call t_startf( 'PhaseChangebeta' )
call t_startf( 'PhaseChange' )

! Enforce expected array sizes
SHR_ASSERT_ALL_FL((ubound(dhsdT) == (/bounds%endc/)), sourcefile, __LINE__)
Expand Down Expand Up @@ -1279,7 +1279,6 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &
! If ice exists above melt point, melt some to liquid.
if (h2osoi_ice(c,j) > 0._r8 .and. t_soisno(c,j) > tfrz) then
imelt(c,j) = 1
! tinc(c,j) = t_soisno(c,j) - tfrz
tinc(c,j) = tfrz - t_soisno(c,j)
t_soisno(c,j) = tfrz
endif
Expand All @@ -1288,7 +1287,6 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &
! If liquid exists below melt point, freeze some to ice.
if (h2osoi_liq(c,j) > 0._r8 .AND. t_soisno(c,j) < tfrz) then
imelt(c,j) = 2
! tinc(c,j) = t_soisno(c,j) - tfrz
tinc(c,j) = tfrz - t_soisno(c,j)
t_soisno(c,j) = tfrz
endif
Expand All @@ -1310,7 +1308,6 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &

if (h2osoi_ice(c,j) > 0. .AND. t_soisno(c,j) > tfrz) then
imelt(c,j) = 1
! tinc(c,j) = t_soisno(c,j) - tfrz
tinc(c,j) = tfrz - t_soisno(c,j)
t_soisno(c,j) = tfrz
endif
Expand All @@ -1334,7 +1331,6 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &

if (h2osoi_liq(c,j) > supercool(c,j) .AND. t_soisno(c,j) < tfrz) then
imelt(c,j) = 2
! tinc(c,j) = t_soisno(c,j) - tfrz
tinc(c,j) = tfrz - t_soisno(c,j)
t_soisno(c,j) = tfrz
endif
Expand All @@ -1343,7 +1339,6 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &
if (h2osno_no_layers(c) > 0._r8 .AND. j == 1) then
if (t_soisno(c,j) > tfrz) then
imelt(c,j) = 1
! tincc,j) = t_soisno(c,j) - tfrz
tinc(c,j) = tfrz - t_soisno(c,j)
t_soisno(c,j) = tfrz
endif
Expand Down Expand Up @@ -1438,14 +1433,16 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &
heatr = 0._r8
if (xm(c,j) > 0._r8) then !if there is excess heat to melt the ice
h2osoi_ice(c,j) = max(0._r8, wice0(c,j)-xm(c,j))
heatr = hm(c,j) - hfus*(wice0(c,j)-h2osoi_ice(c,j))/dtime
xm2(c,j) = xm(c,j) - h2osoi_ice(c,j) !excess ice melting
if (h2osoi_ice(c,j) == 0._r8) then ! this might be redundant
if (excess_ice(c,j) >= 0._r8 .and. xm2(c,j)>0._r8 .and. j>=2) then ! if there is excess ice to melt
excess_ice(c,j) = max(0._r8,wexice0(c,j) - xm2(c,j))
heatr = hm(c,j) - hfus * (wexice0(c,j)-excess_ice(c,j)+wice0(c,j)-h2osoi_ice(c,j)) / dtime
xm2(c,j) = xm(c,j) - wice0(c,j) ! Leftover melt
if (j>=1) then ! soil
if (excess_ice(c,j) >= 0._r8 .and. xm2(c,j)>0._r8) then ! if there is excess ice to melt
excess_ice(c,j) = max(0._r8,wexice0(c,j) - xm2(c,j))
endif
endif !end of excess ice block
heatr = hm(c,j) - hfus * (wexice0(c,j)-excess_ice(c,j)+ &
wice0(c,j)-h2osoi_ice(c,j)) / dtime
else !snow
heatr = hm(c,j) - hfus * (wice0(c,j)-h2osoi_ice(c,j)) / dtime
endif
else if (xm(c,j) < 0._r8) then
if (j <= 0) then
h2osoi_ice(c,j) = min(wmass0(c,j), wice0(c,j)-xm(c,j)) ! snow
Expand Down Expand Up @@ -1535,10 +1532,10 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &
end if
end do

call t_stopf( 'PhaseChangebeta' )
call t_stopf( 'PhaseChange' )
end associate

end subroutine Phasechange_beta
end subroutine Phasechange

!-----------------------------------------------------------------------
subroutine ComputeGroundHeatFluxAndDeriv(bounds, &
Expand Down