Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

soilsnow developments from AM3 #497

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
77 changes: 0 additions & 77 deletions src/coupled/AM3/control/cable/util/init/cbl_um_init_veg.F90.manual

This file was deleted.

16 changes: 8 additions & 8 deletions src/science/soilsnow/cbl_hyd_redistrib.F90
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,8 @@ MODULE hydraulic_redistribution_mod
SUBROUTINE hydraulic_redistribution(dels, soil, ssnow, canopy, veg, met)

USE cable_common_module, ONLY : wiltParam, satuParam
!data
USE cable_surface_types_mod, ONLY: evergreen_broadleaf, c4_grassland

IMPLICIT NONE
REAL, INTENT(IN) :: dels ! integration time step (s)
Expand Down Expand Up @@ -97,10 +99,9 @@ SUBROUTINE hydraulic_redistribution(dels, soil, ssnow, canopy, veg, met)
hr_perTime(:,k,j) = hr_perTime(:,k,j)/soil%zse(k)
hr_perTime(:,j,k) = hr_perTime(:,j,k)/soil%zse(j)

! Overwrite to give zero redistribution for all types except
! evergreen broadleaf (2) and c4 grass (7)
! NB: Hard-wired numbers should be removed in future version
WHERE( .NOT.(veg%iveg == 2 .OR. veg%iveg == 7 ) )
! Zero redistribution: all types except evergreen broadleaf, c4 grass
WHERE( .NOT. ( veg%iveg == evergreen_broadleaf .OR. &
veg%iveg == c4_grassland ) )
hr_perTime(:,k,j) = 0.0
hr_perTime(:,j,k) = 0.0
ENDWHERE
Expand Down Expand Up @@ -175,10 +176,9 @@ SUBROUTINE hydraulic_redistribution(dels, soil, ssnow, canopy, veg, met)
hr_perTime(:,k,j) = hr_perTime(:,k,j)/soil%zse(k)
hr_perTime(:,j,k) = hr_perTime(:,j,k)/soil%zse(j)

! Overwrite to give zero redistribution for all types except
! evergreen broadleaf (2) and c4 grass (7)
! NB: Hard-wired numbers should be removed in future version
WHERE( .NOT.( veg%iveg == 2 .OR. veg%iveg == 7 ) )
! Zero redistribution: all types except evergreen broadleaf, c4 grass
WHERE( .NOT. ( veg%iveg == evergreen_broadleaf .OR. &
veg%iveg == c4_grassland ) )
hr_perTime(:,k,j) = 0.0
hr_perTime(:,j,k) = 0.0
ENDWHERE
Expand Down
11 changes: 4 additions & 7 deletions src/science/soilsnow/cbl_remove_trans.F90
Original file line number Diff line number Diff line change
Expand Up @@ -31,14 +31,14 @@ SUBROUTINE remove_trans(dels, soil, ssnow, canopy, veg)
! which can be removed:
xx = canopy%fevc * dels / CHL * veg%froot(:,k) + diff(:,k-1) ! kg/m2
diff(:,k) = MAX( 0.0_r_2, ssnow%wb(:,k) - soil%swilt) & ! m3/m3
* soil%zse(k)*1000.0
* soil%zse(k)*Cdensity_liq
xxd = xx - diff(:,k)

WHERE ( xxd .GT. 0.0 )
ssnow%wb(:,k) = ssnow%wb(:,k) - diff(:,k) / (soil%zse(k)*1000.0)
ssnow%wb(:,k) = ssnow%wb(:,k) - diff(:,k) / (soil%zse(k)*Cdensity_liq)
diff(:,k) = xxd
ELSEWHERE
ssnow%wb(:,k) = ssnow%wb(:,k) - xx / (soil%zse(k)*1000.0)
ssnow%wb(:,k) = ssnow%wb(:,k) - xx / (soil%zse(k)*Cdensity_liq)
diff(:,k) = 0.0
ENDWHERE

Expand All @@ -52,10 +52,7 @@ SUBROUTINE remove_trans(dels, soil, ssnow, canopy, veg)
canopy%fevc = 0.0_r_2
END WHERE
DO k = 1,ms
ssnow%wb(:,k) = ssnow%wb(:,k) - ssnow%evapfbl(:,k)/(soil%zse(k)*1000.0)

! write(59,*) k, ssnow%wb(:,k), ssnow%evapfbl(:,k)/(soil%zse(k)*1000.0)
! write(59,*)
ssnow%wb(:,k) = ssnow%wb(:,k) - ssnow%evapfbl(:,k)/(soil%zse(k)*Cdensity_liq)
ENDDO


Expand Down
31 changes: 21 additions & 10 deletions src/science/soilsnow/cbl_smoisturev.F90
Original file line number Diff line number Diff line change
Expand Up @@ -60,9 +60,9 @@ SUBROUTINE smoisturev (dels,ssnow,soil,veg)
wbl_kp, & !
wh, & !
z3_k, & !
pwb_wbh !

pwb_wbh

REAL :: dfactor, sicemelt
REAL, DIMENSION(mp,ms+1) :: z1mult

! change dimension of at,bt,ct from 3*ms to ms (BP Jun2010)
Expand All @@ -83,10 +83,7 @@ SUBROUTINE smoisturev (dels,ssnow,soil,veg)
delt, & !
dtt !

LOGICAL :: is_open ! Is file open?
INTEGER :: &
u, & ! I/O unit
k
INTEGER :: i, k

at = 0.0
bt = 1.0
Expand Down Expand Up @@ -225,7 +222,7 @@ SUBROUTINE smoisturev (dels,ssnow,soil,veg)

END DO

ssnow%rnof2 = dels * REAL( fluxh(:,ms) ) * 1000.0
ssnow%rnof2 = dels * REAL( fluxh(:,ms) ) * Cdensity_liq

! wbh_k represents wblf(k-.5)
DO k = 2, ms
Expand Down Expand Up @@ -425,9 +422,23 @@ SUBROUTINE smoisturev (dels,ssnow,soil,veg)
CALL trimb(at, bt, ct, ssnow%wblf, ms)

DO k = 1, ms
ssatcurr(:,k) = soil%ssat - ssnow%wbice(:,k)
ssnow%wb(:,k) = ssnow%wblf(:,k) * ssatcurr(:,k) + ssnow%wbice(:,k)
ssnow%wbice(:,k) = MIN( ssnow%wbice(:,k), frozen_limit * ssnow%wb(:,k) )
ssatcurr(:,k) = soil%ssat - ssnow%wbice(:,k)
ssnow%wb(:,k) = ssnow%wblf(:,k) * ssatcurr(:,k) + ssnow%wbice(:,k)
END DO
! If redistribution of liquid has caused a layer to exceed frozen limit
! melt the excess and account for latent heat.
dfactor = 1.0 - Cdensity_ice / Cdensity_liq
DO k = 1, ms
DO i = 1, mp
IF (ssnow%wbice(i,k) > frozen_limit * ssnow%wb(i,k)) THEN
sicemelt = (ssnow%wbice(i,k) - frozen_limit * ssnow%wb(i,k)) / &
(1.0 - frozen_limit*dfactor)
ssnow%wbice(i,k) = ssnow%wbice(i,k) - sicemelt
ssnow%wb(i,k) = ssnow%wb(i,k) - dfactor*sicemelt
ssnow%tgg(i,k) = ssnow%tgg(i,k) - &
sicemelt*soil%zse(k)*Cdensity_ice*CHLF / REAL(ssnow%gammzz(i,k) )
END IF
END DO
END DO

END SUBROUTINE smoisturev
Expand Down
11 changes: 6 additions & 5 deletions src/science/soilsnow/cbl_snowAccum.F90
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,8 @@ SUBROUTINE snow_accum ( dels, canopy, met, ssnow, soil )

ssnow%tgg(i,1) = ssnow%tgg(i,1) + canopy%precis(i) * CHLF &
/ ( REAL( ssnow%gammzz(i,1) ) + Ccswat *canopy%precis(i) )
ssnow%dtmlt(i,1) = ssnow%dtmlt(i,1) + canopy%precis(i) * CHLF &
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

these are new terms from AM3 development. they are the counterpoint to the melt terms to %dtmlt and come about by the freezing of liquid precip onto frozen ground or snow pack.

good that they've been picked up - but not that surprising given they sit in a different part of the code base to the rest of the ice fixes

/ ( REAL( ssnow%gammzz(i,1) ) + Ccswat *canopy%precis(i) )
! change density due to water being added
ssnow%ssdn(i,1) = MIN( max_ssdn, MAX( 120.0, ssnow%ssdn(i,1) &
* ssnow%osnowd(i) / MAX( 0.01, ssnow%snowd(i) ) + Cdensity_liq &
Expand Down Expand Up @@ -87,6 +89,8 @@ SUBROUTINE snow_accum ( dels, canopy, met, ssnow, soil )

ssnow%tggsn(i,1) = ssnow%tggsn(i,1) + canopy%precis(i) * CHLF &
* osm(i) / (sgamm(i) * ssnow%osnowd(i) )
ssnow%dtmlt(i,1) = ssnow%dtmlt(i,1) + canopy%precis(i) * CHLF &
* osm(i) / (sgamm(i) * ssnow%osnowd(i) )
ssnow%smass(i,1) = ssnow%smass(i,1) + canopy%precis(i) &
* osm(i)/ssnow%osnowd(i)

Expand Down Expand Up @@ -173,11 +177,8 @@ SUBROUTINE snow_accum ( dels, canopy, met, ssnow, soil )
ssnow%sdepth(i,1) = MAX( 0.02, ssnow%smass(i,1) / ssnow%ssdn(i,1) )
ENDIF

canopy%segg(i) = 0.0

!INH: cls package
!we still need to conserve moisture/energy when evapsn is limited
!this is a key point of moisture non-conservation
! Use any remaining energy for evaporation of soil water
canopy%segg(i) = (CHL + CHLF) * (xxx(i) - ssnow%evapsn(i)) / CHL / dels

ENDIF

Expand Down
24 changes: 16 additions & 8 deletions src/science/soilsnow/cbl_soilfreeze.F90
Original file line number Diff line number Diff line change
Expand Up @@ -26,16 +26,22 @@ SUBROUTINE soilfreeze(dels, soil, ssnow,heat_cap_lower_limit)
IF(ssnow%tgg(i,k) < CTFRZ &
& .AND. frozen_limit*ssnow%wb(i,k) - ssnow%wbice(i,k) > .001) THEN

sicefreeze(i) = MIN( MAX( 0.0_r_2, ( frozen_limit * ssnow%wb(i,k) - &
ssnow%wbice(i,k) ) ) * soil%zse(k) * 1000.0, &
( CTFRZ - ssnow%tgg(i,k) ) * ssnow%gammzz(i,k) / CHLF )
! Limits on water mass that can be converted to ice due frozen_limit and saturation
sicefreeze(i) = MIN( frozen_limit * ssnow%wb(i,k) - ssnow%wbice(i,k) , &
(soil%ssat(i) - ssnow%wb(i,k)) / MAX((1.0-Cdensity_ice/Cdensity_liq),1.0E-3) )
! Apply limits on watermass to freeze by zero and the energy limit
sicefreeze(i) = MIN( MAX( 0.0_r_2, sicefreeze(i) ) * soil%zse(k) * Cdensity_ice, &
( CTFRZ - ssnow%tgg(i,k) ) * ssnow%gammzz(i,k) / CHLF )

ssnow%wbice(i,k) = MIN( ssnow%wbice(i,k) + sicefreeze(i) / (soil%zse(k) &
* 1000.0), frozen_limit * ssnow%wb(i,k) )
* Cdensity_ice), frozen_limit * ssnow%wb(i,k) )
ssnow%wb(i,k) = ssnow%wb(i,k) + sicefreeze(i) / (soil%zse(k) * Cdensity_ice) &
- sicefreeze(i) / (soil%zse(k) * Cdensity_liq)
xx(i) = soil%css(i) * soil%rhosoil(i)
max_arg1 = heat_cap_lower_limit(i,k)
max_arg2 = REAL((1.0 - soil%ssat(i)) * soil%css(i) * soil%rhosoil(i) ,r_2) &
+ (ssnow%wb(i,k) - ssnow%wbice(i,k)) * REAL(Ccswat * Cdensity_liq,r_2) &
+ ssnow%wbice(i,k) * REAL(Ccsice * Cdensity_liq * 0.9,r_2)
+ ssnow%wbice(i,k) * REAL(Ccsice * Cdensity_ice,r_2)
ssnow%gammzz(i,k) = MAX( max_arg1, max_arg2 ) * REAL( soil%zse(k),r_2 )

IF (k == 1 .AND. ssnow%isflag(i) == 0) &
Expand All @@ -46,16 +52,18 @@ SUBROUTINE soilfreeze(dels, soil, ssnow,heat_cap_lower_limit)

ELSEIF( ssnow%tgg(i,k) > CTFRZ .AND. ssnow%wbice(i,k) > 0. ) THEN

sicemelt(i) = MIN( ssnow%wbice(i,k) * soil%zse(k) * 1000.0, &
sicemelt(i) = MIN( ssnow%wbice(i,k) * soil%zse(k) * Cdensity_ice, &
( ssnow%tgg(i,k) - CTFRZ ) * ssnow%gammzz(i,k) / CHLF )

ssnow%wbice(i,k) = MAX( 0.0_r_2, ssnow%wbice(i,k) - sicemelt(i) &
/ (soil%zse(k) * 1000.0) )
/ (soil%zse(k) * Cdensity_ice) )
ssnow%wb(i,k) = ssnow%wb(i,k) - sicemelt(i) / (soil%zse(k) * Cdensity_ice) &
+ sicemelt(i) / (soil%zse(k) * Cdensity_liq)
xx(i) = soil%css(i) * soil%rhosoil(i)
max_arg1 = heat_cap_lower_limit(i,k)
max_arg2 = REAL((1.0 - soil%ssat(i)) * soil%css(i) * soil%rhosoil(i) ,r_2) &
+ (ssnow%wb(i,k) - ssnow%wbice(i,k)) * REAL(Ccswat * Cdensity_liq,r_2) &
+ ssnow%wbice(i,k) * REAL(Ccsice * Cdensity_liq * 0.9,r_2)
+ ssnow%wbice(i,k) * REAL(Ccsice * Cdensity_ice,r_2)

ssnow%gammzz(i,k) = MAX( max_arg1, max_arg2 ) * REAL( soil%zse(k),r_2 )

Expand Down
1 change: 1 addition & 0 deletions src/science/soilsnow/cbl_soilsnow_data.F90
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ MODULE cbl_ssnow_data_mod
USE cable_phys_constants_mod, ONLY : CHL => HL
USE cable_phys_constants_mod, ONLY : CHLF => HLF
USE cable_phys_constants_mod, ONLY : Cdensity_liq => density_liq
USE cable_phys_constants_mod, ONLY : Cdensity_ice => density_ice
USE cable_phys_constants_mod, ONLY : Ccgsnow => cgsnow
USE cable_phys_constants_mod, ONLY : Ccswat => cswat
USE cable_phys_constants_mod, ONLY : Ccsice => csice
Expand Down
16 changes: 8 additions & 8 deletions src/science/soilsnow/cbl_soilsnow_init_special.F90
Original file line number Diff line number Diff line change
Expand Up @@ -64,25 +64,25 @@ SUBROUTINE spec_init_soil_snow(dels, soil, ssnow, canopy, met, bal, veg,heat_cap
ssnow%wb(:,4) = 0.95 * soil%ssat
ssnow%wb(:,5) = 0.95 * soil%ssat
ssnow%wb(:,6) = 0.95 * soil%ssat
ssnow%wbice(:,1) = 0.90 * ssnow%wb(:,1)
ssnow%wbice(:,2) = 0.90 * ssnow%wb(:,2)
ssnow%wbice(:,3) = 0.90 * ssnow%wb(:,3)
ssnow%wbice(:,4) = 0.90 * ssnow%wb(:,4)
ssnow%wbice(:,5) = 0.90 * ssnow%wb(:,5)
ssnow%wbice(:,6) = 0.90 * ssnow%wb(:,6)
ssnow%wbice(:,1) = frozen_limit * ssnow%wb(:,1)
ssnow%wbice(:,2) = frozen_limit * ssnow%wb(:,2)
ssnow%wbice(:,3) = frozen_limit * ssnow%wb(:,3)
ssnow%wbice(:,4) = frozen_limit * ssnow%wb(:,4)
ssnow%wbice(:,5) = frozen_limit * ssnow%wb(:,5)
ssnow%wbice(:,6) = frozen_limit * ssnow%wb(:,6)
ENDWHERE
xx=REAL(heat_cap_lower_limit(:,1))
ssnow%gammzz(:,1) = MAX( (1.0 - soil%ssat) * soil%css * soil%rhosoil &
& + (ssnow%wb(:,1) - ssnow%wbice(:,1) ) * Ccswat * Cdensity_liq &
& + ssnow%wbice(:,1) * Ccsice * Cdensity_liq * .9, xx ) * soil%zse(1)
& + ssnow%wbice(:,1) * Ccsice * Cdensity_ice, xx ) * soil%zse(1)
END IF
ENDIF ! if(.NOT.cable_runtime_coupled)

IF (ktau <= 1) THEN
xx=heat_cap_lower_limit(:,1)
ssnow%gammzz(:,1) = MAX( (1.0 - soil%ssat) * soil%css * soil%rhosoil &
& + (ssnow%wb(:,1) - ssnow%wbice(:,1) ) * Ccswat * Cdensity_liq &
& + ssnow%wbice(:,1) * Ccsice * Cdensity_liq * .9, xx ) * soil%zse(1) + &
& + ssnow%wbice(:,1) * Ccsice * Cdensity_ice, xx ) * soil%zse(1) + &
& (1. - ssnow%isflag) * Ccgsnow * ssnow%snowd
END IF

Expand Down
Loading
Loading