Changeset 13810

Timestamp:

2020-11-18T11:17:25+01:00 (4 years ago)

Author:

vancop

Message:

Some quick fixes

File:

• NEMO/branches/2020/SI3-05_MP/src/ICE/icethd_pnd.F90 (modified) (14 diffs)

NEMO/branches/2020/SI3-05_MP/src/ICE/icethd_pnd.F90

@@ -277 +277 @@
             
             ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010)
-            ! MV the expression for zsbr has wrong sign!!!!
-            ! MV note - the sign-corrected expression is inconsistent 
-            ! with the rest of the SI3 code which assumes linear liquidus
-            ! best expression (most consistent for SI3 should be)
             ! zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt
             DO jk = 1, nlay_i
@@ -287 +283 @@
                   &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 &
                   &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3
-               
-               ztmp(jk) = sz_i_1d(ji,jk) / zsbr ! MV -> this is brine fraction and as it reads is always negative
+!               zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt
+               ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
             END DO
-            zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) ! MV at present since ztmp is negative, this is always zero
+            
+            zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 
             
             ! Do the drainage using Darcy's law
@@ -504 +501 @@
 
       !js 03/05/19: we truncate negative values after calculating zvolp, in a
-      ! similar manner to the subroutine lim_mp_cesm. Variation dh_i_pnd and
+      ! similar manner to the subroutine ice_thd_pnd_cesm. Variation dh_i_pnd and
       ! dh_s_pnd are negative, indicating a loss of ice or snow. But we can expect them
       ! to be negative for some reasons. We keep this behaviour as it is, for
@@ -582 +579 @@
           ! OLI 07/2017: Done like for empirical melt pond scheme (CESM).
           ! Therefore, I chose to put this part of the code before the main
-          ! routines lim_mp_area/depth (contrary to the original code),
+          ! routines ice_thd_pnd_area/depth (contrary to the original code),
           ! seeing the freeze-up as a global sink of
           ! freshwater for melt ponds in the whole grid cell. If this was done
@@ -626 +623 @@
           !--------------------------------------------------------------
           zdvn = 0._wp
-          CALL lim_mp_area(aice(ji,jj),vice(ji,jj), &
+          CALL ice_thd_pnd_area(aice(ji,jj),vice(ji,jj), &
                     aicen(ji,jj,:), vicen(ji,jj,:), vsnon(ji,jj,:), &
                     ticen(ji,jj,:,:), salin(ji,jj,:,:), &
@@ -663 +660 @@
    END SUBROUTINE pnd_TOPO
    
-       SUBROUTINE lim_mp_area(aice,  vice, &
+       SUBROUTINE ice_thd_pnd_area(aice,  vice, &
                            aicen, vicen, vsnon, ticen, &
                            salin, zvolpn, zvolp,         &
@@ -669 +666 @@
 
        !!-------------------------------------------------------------------
-       !!                ***  ROUTINE lim_mp_area ***
+       !!                ***  ROUTINE ice_thd_pnd_area ***
        !!
        !! ** Purpose : Given the total volume of meltwater, update
@@ -874 +871 @@
        ! height and area corresponding to the remaining volume
 
-      CALL lim_mp_depth(reduced_aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index)
+      CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index)
 
        DO n=1, m_index
@@ -910 +907 @@
              !IF (hicen(n) > z0) THEN           !js: from CICE 5.1.2
                 perm = 0._wp ! MV ugly dummy patch
-                CALL lim_mp_perm(ticen(:,n), salin(:,n), perm)
+                CALL ice_thd_pnd_perm(ticen(:,n), salin(:,n), perm)
                 IF (perm > z0) permflag = 1
 
@@ -927 +924 @@
           IF (permflag > 0) THEN
              ! recompute pond depth
-            CALL lim_mp_depth(reduced_aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index)
+            CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index)
              DO n=1, m_index
                 zhpondn(n) = hpond - alfan(n) + alfan(1)
@@ -989 +986 @@
        END DO
 
-    END SUBROUTINE lim_mp_area
-
-
-    SUBROUTINE lim_mp_depth(aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index)
+    END SUBROUTINE ice_thd_pnd_area
+
+
+    SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index)
        !!-------------------------------------------------------------------
-       !!                ***  ROUTINE lim_mp_depth  ***
+       !!                ***  ROUTINE ice_thd_pnd_depth  ***
        !!
        !! ** Purpose :   Compute melt pond depth
@@ -1154 +1151 @@
        END IF
 
-    END SUBROUTINE lim_mp_depth
-
-
-    SUBROUTINE lim_mp_perm(ticen, salin, perm)
+    END SUBROUTINE ice_thd_pnd_depth
+
+
+    SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm)
        !!-------------------------------------------------------------------
-       !!                ***  ROUTINE lim_mp_perm ***
+       !!                ***  ROUTINE ice_thd_pnd_perm ***
        !!
        !! ** Purpose :   Determine the liquid fraction of brine in the ice
@@ -1193 +1190 @@
        !-----------------------------------------------------------------
 
-       IF (maxval(Tin) <= -2.0_wp) THEN
-
-          ! Assur 1958
-          DO k = 1,nlay_i
-             Sbr = - 1.2_wp                 &
-                   - 21.8_wp    * Tin(k)    &
-                   - 0.919_wp   * Tin(k)**2 &
-                   - 0.01878_wp * Tin(k)**3
-             phi(k) = salin(k) / Sbr ! liquid fraction
-
-             !js: Sbr must not be zero
-         END DO ! k
-
-       ELSE
-
-          ! Notz 2005 thesis eq. 3.2
-          !js: "interstitial brine Sbr (in ppt) as a function of temperature T (in degC)".
-          DO k = 1,nlay_i
-             Sbr = - 17.6_wp   * Tin(k)    &
-                   - 0.389_wp  * Tin(k)**2 &
-                   - 0.00362_wp* Tin(k)**3
-             phi(k) = salin(k) / Sbr ! liquid fraction
-
-             !js: Sbr must not be zero
-          END DO
-
-       END IF
+       DO k = 1, nlay_i
+       
+          Sbr    = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus)
+          ! Best expression to date is that one
+          ! Sbr  = - 18.7 * Tin(k) − 0.519 * Tin(k)**2 − 0.00535 * Tin(k) ***3
+          phi(k) = salin(k) / Sbr
+          
+       END DO
 
        !-----------------------------------------------------------------
@@ -1227 +1205 @@
        perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007)
 
-   END SUBROUTINE lim_mp_perm
+   END SUBROUTINE ice_thd_pnd_perm