494 redux local ice fixes

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

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
 
 

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

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)
 
@@ -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
 

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 )
 

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

src/science/soilsnow/cbl_soilsnow_init_special.F90

@@ -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