From 8fe05fd23ca3bff38397d9df62019121299ccebb Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Wed, 27 Nov 2024 17:22:40 +1100 Subject: [PATCH 01/16] developments from AM3 --- src/science/soilsnow/cbl_hyd_redistrib.F90 | 16 +++++----- src/science/soilsnow/cbl_remove_trans.F90 | 11 +++---- src/science/soilsnow/cbl_smoisturev.F90 | 31 +++++++++++++------ src/science/soilsnow/cbl_snowAccum.F90 | 7 ++--- src/science/soilsnow/cbl_soilfreeze.F90 | 24 +++++++++----- src/science/soilsnow/cbl_soilsnow_data.F90 | 1 + .../soilsnow/cbl_soilsnow_init_special.F90 | 16 +++++----- src/science/soilsnow/cbl_soilsnow_main.F90 | 24 +++++++------- src/science/soilsnow/cbl_stempv.F90 | 7 ++--- src/science/soilsnow/cbl_surfbv.F90 | 24 +++++++------- src/science/soilsnow/cbl_thermal.F90 | 12 ++----- 11 files changed, 90 insertions(+), 83 deletions(-) diff --git a/src/science/soilsnow/cbl_hyd_redistrib.F90 b/src/science/soilsnow/cbl_hyd_redistrib.F90 index 0c63f61e3..ad63b8ef3 100644 --- a/src/science/soilsnow/cbl_hyd_redistrib.F90 +++ b/src/science/soilsnow/cbl_hyd_redistrib.F90 @@ -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) @@ -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 @@ -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 diff --git a/src/science/soilsnow/cbl_remove_trans.F90 b/src/science/soilsnow/cbl_remove_trans.F90 index 6274cedbc..50a438ad8 100644 --- a/src/science/soilsnow/cbl_remove_trans.F90 +++ b/src/science/soilsnow/cbl_remove_trans.F90 @@ -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 @@ -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 diff --git a/src/science/soilsnow/cbl_smoisturev.F90 b/src/science/soilsnow/cbl_smoisturev.F90 index ac38c482a..635e4f06b 100644 --- a/src/science/soilsnow/cbl_smoisturev.F90 +++ b/src/science/soilsnow/cbl_smoisturev.F90 @@ -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) @@ -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 @@ -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 @@ -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 diff --git a/src/science/soilsnow/cbl_snowAccum.F90 b/src/science/soilsnow/cbl_snowAccum.F90 index 28b6a3dd6..55af4ac17 100644 --- a/src/science/soilsnow/cbl_snowAccum.F90 +++ b/src/science/soilsnow/cbl_snowAccum.F90 @@ -173,11 +173,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 diff --git a/src/science/soilsnow/cbl_soilfreeze.F90 b/src/science/soilsnow/cbl_soilfreeze.F90 index ec068b2b0..a1b1a3cde 100644 --- a/src/science/soilsnow/cbl_soilfreeze.F90 +++ b/src/science/soilsnow/cbl_soilfreeze.F90 @@ -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) & @@ -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 ) diff --git a/src/science/soilsnow/cbl_soilsnow_data.F90 b/src/science/soilsnow/cbl_soilsnow_data.F90 index a3a34060c..88ac8adc5 100644 --- a/src/science/soilsnow/cbl_soilsnow_data.F90 +++ b/src/science/soilsnow/cbl_soilsnow_data.F90 @@ -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 diff --git a/src/science/soilsnow/cbl_soilsnow_init_special.F90 b/src/science/soilsnow/cbl_soilsnow_init_special.F90 index 816ca797b..ea6fafbd2 100644 --- a/src/science/soilsnow/cbl_soilsnow_init_special.F90 +++ b/src/science/soilsnow/cbl_soilsnow_init_special.F90 @@ -64,17 +64,17 @@ 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) @@ -82,7 +82,7 @@ SUBROUTINE spec_init_soil_snow(dels, soil, ssnow, canopy, met, bal, veg,heat_cap 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 diff --git a/src/science/soilsnow/cbl_soilsnow_main.F90 b/src/science/soilsnow/cbl_soilsnow_main.F90 index b7c100480..1ed8af8ae 100644 --- a/src/science/soilsnow/cbl_soilsnow_main.F90 +++ b/src/science/soilsnow/cbl_soilsnow_main.F90 @@ -40,6 +40,8 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) USE snow_accum_mod, ONLY: snow_accum USE snowdensity_mod, ONLY: snowDensity +USE cable_phys_constants_mod, ONLY: density_liq, density_ice + REAL, INTENT(IN) :: dels ! integration time step (s) TYPE(soil_parameter_type), INTENT(INOUT) :: soil TYPE(soil_snow_type), INTENT(INOUT) :: ssnow @@ -51,21 +53,20 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) REAL, DIMENSION(mp) :: snowmlt REAL, DIMENSION(mp) :: totwet REAL, DIMENSION(mp) :: weting - REAL, DIMENSION(mp) :: xx + REAL(r_2), DIMENSION(mp) :: xx REAL(r_2), DIMENSION(mp) :: xxx REAL(r_2), DIMENSION(mp) :: deltat,sinfil1,sinfil2,sinfil3 REAL :: zsetot INTEGER, SAVE :: ktau =0 -REAL :: wbliq(mp,ms) ktau = ktau +1 !this is the value it is initialized with in cable_common anyway max_glacier_snowd = 1100.0 ! for ACCESS1.3 onwards. = 50000.0 for ACCESS1.0 zsetot = SUM(soil%zse) - ssnow%tggav = 0. + ssnow%tggav(:) = 0. DO k = 1, ms - ssnow%tggav = ssnow%tggav + soil%zse(k)*ssnow%tgg(:,k)/zsetot + ssnow%tggav(:) = ssnow%tggav(:) + ( (soil%zse(k)/zsetot) * ssnow%tgg(:,k) ) soil%heat_cap_lower_limit(:,k) = MAX( 0.01, soil%css(:) * soil%rhosoil(:) ) END DO @@ -83,8 +84,7 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) ssnow%dtmlt = 0.0 ssnow%osnowd = ssnow%snowd - - wbliq = ssnow%wb - ssnow%wbice + ssnow%wbliq = ssnow%wb - ssnow%wbice !%cable_runtime_coupled special initalizations in um_init NA for ESM1.5 @@ -92,7 +92,7 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) IF (ktau <= 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 @@ -122,7 +122,7 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) ! snow aging etc... CALL snowl_adjust(dels, ssnow, canopy ) - CALL stempv(dels, canopy, ssnow, soil, soil%heat_cap_lower_limit ) + CALL stempv(dels, canopy, ssnow, soil, REAL(soil%heat_cap_lower_limit) ) ssnow%tss = (1-ssnow%isflag)*ssnow%tgg(:,1) + ssnow%isflag*ssnow%tggsn(:,1) @@ -133,7 +133,7 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) CALL remove_trans(dels, soil, ssnow, canopy, veg) - CALL soilfreeze(dels, soil, ssnow, soil%heat_cap_lower_limit) + CALL soilfreeze(dels, soil, ssnow, REAL(soil%heat_cap_lower_limit) ) totwet = canopy%precis + ssnow%smelt @@ -191,11 +191,13 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) ! Set weighted soil/snow surface temperature ssnow%tss=(1-ssnow%isflag)*ssnow%tgg(:,1) + ssnow%isflag*ssnow%tggsn(:,1) - wbliq = ssnow%wb - ssnow%wbice + ssnow%wbliq = ssnow%wb - ssnow%wbice ssnow%wbtot = 0.0 DO k = 1, ms - ssnow%wbtot = ssnow%wbtot + REAL(ssnow%wb(:,k)*1000.0*soil%zse(k),r_2) + ! tot moisture this timestep + ssnow%wbtot(:) = ssnow%wbtot(:) + & + (ssnow%wbliq(:,k)*density_liq + ssnow%wbice(:,k)*density_ice) * soil%zse(k) END DO diff --git a/src/science/soilsnow/cbl_stempv.F90 b/src/science/soilsnow/cbl_stempv.F90 index 2111c0d9a..cd5a280b2 100644 --- a/src/science/soilsnow/cbl_stempv.F90 +++ b/src/science/soilsnow/cbl_stempv.F90 @@ -90,7 +90,7 @@ SUBROUTINE stempv(dels, canopy, ssnow, soil,heat_cap_lower_limit) ssnow%gammzz(:,k) = MAX( heat_cap_lower_limit(:,k), & ( 1.0 - soil%ssat ) * soil%css * soil%rhosoil & + soil%ssat * ( wblfsp * Ccswat * Cdensity_liq + & - ssnow%wbfice(:,k) * Ccsice * Cdensity_liq * 0.9 ) ) & + ssnow%wbfice(:,k) * Ccsice * Cdensity_ice ) ) & * soil%zse(k) ssnow%gammzz(:,k) = ssnow%gammzz(:,k) + Ccgsnow * ssnow%snowd @@ -113,7 +113,7 @@ SUBROUTINE stempv(dels, canopy, ssnow, soil,heat_cap_lower_limit) ssnow%gammzz(:,k) = MAX( REAL(heat_cap_lower_limit(:,k)), & ( 1.0 - soil%ssat ) * soil%css * soil%rhosoil & + soil%ssat * ( wblfsp * Ccswat * Cdensity_liq + & - ssnow%wbfice(:,k) * Ccsice * Cdensity_liq * 0.9 ) ) & + ssnow%wbfice(:,k) * Ccsice * Cdensity_ice ) ) & * soil%zse(k) dtg = dels / ssnow%gammzz(:,k) @@ -183,8 +183,7 @@ SUBROUTINE stempv(dels, canopy, ssnow, soil,heat_cap_lower_limit) ssnow%gammzz(:,k) = MAX( ( 1.0 - soil%ssat ) * soil%css * & soil%rhosoil + soil%ssat * ( wblfsp * Ccswat * & - Cdensity_liq + ssnow%wbfice(:,k) * Ccsice * Cdensity_liq * & - 0.9) , & + Cdensity_liq + ssnow%wbfice(:,k) * Ccsice * Cdensity_ice) , & heat_cap_lower_limit(:,k) ) * soil%zse(k) dtg = dels / ssnow%gammzz(:,k) diff --git a/src/science/soilsnow/cbl_surfbv.F90 b/src/science/soilsnow/cbl_surfbv.F90 index a63a32250..b6117bb70 100644 --- a/src/science/soilsnow/cbl_surfbv.F90 +++ b/src/science/soilsnow/cbl_surfbv.F90 @@ -7,14 +7,12 @@ MODULE surfbv_mod CONTAINS SUBROUTINE surfbv (dels, met, ssnow, soil, veg, canopy ) -! subrs -USE smoisturev_mod, ONLY: smoisturev +USE smoisturev_mod, ONLY: smoisturev + USE cable_common_module ! data -USE cable_surface_types_mod, ONLY: lakes_cable -USE cable_common_module - -USE cable_common_module +USE cable_surface_types_mod, ONLY: lakes_cable +USE cable_other_constants_mod, ONLY : wilt_limitfactor IMPLICIT NONE @@ -53,15 +51,15 @@ SUBROUTINE surfbv (dels, met, ssnow, soil, veg, canopy ) DO k = 1, ms xxx = REAL( soil%ssat,r_2 ) ssnow%rnof1 = ssnow%rnof1 + REAL( MAX( ssnow%wb(:,k) - xxx, 0.0_r_2 ) & - * 1000.0 ) * soil%zse(k) - ssnow%wb(:,k) = MAX( 0.01_r_2, MIN( ssnow%wb(:,k), xxx ) ) + * Cdensity_liq ) * soil%zse(k) + ssnow%wb(:,k) = MAX( soil%swilt(:)/(2.*wilt_limitfactor), MIN( ssnow%wb(:,k), xxx ) ) END DO ! for deep runoff use wb-sfc, but this value not to exceed .99*wb-wbice ! account for soil/ice cracking ! fracm = MIN(0.2, 1. - MIN(1., ssnow%wb(:,ms) / soil%sfc ) ) ! ssnow%wb(:,ms) = ssnow%wb(:,ms) & - ! + fracm*ssnow%rnof1/(1000.0*soil%zse(ms)) + ! + fracm*ssnow%rnof1/(Cdensity_liq*soil%zse(ms)) ! ssnow%rnof1 = (1. - fracm) * ssnow%rnof1 ! Scaling runoff to kg/m^2/s to match rest of the model @@ -114,13 +112,13 @@ SUBROUTINE surfbv (dels, met, ssnow, soil, veg, canopy ) ssnow%sinfil = MIN( ssnow%rnof2, ssnow%wb_lake ) ! water that can be extracted from the rnof2 ssnow%rnof2 = MAX( 0.0, ssnow%rnof2 - ssnow%sinfil ) ssnow%wb_lake = MAX( 0.0, ssnow%wb_lake - ssnow%sinfil) - xxx = MAX(0.0_r_2, (ssnow%wb(:,ms) - REAL(soil%sfc(:),r_2))*soil%zse(ms)*1000.0) + xxx = MAX(0.0_r_2, (ssnow%wb(:,ms) - REAL(soil%sfc(:),r_2))*soil%zse(ms)*Cdensity_liq) ssnow%sinfil = MIN( REAL(xxx), ssnow%wb_lake ) - ssnow%wb(:,ms) = ssnow%wb(:,ms) - REAL(ssnow%sinfil / (soil%zse(ms)*1000.0), r_2) + ssnow%wb(:,ms) = ssnow%wb(:,ms) - REAL(ssnow%sinfil / (soil%zse(ms)*Cdensity_liq), r_2) ssnow%wb_lake = MAX( 0.0, ssnow%wb_lake - ssnow%sinfil) - xxx = MAX(0.0_r_2, (ssnow%wb(:,ms) - 0.5*(soil%sfc + soil%swilt))*soil%zse(ms)*1000.0) + xxx = MAX(0.0_r_2, (ssnow%wb(:,ms) - 0.5*(soil%sfc + soil%swilt))*soil%zse(ms)*Cdensity_liq) ssnow%sinfil = MIN( REAL(xxx), ssnow%wb_lake ) - ssnow%wb(:,ms) = ssnow%wb(:,ms) - ssnow%sinfil / (soil%zse(ms)*1000.0) + ssnow%wb(:,ms) = ssnow%wb(:,ms) - ssnow%sinfil / (soil%zse(ms)*Cdensity_liq) ssnow%wb_lake = MAX( 0.0, ssnow%wb_lake - ssnow%sinfil) ENDWHERE diff --git a/src/science/soilsnow/cbl_thermal.F90 b/src/science/soilsnow/cbl_thermal.F90 index cb682df1f..163b9266c 100644 --- a/src/science/soilsnow/cbl_thermal.F90 +++ b/src/science/soilsnow/cbl_thermal.F90 @@ -6,11 +6,7 @@ MODULE snow_processes_soil_thermal_mod CONTAINS -!SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal,snowmlt) ! replaced by rk4417 - phase2 -! snowmlt is not used in cable_gw_hydro_module and as far as I can tell serves no purpose in the code - rk4417 - -SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal) -!* calculate snow processes and thermal soil +SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal,snowmlt) USE snowl_adjust_mod, ONLY: snowl_adjust USE snowCheck_mod, ONLY: snowCheck @@ -28,12 +24,10 @@ SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal) TYPE(veg_parameter_type), INTENT(INOUT) :: veg TYPE(met_type), INTENT(INOUT) :: met ! all met forcing TYPE (balances_type), INTENT(INOUT) :: bal - REAL, DIMENSION(mp) :: snowmlt !track snow melt -! REAL, DIMENSION(:), INTENT(INOUT) :: snowmlt ! replaced by rk4417 - phase2 + REAL, DIMENSION(:), INTENT(INOUT) :: snowmlt + INTEGER :: k,i - snowmlt = 0.0 ! inserted by rk4417 - phase2 - CALL snowcheck (dels, ssnow, soil, met ) CALL snowdensity (dels, ssnow, soil) From e0c0617d133c2ff65c6f786bdf36e244978731c7 Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Wed, 27 Nov 2024 18:12:40 +1100 Subject: [PATCH 02/16] rollback --- src/science/soilsnow/cbl_thermal.F90 | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/science/soilsnow/cbl_thermal.F90 b/src/science/soilsnow/cbl_thermal.F90 index 163b9266c..cb682df1f 100644 --- a/src/science/soilsnow/cbl_thermal.F90 +++ b/src/science/soilsnow/cbl_thermal.F90 @@ -6,7 +6,11 @@ MODULE snow_processes_soil_thermal_mod CONTAINS -SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal,snowmlt) +!SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal,snowmlt) ! replaced by rk4417 - phase2 +! snowmlt is not used in cable_gw_hydro_module and as far as I can tell serves no purpose in the code - rk4417 + +SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal) +!* calculate snow processes and thermal soil USE snowl_adjust_mod, ONLY: snowl_adjust USE snowCheck_mod, ONLY: snowCheck @@ -24,10 +28,12 @@ SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal,snowml TYPE(veg_parameter_type), INTENT(INOUT) :: veg TYPE(met_type), INTENT(INOUT) :: met ! all met forcing TYPE (balances_type), INTENT(INOUT) :: bal - REAL, DIMENSION(:), INTENT(INOUT) :: snowmlt - + REAL, DIMENSION(mp) :: snowmlt !track snow melt +! REAL, DIMENSION(:), INTENT(INOUT) :: snowmlt ! replaced by rk4417 - phase2 INTEGER :: k,i + snowmlt = 0.0 ! inserted by rk4417 - phase2 + CALL snowcheck (dels, ssnow, soil, met ) CALL snowdensity (dels, ssnow, soil) From 62361981aa60c99772616e8dde2bd2e33ae8d6be Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Fri, 13 Dec 2024 12:42:05 +1100 Subject: [PATCH 03/16] remove reduntant file --- .../util/init/cbl_um_init_veg.F90.manual | 77 ------------------- 1 file changed, 77 deletions(-) delete mode 100644 src/coupled/AM3/control/cable/util/init/cbl_um_init_veg.F90.manual diff --git a/src/coupled/AM3/control/cable/util/init/cbl_um_init_veg.F90.manual b/src/coupled/AM3/control/cable/util/init/cbl_um_init_veg.F90.manual deleted file mode 100644 index b1548b0d8..000000000 --- a/src/coupled/AM3/control/cable/util/init/cbl_um_init_veg.F90.manual +++ /dev/null @@ -1,77 +0,0 @@ -MODULE cbl_um_init_veg_mod - -IMPLICIT NONE - -PUBLIC initialize_veg - -CONTAINS - -SUBROUTINE initialize_veg( SurfaceType, SoilType, mp,& - ms, nrb, npft, nsurft, land_pts, l_tile_pts, & - ICE_SurfaceType, ICE_SoilType, LAI_thresh, & - soil_zse, surft_pts, surft_index, tile_frac, & - veg, soil, pars ) - -USE init_cable_parms_mod, ONLY: init_cable_parms_expl -USE params_io_mod_cbl, ONLY: params_io_data_type -USE cable_def_types_mod, ONLY: veg_parameter_type -USE cable_def_types_mod, ONLY: soil_parameter_type, r_2 - -IMPLICIT NONE - -INTEGER, INTENT(IN) :: mp -INTEGER, INTENT(IN) :: ms -INTEGER, INTENT(IN) :: nrb -INTEGER, INTENT(IN) :: nsurft -INTEGER, INTENT(IN) :: npft -INTEGER, INTENT(IN) :: land_pts -INTEGER, INTENT(IN) :: ICE_SurfaceType ! index ICE surface type -INTEGER, INTENT(IN) :: ICE_SoilType ! index soil type -REAL, INTENT(IN) :: lai_thresh -REAL, INTENT(IN) :: soil_zse(ms) ! soil depth per layer -INTEGER, INTENT(IN) :: surft_pts(nsurft) ! # land points per tile -INTEGER, INTENT(IN) :: surft_index(land_pts, nsurft) ! land_pt index of point -LOGICAL, INTENT(IN) :: L_tile_pts(land_pts, nsurft) ! TRUE if active tile -REAL, INTENT(IN) :: tile_frac(land_pts,nsurft) - -INTEGER, INTENT(OUT) :: SurfaceType(mp) ! surface tile PFT/nveg -INTEGER, INTENT(OUT) :: SoilType(mp) ! soil type per tile - -TYPE(veg_parameter_type), INTENT(OUT) :: veg ! vegetation parameters -TYPE(soil_parameter_type), INTENT(OUT) :: soil ! soil parameters -TYPE(params_io_data_type), INTENT(IN) :: pars - -! local vars -INTEGER :: i -REAL :: totdepth - -CALL init_cable_parms_expl( SurfaceType, SoilType, mp, ms, nrb, land_pts, & - nsurft, l_tile_pts, ICE_SurfaceType, & - ICE_SoilType, soil_zse, veg, soil, pars, tile_frac ) - -veg%meth = 1 ! review !should be namelist config if even used -veg%ejmax = 2.*veg%vcmax -veg%gamma = 3.e-2 ! for Haverd2013 switch ! offline set init _parameters -veg%F10 = 0.85 ! offline set init _parameters - -! calculate veg%froot from using rootbeta and soil depth -! (Jackson et al. 1996, Oceologica, 108:389-411) -totdepth = 0.0 -DO i = 1, ms-1 - totdepth = totdepth + soil_zse(i) * 100.0 ! unit in centimetres - veg%froot(:, i) = MIN( 1.0_r_2, 1.0-veg%rootbeta(:)**totdepth ) -END DO -veg%froot(:, ms) = 1.0 - veg%froot(:, ms-1) -DO i = ms-1, 2, -1 - veg%froot(:, i) = veg%froot(:, i)-veg%froot(:,i-1) -END DO - - -RETURN -END SUBROUTINE initialize_veg - -END MODULE cbl_um_init_veg_mod - - - - From 5983ee935b4cfd5f8a9cf354224ae9c2e9e79b63 Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Fri, 27 Dec 2024 22:38:41 +1100 Subject: [PATCH 04/16] how did this happen - droped a whole calc --- src/science/soilsnow/cbl_snowAccum.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/science/soilsnow/cbl_snowAccum.F90 b/src/science/soilsnow/cbl_snowAccum.F90 index 55af4ac17..27ed28d03 100644 --- a/src/science/soilsnow/cbl_snowAccum.F90 +++ b/src/science/soilsnow/cbl_snowAccum.F90 @@ -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 & + / ( 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 & @@ -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) From 319dc6a1eb0e40d04a1d0486da1c15b4d38ba3c6 Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Mon, 13 Jan 2025 09:07:56 +1100 Subject: [PATCH 05/16] run with @Whyborn's revised grid_info completes IF this error check is switched off --- src/offline/cable_mpimaster.F90 | 34 ++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/offline/cable_mpimaster.F90 b/src/offline/cable_mpimaster.F90 index 83ea7d813..1d3054c1e 100644 --- a/src/offline/cable_mpimaster.F90 +++ b/src/offline/cable_mpimaster.F90 @@ -796,25 +796,25 @@ SUBROUTINE mpidrv_master (comm, trunk_sumbal, dels, koffset, kend, PLUME, CRU) !$ if (ktau == kend-1) PRINT*, "sum_fe[Wm-2], sum_fpn[umol/m2/s]", & !$ new_sumfe/count_bal, new_sumfpn/count_bal - ! check for Nans in biophysical outputs and abort if there are any - IF (ANY( canopy%fe.NE. canopy%fe)) THEN - DO kk=1,mp - - IF (canopy%fe(kk).NE. canopy%fe(kk)) THEN - WRITE(*,*) 'Nan in evap flux,', kk, patch(kk)%latitude, patch(kk)%longitude - WRITE(*,*) 'fe nan', kk, ktau,met%qv(kk), met%precip(kk),met%precip_sn(kk), & - met%fld(kk), met%fsd(kk,:), met%tk(kk), met%ua(kk), & - ssnow%potev(kk), met%pmb(kk), & - canopy%ga(kk), ssnow%tgg(kk,:), canopy%fwsoil(kk), & - rad%fvlai(kk,:) , rad%fvlai(kk,1), & - rad%fvlai(kk,2), canopy%vlaiw(kk) - - CALL MPI_Abort(comm, 0, ierr) - ENDIF + !jhan:test!! check for Nans in biophysical outputs and abort if there are any + !jhan:test!IF (ANY( canopy%fe.NE. canopy%fe)) THEN + !jhan:test! DO kk=1,mp - ENDDO + !jhan:test! IF (canopy%fe(kk).NE. canopy%fe(kk)) THEN + !jhan:test! WRITE(*,*) 'Nan in evap flux,', kk, patch(kk)%latitude, patch(kk)%longitude + !jhan:test! WRITE(*,*) 'fe nan', kk, ktau,met%qv(kk), met%precip(kk),met%precip_sn(kk), & + !jhan:test! met%fld(kk), met%fsd(kk,:), met%tk(kk), met%ua(kk), & + !jhan:test! ssnow%potev(kk), met%pmb(kk), & + !jhan:test! canopy%ga(kk), ssnow%tgg(kk,:), canopy%fwsoil(kk), & + !jhan:test! rad%fvlai(kk,:) , rad%fvlai(kk,1), & + !jhan:test! rad%fvlai(kk,2), canopy%vlaiw(kk) - ENDIF + !jhan:test! CALL MPI_Abort(comm, 0, ierr) + !jhan:test! ENDIF + + !jhan:test! ENDDO + + !jhan:test!ENDIF IF(ktau==(kend-1)) THEN nkend = nkend+1 From 0e8496236666162d872214238124e1752c58983b Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Tue, 14 Jan 2025 13:12:43 +1100 Subject: [PATCH 06/16] clobber ice point soil parameters --- src/offline/cable_parameters.F90 | 62 ++++++++++++++++++++++---------- 1 file changed, 44 insertions(+), 18 deletions(-) diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index 0c553c89b..cb1ffd661 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -736,24 +736,39 @@ SUBROUTINE spatialSoil(nlon, nlat, logn) ok = NF90_CLOSE(ncid) IF (ok /= NF90_NOERR) CALL nc_abort(ok, 'Error closing IGBP soil map.') - ! Code if using UM soil file - ! unit change and glacial-point check were done in preprocessing - ! ! change unit to m/s - ! inhyds = inhyds * 1.0E-3 - ! ! Assign values to glacial points which are zeroes - ! WHERE(inswilt==0.) inswilt = 0.216 - ! WHERE( insfc ==0.) insfc = 0.301 - ! WHERE(inssat ==0.) inssat = 0.479 - ! WHERE( inbch ==0.) inbch = 7.1 - ! WHERE(inhyds ==0.) inhyds = 1.E-3 - ! WHERE(insucs ==0.) insucs = 0.153 - ! WHERE(inrhosoil==0.) inrhosoil = 1455 - ! WHERE(incnsd ==0.) incnsd = 0.272 - ! WHERE(incss > 630000.0) - ! incss = incss / inrhosoil ! normal points need unit conversion - ! ELSEWHERE - ! incss = 2100.0 ! glacial points - ! ENDWHERE +!ice values from nml file +!!soilin%bch = 7.1 +!!soilin%clay = 0.3 +!!soilin%css = 2100 +!!soilin%hyds = 0.000001 +!!soilin%rhosoil = 917 +!!soilin%sand = 0.37 +!!soilin%sfc = 0.301 +!!soilin%silt = 0.33 +!!soilin%ssat = 0.479 +!!soilin%sucs = -0.153 +!!soilin%swilt = 0.216 + + !IDK where this code came from. Some are == nml values (for ice) + !units change might be reasonable but others !!!!! e.g. density=1455?? + !Code if using UM soil file + !unit change and glacial-point check were done in preprocessing + ! change unit to m/s + inhyds = inhyds * 1.0E-3 + ! Assign values to glacial points which are zeroes + WHERE(inswilt==0.) inswilt = 0.216 + WHERE( insfc ==0.) insfc = 0.301 + WHERE(inssat ==0.) inssat = 0.479 + WHERE( inbch ==0.) inbch = 7.1 + WHERE(inhyds ==0.) inhyds = 1.E-3 + WHERE(insucs ==0.) insucs = 0.153 + WHERE(inrhosoil==0.) inrhosoil = 1455 + WHERE(incnsd ==0.) incnsd = 0.272 + WHERE(incss > 630000.0) + incss = incss / inrhosoil ! normal points need unit conversion + ELSEWHERE + incss = 2100.0 ! glacial points + ENDWHERE ! Calculate albedo for radiation bands and overwrite previous ! initialization @@ -1656,6 +1671,17 @@ SUBROUTINE write_default_params(met, air, ssnow, veg, bgc, & ! SLI specific initialisations: ! IF(cable_user%SOIL_STRUC=='sli') THEN ssnow%h0(:) = 0.0 +!! DO i=1,mp +!! !IF(i == 1852) THEN +!! IF( soil%ssat(i) == 0.0 ) THEN +!! write(*,*) "jh:soil%isoilm ", soil%isoilm(i) +!! !write(*,*) "jh:ssnow%S ", ssnow%S(i,:) +!! !write(*,*) "jh:ssnow%wb ", ssnow%wb(i,:) +!! !write(*,*) "jh:soil%ssat ", soil%ssat(i) +!! !write(*,*) "jh:SPREAD ", SPREAD(soil%ssat,2,ms) +!! END IF +!! END DO +!!!stop ssnow%S(:,:) = ssnow%wb(:,:)/SPREAD(soil%ssat,2,ms) ssnow%snowliq(:,:) = 0.0 ssnow%Tsurface = 25.0 From a3ba35b50ff372c6f14b1b7bfe09642a37178165 Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Wed, 27 Nov 2024 17:22:40 +1100 Subject: [PATCH 07/16] developments from AM3 --- src/science/soilsnow/cbl_thermal.F90 | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/src/science/soilsnow/cbl_thermal.F90 b/src/science/soilsnow/cbl_thermal.F90 index cb682df1f..163b9266c 100644 --- a/src/science/soilsnow/cbl_thermal.F90 +++ b/src/science/soilsnow/cbl_thermal.F90 @@ -6,11 +6,7 @@ MODULE snow_processes_soil_thermal_mod CONTAINS -!SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal,snowmlt) ! replaced by rk4417 - phase2 -! snowmlt is not used in cable_gw_hydro_module and as far as I can tell serves no purpose in the code - rk4417 - -SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal) -!* calculate snow processes and thermal soil +SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal,snowmlt) USE snowl_adjust_mod, ONLY: snowl_adjust USE snowCheck_mod, ONLY: snowCheck @@ -28,12 +24,10 @@ SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal) TYPE(veg_parameter_type), INTENT(INOUT) :: veg TYPE(met_type), INTENT(INOUT) :: met ! all met forcing TYPE (balances_type), INTENT(INOUT) :: bal - REAL, DIMENSION(mp) :: snowmlt !track snow melt -! REAL, DIMENSION(:), INTENT(INOUT) :: snowmlt ! replaced by rk4417 - phase2 + REAL, DIMENSION(:), INTENT(INOUT) :: snowmlt + INTEGER :: k,i - snowmlt = 0.0 ! inserted by rk4417 - phase2 - CALL snowcheck (dels, ssnow, soil, met ) CALL snowdensity (dels, ssnow, soil) From be834b8b90d6a1ea34d9f91c8dbe6080ea2adfa9 Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Wed, 27 Nov 2024 18:12:40 +1100 Subject: [PATCH 08/16] rollback --- src/science/soilsnow/cbl_thermal.F90 | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/science/soilsnow/cbl_thermal.F90 b/src/science/soilsnow/cbl_thermal.F90 index 163b9266c..cb682df1f 100644 --- a/src/science/soilsnow/cbl_thermal.F90 +++ b/src/science/soilsnow/cbl_thermal.F90 @@ -6,7 +6,11 @@ MODULE snow_processes_soil_thermal_mod CONTAINS -SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal,snowmlt) +!SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal,snowmlt) ! replaced by rk4417 - phase2 +! snowmlt is not used in cable_gw_hydro_module and as far as I can tell serves no purpose in the code - rk4417 + +SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal) +!* calculate snow processes and thermal soil USE snowl_adjust_mod, ONLY: snowl_adjust USE snowCheck_mod, ONLY: snowCheck @@ -24,10 +28,12 @@ SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal,snowml TYPE(veg_parameter_type), INTENT(INOUT) :: veg TYPE(met_type), INTENT(INOUT) :: met ! all met forcing TYPE (balances_type), INTENT(INOUT) :: bal - REAL, DIMENSION(:), INTENT(INOUT) :: snowmlt - + REAL, DIMENSION(mp) :: snowmlt !track snow melt +! REAL, DIMENSION(:), INTENT(INOUT) :: snowmlt ! replaced by rk4417 - phase2 INTEGER :: k,i + snowmlt = 0.0 ! inserted by rk4417 - phase2 + CALL snowcheck (dels, ssnow, soil, met ) CALL snowdensity (dels, ssnow, soil) From da2aef5673a59ab0ecc85a0e87c4ca72d8a79351 Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Tue, 14 Jan 2025 22:33:35 +1100 Subject: [PATCH 09/16] changes to accomodate new grid info --- src/offline/cable_parameters.F90 | 56 ++++++++++++++++++++++---------- 1 file changed, 38 insertions(+), 18 deletions(-) diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index 0c553c89b..b4bccfae5 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -736,24 +736,44 @@ SUBROUTINE spatialSoil(nlon, nlat, logn) ok = NF90_CLOSE(ncid) IF (ok /= NF90_NOERR) CALL nc_abort(ok, 'Error closing IGBP soil map.') - ! Code if using UM soil file - ! unit change and glacial-point check were done in preprocessing - ! ! change unit to m/s - ! inhyds = inhyds * 1.0E-3 - ! ! Assign values to glacial points which are zeroes - ! WHERE(inswilt==0.) inswilt = 0.216 - ! WHERE( insfc ==0.) insfc = 0.301 - ! WHERE(inssat ==0.) inssat = 0.479 - ! WHERE( inbch ==0.) inbch = 7.1 - ! WHERE(inhyds ==0.) inhyds = 1.E-3 - ! WHERE(insucs ==0.) insucs = 0.153 - ! WHERE(inrhosoil==0.) inrhosoil = 1455 - ! WHERE(incnsd ==0.) incnsd = 0.272 - ! WHERE(incss > 630000.0) - ! incss = incss / inrhosoil ! normal points need unit conversion - ! ELSEWHERE - ! incss = 2100.0 ! glacial points - ! ENDWHERE +!ice values from nml file +!!dont even esxist in what we have +!!soilin%clay = 0.3 +!!soilin%silt = 0.33 + +!!agreement with what we have +!!soilin%css = 2100 +!!soilin%hyds = 0.000001 (*1e3?) +!!soilin%rhosoil = 917 (correct ice density) +!!soilin%sucs = -0.153 (sign) + +!!agreement with what we have +!!soilin%bch = 7.1 +!!soilin%sand = 0.37 +!!soilin%sfc = 0.301 +!!soilin%ssat = 0.479 +!!soilin%swilt = 0.216 + +!IDK where this code came from. Some are == nml values (for ice) +!units change might be reasonable but others !!!!! e.g. density=1455?? +!Code if using UM soil file +!unit change and glacial-point check were done in preprocessing +! change unit to m/s +! Assign values to glacial points which are zeroes +WHERE( inswilt == 0.0 ) inswilt = 0.216 +WHERE( insfc == 0.0 ) insfc = 0.301 +WHERE( inssat == 0.0 ) inssat = 0.479 +WHERE( inbch == 0.0 ) inbch = 7.1 +WHERE( inhyds == 0.0 ) inhyds = 1.E-3 +WHERE( insucs == 0.0 ) insucs = -0.153 +WHERE( inrhosoil == 0.0 ) inrhosoil = 917 +WHERE( incnsd == 0.0 ) incnsd = 0.272 + +WHERE(incss > 630000.0) + incss = incss / inrhosoil ! normal points need unit conversion +ELSEWHERE + incss = 2100.0 ! glacial points +ENDWHERE ! Calculate albedo for radiation bands and overwrite previous ! initialization From 1e508a5ea460bacb8dddac2b8933bc78ba210a30 Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Wed, 15 Jan 2025 07:55:23 +1100 Subject: [PATCH 10/16] update commments --- src/offline/cable_parameters.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index b4bccfae5..5a6a80e70 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -736,13 +736,13 @@ SUBROUTINE spatialSoil(nlon, nlat, logn) ok = NF90_CLOSE(ncid) IF (ok /= NF90_NOERR) CALL nc_abort(ok, 'Error closing IGBP soil map.') -!ice values from nml file -!!dont even esxist in what we have +!comparing mysterious existing loop with ice values from nml file +!!dont even exist in what we have !!soilin%clay = 0.3 !!soilin%silt = 0.33 -!!agreement with what we have -!!soilin%css = 2100 +!!quasi-agreement with what we have. Below have adpopted these values +!!soilin%css = 2100 (see WHERE() below ) !!soilin%hyds = 0.000001 (*1e3?) !!soilin%rhosoil = 917 (correct ice density) !!soilin%sucs = -0.153 (sign) From 535a3342ba5d6c8b8fb118a9b101d7dd25c8ac52 Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Thu, 16 Jan 2025 16:09:37 +1100 Subject: [PATCH 11/16] rollback to main version --- src/offline/cable_parameters.F90 | 56 ++++++++++---------------------- 1 file changed, 18 insertions(+), 38 deletions(-) diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index 5a6a80e70..0c553c89b 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -736,44 +736,24 @@ SUBROUTINE spatialSoil(nlon, nlat, logn) ok = NF90_CLOSE(ncid) IF (ok /= NF90_NOERR) CALL nc_abort(ok, 'Error closing IGBP soil map.') -!comparing mysterious existing loop with ice values from nml file -!!dont even exist in what we have -!!soilin%clay = 0.3 -!!soilin%silt = 0.33 - -!!quasi-agreement with what we have. Below have adpopted these values -!!soilin%css = 2100 (see WHERE() below ) -!!soilin%hyds = 0.000001 (*1e3?) -!!soilin%rhosoil = 917 (correct ice density) -!!soilin%sucs = -0.153 (sign) - -!!agreement with what we have -!!soilin%bch = 7.1 -!!soilin%sand = 0.37 -!!soilin%sfc = 0.301 -!!soilin%ssat = 0.479 -!!soilin%swilt = 0.216 - -!IDK where this code came from. Some are == nml values (for ice) -!units change might be reasonable but others !!!!! e.g. density=1455?? -!Code if using UM soil file -!unit change and glacial-point check were done in preprocessing -! change unit to m/s -! Assign values to glacial points which are zeroes -WHERE( inswilt == 0.0 ) inswilt = 0.216 -WHERE( insfc == 0.0 ) insfc = 0.301 -WHERE( inssat == 0.0 ) inssat = 0.479 -WHERE( inbch == 0.0 ) inbch = 7.1 -WHERE( inhyds == 0.0 ) inhyds = 1.E-3 -WHERE( insucs == 0.0 ) insucs = -0.153 -WHERE( inrhosoil == 0.0 ) inrhosoil = 917 -WHERE( incnsd == 0.0 ) incnsd = 0.272 - -WHERE(incss > 630000.0) - incss = incss / inrhosoil ! normal points need unit conversion -ELSEWHERE - incss = 2100.0 ! glacial points -ENDWHERE + ! Code if using UM soil file + ! unit change and glacial-point check were done in preprocessing + ! ! change unit to m/s + ! inhyds = inhyds * 1.0E-3 + ! ! Assign values to glacial points which are zeroes + ! WHERE(inswilt==0.) inswilt = 0.216 + ! WHERE( insfc ==0.) insfc = 0.301 + ! WHERE(inssat ==0.) inssat = 0.479 + ! WHERE( inbch ==0.) inbch = 7.1 + ! WHERE(inhyds ==0.) inhyds = 1.E-3 + ! WHERE(insucs ==0.) insucs = 0.153 + ! WHERE(inrhosoil==0.) inrhosoil = 1455 + ! WHERE(incnsd ==0.) incnsd = 0.272 + ! WHERE(incss > 630000.0) + ! incss = incss / inrhosoil ! normal points need unit conversion + ! ELSEWHERE + ! incss = 2100.0 ! glacial points + ! ENDWHERE ! Calculate albedo for radiation bands and overwrite previous ! initialization From a2392909042b071ed91ebb360c5cda9de5c013dc Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Thu, 16 Jan 2025 17:50:31 +1100 Subject: [PATCH 12/16] rm legacy type casting --- src/science/soilsnow/cbl_soilsnow_main.F90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/science/soilsnow/cbl_soilsnow_main.F90 b/src/science/soilsnow/cbl_soilsnow_main.F90 index 1ed8af8ae..c7fd4d7db 100644 --- a/src/science/soilsnow/cbl_soilsnow_main.F90 +++ b/src/science/soilsnow/cbl_soilsnow_main.F90 @@ -122,7 +122,7 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) ! snow aging etc... CALL snowl_adjust(dels, ssnow, canopy ) - CALL stempv(dels, canopy, ssnow, soil, REAL(soil%heat_cap_lower_limit) ) + CALL stempv(dels, canopy, ssnow, soil, soil%heat_cap_lower_limit ) ssnow%tss = (1-ssnow%isflag)*ssnow%tgg(:,1) + ssnow%isflag*ssnow%tggsn(:,1) @@ -133,8 +133,7 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) CALL remove_trans(dels, soil, ssnow, canopy, veg) - CALL soilfreeze(dels, soil, ssnow, REAL(soil%heat_cap_lower_limit) ) - + CALL soilfreeze(dels, soil, ssnow, soil%heat_cap_lower_limit ) totwet = canopy%precis + ssnow%smelt From 07204418fb146149b758ed00de04042c165b67e1 Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Fri, 17 Jan 2025 15:26:04 +1100 Subject: [PATCH 13/16] rm redundant file --- src/params/cable_params_mod.F90 | 255 -------------------------------- 1 file changed, 255 deletions(-) delete mode 100644 src/params/cable_params_mod.F90 diff --git a/src/params/cable_params_mod.F90 b/src/params/cable_params_mod.F90 deleted file mode 100644 index 8bd5f73ec..000000000 --- a/src/params/cable_params_mod.F90 +++ /dev/null @@ -1,255 +0,0 @@ -MODULE cable_params_mod - -!H! Elevate these to namelist definable -USE cable_other_constants_mod, ONLY: nsl -USE cable_other_constants_mod, ONLY: nrb ! # radiation "bANDS" - !dir/dif components in bands VIS/NIR -USE cable_other_constants_mod, ONLY: nscs ! number of soil carbon stores -USE cable_other_constants_mod, ONLY: nvcs ! number of vegetation carbon stores -USE grid_constants_mod_cbl, ONLY : mstype => nsoil_max ! # of soil types [9] -USE grid_constants_mod_cbl, ONLY : ntype_max ! # of PFTs [17] - -IMPLICIT NONE - -PUBLIC :: veg_parameter_type -PUBLIC :: vegin_type -PUBLIC :: soil_parameter_type -PUBLIC :: soilin_type - -PUBLIC :: vegin -PUBLIC :: soilin - -!----------------------------------------------------------------------------- -! Description: -! Defines variable types and variables for CABLE standalone runs. -! Based on cable_def_types_mod.F90 from the CABLE trunk. -! -! Code Owner: Please refer to ModuleLeaders.txt -! This file belongs in CABLE SCIENCE -!----------------------------------------------------------------------------- - -! Vegetation parameters used: -TYPE veg_parameter_type -!jhan: surely %meth doesn't need to be 'mp' here - INTEGER, DIMENSION(:), POINTER :: & - meth, & ! method for calculation of canopy fluxes and temp. - iveg , & ! vegetation type - iLU ! land use type - REAL, DIMENSION(:), POINTER :: & - canst1, & ! max intercepted water by canopy (mm/LAI) - dleaf, & ! chararacteristc legnth of leaf (m) - ejmax, & ! max pot. electron transp rate top leaf(mol/m2/s) - frac4, & ! fraction of c4 plants - hc, & ! roughness height of canopy (veg - snow) - vlai, & ! leaf area index - xalbnir, & - rp20, & ! plant respiration coefficient at 20 C - rpcoef, & ! temperature coef nonleaf plant respiration (1/C) - rs20, & ! soil respiration at 20 C [mol m-2 s-1] - shelrb, & ! sheltering factor (dimensionless) - vegcf, & ! kdcorbin, 08/10 - tminvj, & ! min temperature of the start of photosynthesis - toptvj, & ! opt temperature of the start of photosynthesis - tmaxvj, & ! max temperature of the start of photosynthesis - vbeta, & ! - vcmax, & ! max RuBP carboxylation rate top leaf (mol/m2/s) - xfang, & ! leaf angle PARAMETER - extkn, & ! extinction coef for vertical - vlaimax, & ! extinction coef for vertical - wai, & ! wood area index (stem+branches+twigs) - a1gs, & ! a1 parameter in stomatal conductance model - d0gs, & ! d0 in stomatal conductance model - alpha, & ! initial slope of J-Q response curve - convex, & ! convexity of J-Q response curve - cfrd, & ! ratio of day respiration to vcmax - gswmin, & ! minimal stomatal conductance - conkc0, & ! Michaelis-menton constant for carboxylase - conko0, & ! Michaelis-menton constant for oxygenase - ekc, & ! activation energy for caroxylagse - eko, & ! acvtivation enegery for oxygenase - g0, & ! Belinda's stomatal model intercept, Ticket #56. - g1 ! Belinda's stomatal model slope, Ticket #56. - - LOGICAL, DIMENSION(:), POINTER :: & - deciduous ! flag used for phenology fix - - REAL, DIMENSION(:,:), POINTER :: & - refl, & - taul, & - froot ! fraction of root in each soil layer - - ! Additional veg parameters: - REAL, DIMENSION(:), POINTER :: rootbeta ! parameter for estimating vertical root mass distribution (froot) - REAL, DIMENSION(:), POINTER :: gamma ! parameter in root efficiency function (Lai and Katul 2000) - REAL, DIMENSION(:), POINTER :: ZR ! maximum rooting depth (cm) - REAL, DIMENSION(:), POINTER :: F10 ! fraction of roots in top 10 cm - - REAL, DIMENSION(:), POINTER :: clitt ! - - ! Additional POP veg param - INTEGER, DIMENSION(:,:), POINTER :: disturbance_interval - REAL, DIMENSION(:,:), POINTER :: disturbance_intensity - - REAL, POINTER :: & - soil(:,:), & - csoil(:,:), & - cplant(:,:), & - ratecs(:,:), & - ratecp(:,:) - - END TYPE veg_parameter_type - - -! Vegetation parameters I/O: -TYPE vegin_type - REAL :: & - canst1(ntype_max), & - dleaf(ntype_max), & - length(ntype_max), & - width(ntype_max), & - vcmax(ntype_max), & - ejmax(ntype_max), & - hc(ntype_max), & - xfang(ntype_max), & - rp20(ntype_max), & - rpcoef(ntype_max), & - rs20(ntype_max), & - wai(ntype_max), & - rootbeta(ntype_max), & - shelrb(ntype_max), & - vegcf(ntype_max), & - frac4(ntype_max), & - xalbnir(ntype_max), & - extkn(ntype_max), & - tminvj(ntype_max), & - tmaxvj(ntype_max), & - vbeta(ntype_max), & - a1gs(ntype_max), & - d0gs(ntype_max), & - alpha(ntype_max), & - convex(ntype_max), & - cfrd(ntype_max), & - gswmin(ntype_max), & - conkc0(ntype_max), & - conko0(ntype_max), & - ekc(ntype_max), & - eko(ntype_max), & - g0(ntype_max), & - g1(ntype_max), & - zr(ntype_max), & - clitt(ntype_max), & - froot(nsl,ntype_max), & - csoil(nscs,ntype_max), & - ratecs(nscs,ntype_max), & - cplant(nvcs,ntype_max), & - ratecp(nvcs,ntype_max), & - refl(nrb,ntype_max), & - taul(nrb,ntype_max) -END TYPE vegin_type - -! Soil parameters used: - TYPE soil_parameter_type - - INTEGER, DIMENSION(:), POINTER :: & - isoilm ! integer soil type - - REAL, DIMENSION(:), POINTER :: & - bch, & ! parameter b in Campbell equation - c3, & ! c3 drainage coeff (fraction) - clay, & ! fraction of soil which is clay - css, & ! soil specific heat capacity [kJ/kg/K] - hsbh, & ! difsat * etasat (=hyds*abs(sucs)*bch) - hyds, & ! hydraulic conductivity @ saturation [m/s], Ksat - i2bp3, & ! par. one in K vis suction (=nint(bch)+2) - ibp2, & ! par. two in K vis suction (fn of pbch) - rhosoil, & ! soil density [kg/m3] - sand, & ! fraction of soil which is sand - sfc, & ! vol H2O @ field capacity - silt, & ! fraction of soil which is silt - ssat, & ! vol H2O @ saturation - sucs, & ! suction at saturation (m) - swilt, & ! vol H2O @ wilting - zse, & ! thickness of each soil layer (1=top) in m - zshh, & ! distance between consecutive layer midpoints (m) - ! vars intro for Ticket #27 - soilcol, & ! keep color for all patches/tiles - albsoilf ! soil reflectance - - REAL, DIMENSION(:,:), POINTER :: & - heat_cap_lower_limit - - REAL, DIMENSION(:,:), POINTER :: & - zse_vec,css_vec,cnsd_vec - - REAL, DIMENSION(:), POINTER :: & - cnsd, & ! thermal conductivity of dry soil [W/m/K] - pwb_min ! working variable (swilt/ssat)**ibp2 - - REAL, DIMENSION(:,:), POINTER :: & - albsoil ! soil reflectance (2nd dim. BP 21Oct2009) - !mrd561 - !MD parameters for GW module that vary with soil layer - REAL, DIMENSION(:,:), POINTER :: & - sucs_vec, & !psi at saturation in [mm] - hyds_vec, & !saturated hydraulic conductivity [mm/s] - bch_vec, & !C and H B [none] - clay_vec, & !fraction of soil that is clay [frac] - sand_vec, & !fraction of soil that is sand [frac] - silt_vec, & !fraction of soil that is silt [frac] - org_vec, & !fration of soil made of organic soils [frac] - rhosoil_vec,& !soil density [kg/m3] - ssat_vec, & !volumetric water content at saturation [mm3/mm3] - watr, & !residual water content of the soil [mm3/mm3] - sfc_vec, & !field capcacity (hk = 1 mm/day) - swilt_vec ! wilting point (hk = 0.02 mm/day) - - REAL, DIMENSION(:), POINTER :: & - drain_dens,&! drainage density ( mean dist to rivers/streams ) - elev, & !elevation above sea level - elev_std, & !elevation above sea level - slope, & !mean slope of grid cell - slope_std !stddev of grid cell slope - - !MD parameters for GW module for the aquifer - REAL, DIMENSION(:), POINTER :: & - GWsucs_vec, & !head in the aquifer [mm] - GWhyds_vec, & !saturated hydraulic conductivity of the aquifer [mm/s] - GWbch_vec, & !clapp and horn b of the aquifer [none] - GWssat_vec, & !saturated water content of the aquifer [mm3/mm3] - GWwatr, & !residual water content of the aquifer [mm3/mm3] - GWz, & !node depth of the aquifer [m] - GWdz, & !thickness of the aquifer [m] - GWrhosoil_vec !density of the aquifer substrate [kg/m3] - - ! Additional SLI parameters - INTEGER, DIMENSION(:), POINTER :: nhorizons ! number of soil horizons - INTEGER, DIMENSION(:,:), POINTER :: ishorizon ! horizon number 1:nhorizons - REAL, DIMENSION(:), POINTER :: clitt ! litter (tC/ha) - REAL, DIMENSION(:), POINTER :: zeta ! macropore parameter - REAL, DIMENSION(:), POINTER :: fsatmax ! variably saturated area parameter - - END TYPE soil_parameter_type - - - -! Soil parameters I/O: -TYPE soilin_type - REAL :: & - silt(mstype), & - clay(mstype), & - sand(mstype), & - swilt(mstype), & - sfc(mstype), & - ssat(mstype), & - bch(mstype), & - hyds(mstype), & - sucs(mstype), & - rhosoil(mstype), & - css(mstype) -END TYPE soilin_type - -!Instantiate types -TYPE(vegin_type) :: vegin !read from namelist -TYPE(soilin_type) :: soilin !read from namelist - -END MODULE cable_params_mod From 3178b8b204879d41665136e22c8cc09ed8197072 Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Tue, 21 Jan 2025 12:19:30 +1100 Subject: [PATCH 14/16] this spec init should possibly go completely --- src/science/soilsnow/cbl_soilsnow_init_special.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/science/soilsnow/cbl_soilsnow_init_special.F90 b/src/science/soilsnow/cbl_soilsnow_init_special.F90 index ea6fafbd2..0a7a14014 100644 --- a/src/science/soilsnow/cbl_soilsnow_init_special.F90 +++ b/src/science/soilsnow/cbl_soilsnow_init_special.F90 @@ -31,11 +31,11 @@ SUBROUTINE spec_init_soil_snow(dels, soil, ssnow, canopy, met, bal, veg,heat_cap REAL :: heat_cap_lower_limit(mp,ms) ktau = ktau +1 - -IF( .NOT.cable_user%cable_runtime_coupled ) THEN +! since CMIP5 only ever (potentially) TRUE offline if special initialization +IF( .NOT.cable_user%cable_runtime_coupled ) THEN IF( ktau_gl <= 1 ) THEN - IF (cable_runtime%um) canopy%dgdtg = 0.0 ! RML added um condition + ! after discussion with BP ! N.B. snmin should exceed sum of layer depths, i.e. .11 m ssnow%wbtot = 0.0 @@ -78,6 +78,8 @@ SUBROUTINE spec_init_soil_snow(dels, soil, ssnow, canopy, met, bal, veg,heat_cap END IF ENDIF ! if(.NOT.cable_runtime_coupled) +! this is done in all cases? overwrites gammzz from above and mediates thru +! snowd% (but only for single layer snow) IF (ktau <= 1) THEN xx=heat_cap_lower_limit(:,1) ssnow%gammzz(:,1) = MAX( (1.0 - soil%ssat) * soil%css * soil%rhosoil & From 31a233fbbc39deea28b910a666e5fbff2029b54c Mon Sep 17 00:00:00 2001 From: "Srbinovsky, Jhan (O&A, Aspendale)" Date: Fri, 24 Jan 2025 09:05:04 +1100 Subject: [PATCH 15/16] gammz initialization required --- .../soilsnow/cbl_soilsnow_init_special.F90 | 21 ++++++++----------- 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/src/science/soilsnow/cbl_soilsnow_init_special.F90 b/src/science/soilsnow/cbl_soilsnow_init_special.F90 index 0a7a14014..69e25503e 100644 --- a/src/science/soilsnow/cbl_soilsnow_init_special.F90 +++ b/src/science/soilsnow/cbl_soilsnow_init_special.F90 @@ -27,15 +27,12 @@ SUBROUTINE spec_init_soil_snow(dels, soil, ssnow, canopy, met, bal, veg,heat_cap REAL, DIMENSION(mp) :: xx, tgg_old, tggsn_old REAL(r_2), DIMENSION(mp) :: xxx,deltat,sinfil1,sinfil2,sinfil3 REAL :: zsetot -INTEGER, SAVE :: ktau =0 REAL :: heat_cap_lower_limit(mp,ms) -ktau = ktau +1 ! since CMIP5 only ever (potentially) TRUE offline if special initialization +! requested and then only on first timestep IF( .NOT.cable_user%cable_runtime_coupled ) THEN - IF( ktau_gl <= 1 ) THEN - ! after discussion with BP ! N.B. snmin should exceed sum of layer depths, i.e. .11 m ssnow%wbtot = 0.0 @@ -75,18 +72,18 @@ SUBROUTINE spec_init_soil_snow(dels, soil, ssnow, canopy, met, bal, veg,heat_cap 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_ice, xx ) * soil%zse(1) - END IF + ENDIF ! if(.NOT.cable_runtime_coupled) ! this is done in all cases? overwrites gammzz from above and mediates thru ! snowd% (but only for single layer snow) -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_ice, xx ) * soil%zse(1) + & - & (1. - ssnow%isflag) * Ccgsnow * ssnow%snowd -END IF + +! Evaluated on first timestep offline/ ESM1.5 +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_ice, xx ) * soil%zse(1) + & + & (1. - ssnow%isflag) * Ccgsnow * ssnow%snowd END SUBROUTINE spec_init_soil_snow From c71bee41a7315f5b57a824aae9863e6468326f4b Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Fri, 14 Feb 2025 10:16:04 +1100 Subject: [PATCH 16/16] (#494) - Remove incorrect densities in water balance calc. User ssnow%wbtot to avoid approximating the ice density with the liquid water density. --- src/offline/cable_checks.F90 | 56 +++++++++++++++++------------------- 1 file changed, 27 insertions(+), 29 deletions(-) diff --git a/src/offline/cable_checks.F90 b/src/offline/cable_checks.F90 index 5d1b73a5e..5ad32517c 100644 --- a/src/offline/cable_checks.F90 +++ b/src/offline/cable_checks.F90 @@ -474,48 +474,46 @@ END SUBROUTINE constant_check_range SUBROUTINE mass_balance(dels,ktau, ssnow,soil,canopy,met, & air,bal) + !* Calculates the energy and water balances at each time step as well as + ! the cumulative balance over the simulation. ! Input arguments - REAL,INTENT(IN) :: dels ! time step size - INTEGER, INTENT(IN) :: ktau ! timestep number - TYPE (soil_snow_type),INTENT(IN) :: ssnow ! soil data - TYPE (soil_parameter_type),INTENT(IN) :: soil ! soil data - TYPE (canopy_type),INTENT(IN) :: canopy ! canopy variable data - TYPE(met_type),INTENT(IN) :: met ! met data + REAL,INTENT(IN) :: dels + !! Time step length \(s\) + INTEGER, INTENT(IN) :: ktau + !! Number of time steps since the start of the simulation \(-\) + TYPE (soil_snow_type),INTENT(IN) :: ssnow + !! Soil and snow variables data + TYPE (soil_parameter_type),INTENT(IN) :: soil + !! Soil parameters + TYPE (canopy_type),INTENT(IN) :: canopy + !! Canopy variables data + TYPE(met_type),INTENT(IN) :: met + !! Met forcing data TYPE (air_type),INTENT(IN) :: air + !! Air variables data + TYPE (balances_type),INTENT(INOUT) :: bal + !! Water balance variables data ! Local variables - REAL(r_2), DIMENSION(:,:,:),POINTER, SAVE :: bwb ! volumetric soil moisture + REAL(r_2), DIMENSION(:),POINTER, SAVE :: owb ! volumetric soil moisture + ! at previous time step REAL(r_2), DIMENSION(mp) :: delwb ! change in soilmoisture ! b/w tsteps - REAL, DIMENSION(mp) :: canopy_wbal !canopy water balance - TYPE (balances_type),INTENT(INOUT) :: bal - INTEGER :: j, k ! do loop counter + REAL, DIMENSION(mp) :: canopy_wbal !canopy water balance + INTEGER :: j, k ! do loop counter IF(ktau==1) THEN - ALLOCATE( bwb(mp,ms,2) ) + ALLOCATE( owb(mp) ) ! initial vlaue of soil moisture - bwb(:,:,1)=ssnow%wb - bwb(:,:,2)=ssnow%wb - delwb(:) = 0. - ELSE - ! Calculate change in soil moisture b/w timesteps: - IF(MOD(REAL(ktau),2.0)==1.0) THEN ! if odd timestep - bwb(:,:,1)=ssnow%wb - DO k=1,mp ! current smoist - prev tstep smoist - delwb(k) = SUM((bwb(k,:,1) & - - (bwb(k,:,2)))*soil%zse)*1000.0 - END DO - ELSE IF(MOD(REAL(ktau),2.0)==0.0) THEN ! if even timestep - bwb(:,:,2)=ssnow%wb - DO k=1,mp ! current smoist - prev tstep smoist - delwb(k) = SUM((bwb(k,:,2) & - - (bwb(k,:,1)))*soil%zse)*1000.0 - END DO - END IF + owb=ssnow%wbtot END IF + ! Calculate change in soil moisture b/w timesteps: + delwb = ssnow%wbtot - owb + ! Update old soil moisture to this timestep value. + owb = ssnow%wbtot ! net water into soil (precip-(change in canopy water storage) ! - (change in snow depth) - (surface runoff) - (deep drainage)