Changeset 27 for branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_th_diff.f

Timestamp:

10/06/16 11:31:57 (8 years ago)

Author:

vancop

Message:

First physical fixes, new ice liquid fraction and mushy-layer thermal conductivity

File:

• branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_th_diff.f (modified) (5 diffs)

branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_th_diff.f

@@ -59 +59 @@
       INTEGER    numeqmin, numeqmax, numeq
       DIMENSION  ztcond_i(0:nlay_i),
+     &           ze(0:nlay_i),
+     &           ztc(0:nlay_i),
      &           zkappa_s(0:nlay_s),
      &           zkappa_i(0:nlay_i),ztstemp(0:nlay_s),ztitemp(0:nlay_i),
@@ -122 +124 @@
 !------------------------------------------------------------------------------ 
 !
-      ! Pringle et al., JGR 2007 formula
-      ! 2.11 + 0.09 S/T - 0.011.T
+      nn_thcond = 2
 
       DO 10 ji = kideb, kiut
 
-      ! thermal conductivity in the snow
+      ! thermal conductivity of the snow
       ykn(ji)  =  xkn
       zkimin   =  0.1
+
+      ! thermal conductivity of the sea ice
+
+      SELECT CASE ( nn_thcond )
+
+      !---------------------------------------------------
+      CASE(1) ! Pringle et al., JGR 2007
+              ! 2.11 + 0.09 S/T - 0.011.T
+              ! Empirical, direct measurements in sea ice
+      !---------------------------------------------------
 
       ztcond_i(0)     = xkg + betak1*s_i_b(ji,1) 
      &                   / MIN( -zeps , t_i_b(ji,1) - tpw ) 
      &                   - betak2* ( t_i_b(ji,1) - tpw )
-      ztcond_i(0)     = MAX( ztcond_i(0) , zkimin )
 
       DO layer = 1, nlay_i-1
@@ -142 +152 @@
      &                      - betak2 * 0.5 * ( t_i_b(ji,layer) +  ! bugfix fred dupont add 0.5
      &                        t_i_b(ji,layer+1) - 2.0*tpw )
-         ztcond_i(layer) = MAX( ztcond_i(layer) , zkimin )
       END DO
 
@@ -148 +157 @@
      &                        MIN( -zeps , t_bo_b(ji) - tpw ) 
      &                      - betak2 * ( t_bo_b(ji) - tpw )
+
+      !---------------------------------------------------
+      CASE(2) ! Mushy layer: k = ki . ( 1 - phi) + kbr . phi
+              ! Thermal conductivity of pure ice at interfaces from Yen (1991)
+              ! zki = 2.21 - 1.0e-2*ztc + 3.44e-5*ztc.^2 (Yen 1991) 
+              ! Thermal conductivity of brine at interfaces (Castelli et al DSR 1974)
+              ! 0.55286 + 1.8364e-3 * TT - 3.3058e-7 * TT.^3
+      !---------------------------------------------------
+
+      ! Brine fraction at layer mid-points
+      e_i_b(:) = -tmut * s_i_b(ji,:) / ( t_i_b(ji,:) - tpw )
+
+      ! Brine fraction at interfaces
+      zde_dz       = ( e_i_b(2)   - e_i_b(1) ) / 
+     &               ( z_i_phy(2) - z_i_phy(1)  )
+      ze(0)        = e_i_b(1)     - zde_dz * deltaz_i_phy(1) / 2.
+
+      DO layer = 1, nlay_i - 1
+         zde_dz    = ( e_i_b(layer+1) - e_i_b(layer) ) / 
+     &               ( z_i_phy(layer+1) - z_i_phy(layer)  )
+         ze(layer) = e_i_b(layer)    + zde_dz * deltaz_i_phy(layer) / 2.
+      END DO
+
+      ze(nlay_i) = 1.
+
+      ! Temperature at interfaces
+      !!! ztc(0) could be better interpolated ... use effective conductivity ?
+      zdt_dz       = ( t_i_b(ji,2)     - t_i_b(ji,1) ) / 
+     &               ( z_i_phy(2)      - z_i_phy(1)  )
+      ztc(0)       = t_i_b(ji,1)       - zdt_dz * deltaz_i_phy(1) / 2.
+      ztc(0)       = ztc(0) - tpw
+
+      DO layer = 1, nlay_i - 1
+         zdt_dz    = ( t_i_b(ji,layer+1) - t_i_b(ji,layer) ) / 
+     &               ( z_i_phy(layer+1)  - z_i_phy(layer)  )
+         ztc(layer) = t_i_b(ji,layer) + 
+     &                zdt_dz * deltaz_i_phy(layer) / 2.
+         ztc(layer) = ztc(layer) - tpw
+      END DO
+      ztc(nlay_i)    = t_bo_b(ji) - tpw
+
+      DO layer = 0, nlay_i
+         zt1 = ztc(layer)
+         zt2 = zt1*zt1
+         zt3 = zt1*zt2
+         zki = 2.2156 - 1.0045e-2*zt1 + 3.44520e-5*ztc2
+         zkbr = 0.55286 + 1.8364e-3*zt1 - 3.3058e-7*zt3
+         ztcond_i(layer) = zki * ( 1. - ze(layer) ) + zkbr * ze(layer)
+      END DO
+
+      END SELECT
+
+      ! Cap thermal conductivity to prevent small values
+      ztcond_i(0)     = MAX( ztcond_i(0) , zkimin )
+      DO layer = 0, nlay_i
+         ztcond_i(layer) = MAX( ztcond_i(layer) , zkimin )
+      END DO
       ztcond_i(nlay_i)   = MAX( ztcond_i(nlay_i) , zkimin )
 
@@ -488 +554 @@
          zindterm(numeqmin) =  ztiold(1) + zeta_i(1)*
 !    &                         (radab_phy_i(1) + radab_alg_i(1)
-     &                         ( radab_i(1) + 
+     &                         ( radab_i(1)   
      &                         + zkappa_i(1)*t_bo_b(ji) ) 
      &                         + t_su_b(ji)*zeta_i(1)*zkappa_i(0)*2.0