Changeset 12938 for NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icethd_zdf_bl99.F90

Timestamp:

2020-05-15T18:26:45+02:00 (4 years ago)

Author:

dancopsey

Message:

Merge in missing icethd_zdf_bl99 changes including:

• Make isnow a float
• Turn off 0.1m snow test
• Use conductive heat flux from coupler

File:

• NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icethd_zdf_bl99.F90 (modified) (7 diffs)

NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icethd_zdf_bl99.F90

@@ -101 +101 @@
       REAL(wp) ::   zfac                      ! dummy factor
       !
-      REAL(wp), DIMENSION(jpij) ::   isnow        ! switch for presence (1) or absence (0) of snow
+      REAL(wp), DIMENSION(jpij) ::   isnow        ! fraction of sea ice that is snow covered (for thermodynamic use only)
       REAL(wp), DIMENSION(jpij) ::   ztsub        ! surface temperature at previous iteration
       REAL(wp), DIMENSION(jpij) ::   zh_i, z1_h_i ! ice layer thickness
@@ -147 +147 @@
       !------------------
       DO ji = 1, npti
-         isnow(ji) = 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) )  ! is there snow or not
+
+         ! If the snow thickness drops below zhs_min then reduce the snow fraction instead
+         IF( h_s_1d(ji) < zhs_min ) THEN
+           isnow(ji) = h_s_1d(ji) / zhs_min
+           zh_s(ji) = zhs_min * r1_nlay_s
+         ELSE
+           isnow(ji) = 1.0_wp
+           zh_s(ji) = h_s_1d(ji) * r1_nlay_s
+         END IF
+
          ! layer thickness
          zh_i(ji) = h_i_1d(ji) * r1_nlay_i
-         zh_s(ji) = h_s_1d(ji) * r1_nlay_s
+
       END DO
       !
@@ -156 +165 @@
       ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp
       END WHERE
-      !
-      WHERE( zh_s(1:npti) > 0._wp   )       zh_s(1:npti) = MAX( zhs_min * r1_nlay_s, zh_s(1:npti) )
       !
       WHERE( zh_s(1:npti) > 0._wp   )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
@@ -769 +776 @@
       !
       ! --- calculate conduction fluxes (positive downward)
-
+      !     bottom ice conduction flux
       DO ji = 1, npti
-         !                                ! surface ice conduction flux
-         qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0)      * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  &
-            &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0)      * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) )
-         !                                ! bottom ice conduction flux
-         qcn_ice_bot_1d(ji) =                          - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) )
+         qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) )
       END DO
-      
+      !     surface ice conduction flux
+      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
+         !
+         DO ji = 1, npti
+            qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  &
+               &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) )
+         END DO
+         !
+      ELSEIF( k_cnd == np_cnd_ON ) THEN
+         !
+         DO ji = 1, npti
+            qcn_ice_top_1d(ji) = qcn_ice_1d(ji)
+         END DO
+         !
+      ENDIF
+      !     surface ice temperature
+      IF( k_cnd == np_cnd_ON .AND. ln_cndemulate ) THEN
+         !
+         DO ji = 1, npti
+            t_su_1d(ji) = (  qcn_ice_top_1d(ji) &            ! calculate surface temperature
+               &           +           isnow(ji)   * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) &
+               &           + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * t_i_1d(ji,1) &
+               &          ) / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 )
+            t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp )  ! cap t_su
+         END DO
+         !
+      ENDIF
       !
       ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- !
@@ -787 +816 @@
          END DO
          !
-      ELSEIF( k_cnd == np_cnd_ON ) THEN
-         !
-         DO ji = 1, npti
-            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qcn_ice_top_1d(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji) 
-         END DO
-         !
       ENDIF
-      
       !
       ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- !
@@ -839 +861 @@
       ! effective conductivity and 1st layer temperature (needed by Met Office)
       DO ji = 1, npti
-         IF( h_s_1d(ji) > 0.1_wp ) THEN 
-            cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0)
+         IF( h_i_1d(ji) > 0.1_wp ) THEN
+            cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
          ELSE
-            IF( h_i_1d(ji) > 0.1_wp ) THEN
-               cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
-            ELSE
-               cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp
-            ENDIF
+            cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp
          ENDIF
          t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
@@ -856 +874 @@
          t_i_1d    (1:npti,:) = ztiold        (1:npti,:)
          qcn_ice_1d(1:npti)   = qcn_ice_top_1d(1:npti)
-
-         !!clem
-         ! remettre t_su_1d, qns_ice_1d et dqns_ice_1d comme avant puisqu'on devrait faire comme si on avant conduction = input
-         !clem
       ENDIF
       !