Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8316)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8325)
@@ -104,5 +104,5 @@
       INTEGER , POINTER, DIMENSION(:,:) ::   icount    ! number of layers vanished by melting 
 
-      REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2)
+      REAL(wp), POINTER, DIMENSION(:) ::   zeh_i       ! total ice heat content  (J.m-2)
       REAL(wp), POINTER, DIMENSION(:) ::   zsnw        ! distribution of snow after wind blowing
 
@@ -121,5 +121,5 @@
 
       CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema )
-      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i )
+      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zeh_i )
       CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i )
       CALL wrk_alloc( jpij, nlay_i, icount )
@@ -127,5 +127,5 @@
       zqprec   (:) = 0._wp ; zq_su    (:) = 0._wp ; zq_bo    (:) = 0._wp ; zf_tt(:) = 0._wp
       zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp ; zevap_rema(:) = 0._wp ;
-      zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zqh_i(:) = 0._wp
+      zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zeh_i(:) = 0._wp
 
       zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp       
@@ -134,14 +134,14 @@
       ! Initialize enthalpy at nlay_i+1
       DO ji = kideb, kiut
-         q_i_1d(ji,nlay_i+1) = 0._wp
+         e_i_1d(ji,nlay_i+1) = 0._wp
       END DO
 
       ! initialize layer thicknesses and enthalpies
       h_i_old (:,0:nlay_i+1) = 0._wp
-      qh_i_old(:,0:nlay_i+1) = 0._wp
+      eh_i_old(:,0:nlay_i+1) = 0._wp
       DO jk = 1, nlay_i
          DO ji = kideb, kiut
             h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i
-            qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk)
+            eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk)
          ENDDO
       ENDDO
@@ -167,10 +167,10 @@
          IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting
             ! Contribution to heat flux to the ocean [W.m-2], < 0  
-            hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
+            hfx_res_1d(ji) = hfx_res_1d(ji) + e_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
             ! Contribution to mass flux
             wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
             ! updates
             ht_s_1d(ji)   = 0._wp
-            q_s_1d (ji,1) = 0._wp
+            e_s_1d (ji,1) = 0._wp
             t_s_1d (ji,1) = rt0
          END IF
@@ -184,5 +184,5 @@
          DO ji = kideb, kiut
             zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i
-            zqh_i(ji)   = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk)
+            zeh_i(ji)   = zeh_i(ji) + e_i_1d(ji,jk) * zh_i(ji,jk)
          END DO
       END DO
@@ -248,14 +248,14 @@
             ! thickness change
             rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 
-            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) ) 
-            zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 )
+            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) ) 
+            zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( e_s_1d(ji,jk), epsi20 )
             zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting
             zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)    
             ! heat used to melt snow(W.m-2, >0)
-            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice 
+            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d(ji,jk) * r1_rdtice 
             ! snow melting only = water into the ocean (then without snow precip)
             wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
             ! updates available heat + thickness
-            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )
+            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
             ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) )
          END DO
@@ -274,5 +274,5 @@
          ! Heat flux by sublimation [W.m-2], < 0 (sublimate first snow that had fallen, then pre-existing snow)
          zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
-         hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1)  &
+         hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1)  &
             &                              ) * a_i_1d(ji) * r1_rdtice
          ! Mass flux by sublimation
@@ -298,5 +298,5 @@
          DO ji = kideb,kiut
             rswitch       = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) )
-            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *           &
+            e_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *           &
               &            ( ( zdh_s_pre(ji)               ) * zqprec(ji) +  &
               &              ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
@@ -314,5 +314,5 @@
             IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
 
-               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]       
+               zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]       
                zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0)
                                                                       ! set up at 0 since no energy is needed to melt water...(it is already melted)
@@ -336,5 +336,5 @@
             ELSE                                !!! Surface melting
                
-               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]
+               zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]
                zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0]
                zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
@@ -377,5 +377,5 @@
             sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * sm_i_1d(ji) * r1_rdtice
             ! Heat flux [W.m-2], < 0
-            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * q_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice
+            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice
             ! Mass flux > 0
             wfx_ice_sub_1d(ji) =  wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice
@@ -392,5 +392,5 @@
                         
             ! update heat content (J.m-2) and layer thickness
-            qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
+            eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
             h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
          END DO
@@ -464,5 +464,5 @@
                dh_i_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
 
-               q_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
+               e_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
                
             ENDIF
@@ -502,5 +502,5 @@
 
             ! update heat content (J.m-2) and layer thickness
-            qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_1d(ji,nlay_i+1)
+            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * e_i_1d(ji,nlay_i+1)
             h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
          ENDIF
@@ -519,5 +519,5 @@
                IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
 
-                  zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
+                  zEi               = - e_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
                   zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
                                                                     ! set up at 0 since no energy is needed to melt water...(it is already melted)
@@ -539,10 +539,10 @@
 
                   ! update heat content (J.m-2) and layer thickness
-                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
+                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
                   h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
 
                ELSE                               !!! Basal melting
 
-                  zEi             = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
+                  zEi             = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
                   zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0)
                   zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0)
@@ -575,5 +575,5 @@
 
                   ! update heat content (J.m-2) and layer thickness
-                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
+                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
                   h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
                ENDIF
@@ -599,13 +599,13 @@
          zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
          rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow
-         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) )
-         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 )
+         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) )
+         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 )
          zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting
          dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1)
          ht_s_1d   (ji)  = ht_s_1d(ji)   + zdeltah(ji,1)
         
-         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1)                ! update available heat (J.m-2)
+         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                ! update available heat (J.m-2)
          ! heat used to melt snow
-         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * q_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0)
+         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0)
          ! Contribution to mass flux
          wfx_snw_sum_1d(ji)  =  wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
@@ -661,5 +661,5 @@
 
          ! update heat content (J.m-2) and layer thickness
-         qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw
+         eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw
          h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
          
@@ -679,7 +679,7 @@
             ! mask enthalpy
             rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  )
-            q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk)
-            ! recalculate t_s_1d from q_s_1d
-            t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
+            e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk)
+            ! recalculate t_s_1d from e_s_1d
+            t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
          END DO
       END DO
@@ -689,5 +689,5 @@
       
       CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema )
-      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i )
+      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zeh_i )
       CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i )
       CALL wrk_dealloc( jpij, nlay_i, icount )
