Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 6060)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 6069)
@@ -45,4 +45,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE iom
 
    IMPLICIT NONE
@@ -129,4 +130,10 @@
       INTEGER ::   ioptio = 0      ! temporary integer
       INTEGER ::   ios             ! Local integer output status for namelist read
+      !!
+      INTEGER  ::   ji, jj, jk, ikt    ! dummy loop indices      ISF
+      REAL(wp) ::   znad
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop, zrhd ! hypothesys on isf density
+      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_isf  ! density at bottom of ISF
+      REAL(wp), POINTER, DIMENSION(:,:)     ::  ziceload     ! density at bottom of ISF
       !!
       NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     &
@@ -189,6 +196,51 @@
       IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
       ! 
-      ! initialisation of ice load
-      riceload(:,:)=0.0
+      ! initialisation of ice shelf load
+      IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0
+      IF (       ln_isfcav ) THEN
+         CALL wrk_alloc( jpi,jpj, 2,  ztstop) 
+         CALL wrk_alloc( jpi,jpj,jpk, zrhd  )
+         CALL wrk_alloc( jpi,jpj,     zrhdtop_isf, ziceload) 
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'  
+
+         ! To use density and not density anomaly
+         znad=1._wp
+         
+         ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
+         ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp
+
+         ! compute density of the water displaced by the ice shelf 
+         DO jk = 1, jpk
+            CALL eos(ztstop(:,:,:),fsdept_n(:,:,jk),zrhd(:,:,jk))
+         END DO
+      
+         ! compute rhd at the ice/oce interface (ice shelf side)
+         CALL eos(ztstop,risfdep,zrhdtop_isf)
+
+         ! Surface value + ice shelf gradient
+         ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v)
+         ! divided by 2 later
+         ziceload = 0._wp
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ikt=mikt(ji,jj)
+               ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1))
+               DO jk=2,ikt-1
+                  ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * fse3w(ji,jj,jk) &
+                     &                              * (1._wp - tmask(ji,jj,jk))
+               END DO
+               IF (ikt  >=  2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &
+                                  &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) )
+            END DO
+         END DO
+         riceload(:,:)=ziceload(:,:)  ! need to be saved for diaar5
+
+         CALL wrk_dealloc( jpi,jpj, 2,  ztstop) 
+         CALL wrk_dealloc( jpi,jpj,jpk, zrhd  )
+         CALL wrk_dealloc( jpi,jpj,     zrhdtop_isf, ziceload) 
+      END IF
       !
    END SUBROUTINE dyn_hpg_init
@@ -444,5 +496,5 @@
    SUBROUTINE hpg_isf( kt )
       !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE hpg_sco  ***
+      !!                  ***  ROUTINE hpg_isf  ***
       !!
       !! ** Method  :   s-coordinate case. Jacobian scheme.
@@ -463,82 +515,38 @@
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
-      INTEGER  ::   ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j                 ! dummy loop indices
-      REAL(wp) ::   zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1             ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj, zrhd
+      INTEGER  ::   ji, jj, jk, ikt, iktp1i, iktp1j   ! dummy loop indices
+      REAL(wp) ::   zcoef0, zuap, zvap, znad          ! temporary scalars
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop
-      REAL(wp), POINTER, DIMENSION(:,:)     ::  ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj,  2,   ztstop ) 
-      CALL wrk_alloc( jpi,jpj,jpk,   zhpi, zhpj, zrhd)
-      CALL wrk_alloc( jpi,jpj,       ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj ) 
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used'
-      ENDIF
-      !
+      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_oce
+      !!----------------------------------------------------------------------
+      !
+      CALL wrk_alloc( jpi,jpj,  2, ztstop) 
+      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj)
+      CALL wrk_alloc( jpi,jpj,     zrhdtop_oce )
+      !
+      ! Local constant initialization
       zcoef0 = - grav * 0.5_wp
-      IF( ln_linssh ) THEN   ;   znad = 0._wp      ! Fixed    volume: density anomaly
-      ELSE                   ;   znad = 1._wp      ! Variable volume: density
-      ENDIF
-      zhpi(:,:,:) = 0._wp
-      zhpj(:,:,:) = 0._wp
+  
+      ! To use density and not density anomaly
+      znad=1._wp
+
+      ! iniitialised to 0. zhpi zhpi 
+      zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp
+
+      ! compute rhd at the ice/oce interface (ocean side)
+      ! usefull to reduce residual current in the test case ISOMIP with no melting
+      DO ji=1,jpi
+        DO jj=1,jpj
+          ikt=mikt(ji,jj)
+          ztstop(ji,jj,1)=tsn(ji,jj,ikt,1)
+          ztstop(ji,jj,2)=tsn(ji,jj,ikt,2)
+        END DO
+      END DO
+      CALL eos( ztstop, risfdep, zrhdtop_oce )
 
 !==================================================================================     
-!=====Compute iceload and contribution of the half first wet layer ================= 
-!===================================================================================
-
-      ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
-      ztstop(:,:,jp_tem) = -1.9_wp
-      ztstop(:,:,jp_sal) = 34.4_wp
-
-!!gm I have the feeling that a much simplier and faster computation can be performed...
-!!gm     ====>>>>   We have to discuss !
-
-!!gm below, faster to compute the ISF density in zrhd and remplace rhd value where tmask=0
-!!gm        furthermore, this calculation does not depends on time :  do it at the first time-step only....
-
-      ! compute density of the water displaced by the ice shelf 
-      zrhd(:,:,:) = rhd(:,:,:)     ! save rhd
-      DO jk = 1, jpk
-         zdept(:,:) = gdept_1d(jk)
-         CALL eos( ztstop(:,:,:), zdept(:,:), rhd(:,:,jk) )
-      END DO
-      WHERE( tmask(:,:,:) == 1._wp )
-        rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd
-      END WHERE
-      
-      ! compute rhd at the ice/oce interface (ice shelf side)
-      CALL eos( ztstop, risfdep, zrhdtop_isf )
-
-      ! compute rhd at the ice/oce interface (ocean side)
-      DO ji = 1, jpi
-        DO jj = 1, jpj
-          ikt = mikt(ji,jj)
-          ztstop(ji,jj,jp_tem) = tsn(ji,jj,ikt,jp_tem)
-          ztstop(ji,jj,jp_sal) = tsn(ji,jj,ikt,jp_sal)
-        END DO
-      END DO
-      CALL eos( ztstop, risfdep, zrhdtop_oce )
-      !
-      ! Surface value + ice shelf gradient
-      ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v)
-      ziceload = 0._wp
-      DO jj = 1, jpj
-         DO ji = 1, jpi   ! vector opt.
-            ikt = mikt(ji,jj)
-            ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1))
-            DO jk = 2, ikt-1
-               ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * e3w_n(ji,jj,jk) &
-                  &                              * (1._wp - tmask(ji,jj,jk))
-            END DO
-            IF( ikt >= 2 )   ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) &
-               &                                               * ( risfdep(ji,jj) - gdept_1d(ikt-1) )
-         END DO
-      END DO
-      riceload(:,:) = 0._wp   ;   riceload(:,:) = ziceload(:,:)   ! need to be saved for diaar5
-      ! compute zp from z=0 to first T wet point (correction due to zps not yet applied)
+!===== Compute surface value ===================================================== 
+!==================================================================================
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -548,12 +556,14 @@
             ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure
             ! we assume ISF is in isostatic equilibrium
-            zhpi(ji,jj,1) = zcoef0 * (                                    &
-               &            0.5_wp * e3w_n(ji+1,jj,iktp1i) * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) )   &
-               &          - 0.5_wp * e3w_n(ji  ,jj,ikt   ) * ( 2._wp * znad + rhd(ji  ,jj,ikt   ) + zrhdtop_oce(ji  ,jj) )   &
-               &          + ( ziceload(ji+1,jj) - ziceload(ji,jj) )       ) * r1_e1u(ji,jj)
-            zhpj(ji,jj,1) = zcoef0 * (                                    &
-               &            0.5_wp * e3w_n(ji,jj+1,iktp1j) * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   &
-               &          - 0.5_wp * e3w_n(ji,jj  ,ikt   ) * ( 2._wp * znad + rhd(ji,jj  ,ikt   ) + zrhdtop_oce(ji,jj  ) )   &
-               &          + ( ziceload(ji,jj+1) - ziceload(ji,jj) )       ) * r1_e2v(ji,jj)
+            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj,iktp1i)                                    &
+               &                                    * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) )   &
+               &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         &
+               &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
+               &                                  + ( riceload(ji+1,jj) - riceload(ji,jj))                            ) 
+            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji,jj+1,iktp1j)                                    &
+               &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   &
+               &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         & 
+               &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
+               &                                  + ( riceload(ji,jj+1) - riceload(ji,jj))                            ) 
             ! s-coordinate pressure gradient correction (=0 if z coordinate)
             zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
@@ -567,140 +577,32 @@
       END DO
 !==================================================================================     
-!===== Compute partial cell contribution for the top cell ========================= 
-!==================================================================================
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            iku = miku(ji,jj)
-            zpshpi(ji,jj) = 0._wp
-            zpshpj(ji,jj) = 0._wp
-            ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
-            ! u direction
-            IF( iku > 1 ) THEN
-               ! case iku
-               zhpi(ji,jj,iku) = zcoef0 * r1_e1u(ji,jj) * ze3wu                             &
-                  &                     * ( rhd(ji+1,jj,iku) + rhd(ji,jj,iku)               &
-                  &                        + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )
-               ! corrective term ( = 0 if z coordinate )
-               zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) * r1_e1u(ji,jj)
-               ! zhpi will be added in interior loop
-               ua(ji,jj,iku) = ua(ji,jj,iku) + zuap
-               ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure  
-               IF( mbku(ji,jj) == iku + 1 )   zpshpi(ji,jj) = zhpi(ji,jj,iku)
-
-               ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps)
-               zhpiint = zcoef0 * r1_e1u(ji,jj)                                                              &    
-                  &    * (  e3w_n(ji+1,jj  ,iku+1) * (  (rhd(ji+1,jj,iku+1) + znad)                          &
-                  &                                   + (rhd(ji+1,jj,iku  ) + znad) ) * tmask(ji+1,jj,iku)   &
-                  &       - e3w_n(ji  ,jj  ,iku+1) * (  (rhd(ji  ,jj,iku+1) + znad)                          &
-                  &                                   + (rhd(ji  ,jj,iku  ) + znad) ) * tmask(ji  ,jj,iku)   )
-               zhpi(ji,jj,iku+1) = zcoef0 * r1_e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 
-            END IF
-               
-            ! v direction
-            ikv = mikv(ji,jj)
-            ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
-            IF( ikv > 1 ) THEN
-               ! case ikv
-               zhpj(ji,jj,ikv) =  zcoef0 * r1_e2v(ji,jj) * ze3wv                             &
-                  &                      * ( rhd(ji,jj+1,ikv) + rhd(ji,jj,ikv)               &
-                  &                         + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )
-               ! corrective term ( = 0 if z coordinate )
-               zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) * r1_e2v(ji,jj)
-               ! zhpi will be added in interior loop
-               va(ji,jj,ikv) = va(ji,jj,ikv) + zvap
-               ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure  
-               IF( mbkv(ji,jj) == ikv + 1 )   zpshpj(ji,jj) = zhpj(ji,jj,ikv) 
-               
-               ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps)
-               zhpjint =  zcoef0 * r1_e2v(ji,jj)                                                            &
-                  &    * (  e3w_n(ji  ,jj+1,ikv+1) * (  (rhd(ji,jj+1,ikv+1) + znad)                         &
-                  &                                   + (rhd(ji,jj+1,ikv  ) + znad) ) * tmask(ji,jj+1,ikv)  &
-                  &       - e3w_n(ji  ,jj  ,ikv+1) * (  (rhd(ji,jj  ,ikv+1) + znad)                         &
-                  &                                   + (rhd(ji,jj  ,ikv  ) + znad) ) * tmask(ji,jj  ,ikv)  )
-               zhpj(ji,jj,ikv+1) =  zcoef0 * r1_e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint
-            ENDIF
-         END DO
-      END DO
-
-!==================================================================================     
 !===== Compute interior value ===================================================== 
 !==================================================================================
-
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            DO jk = 2, jpkm1
+      ! interior value (2=<jk=<jpkm1)
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
                ! hydrostatic pressure gradient along s-surfaces
-               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc)
-               zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1)                                                     &
-                  &           + zcoef0 * r1_e1u(ji,jj)                                                                &
-                  &                    * ( e3w_n(ji+1,jj,jk) * ( (rhd(ji+1,jj,jk  ) + znad)                           &
-                  &                                            + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1)   &
-                  &                      - e3w_n(ji  ,jj,jk) * ( (rhd(ji  ,jj,jk  ) + znad)                           &
-                  &                                            + (rhd(ji  ,jj,jk-1) + znad) ) * tmask(ji  ,jj,jk-1)   ) 
+               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
+                  &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   &
+                  &              - e3w_n(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   )
+               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
+                  &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   &
+                  &              - e3w_n(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   )
                ! s-coordinate pressure gradient correction
-               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 
-               zuap = - zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )                    &
-                  &            * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk-1)
-               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
-
-               ! hydrostatic pressure gradient along s-surfaces
-               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc)
-               zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1)                                                     &
-                  &           + zcoef0 * r1_e2v(ji,jj)                                                                &
-                  &                    * ( e3w_n(ji  ,jj+1,jk) * ( (rhd(ji,jj+1,jk  ) + znad)                           &
-                  &                                              + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1)   &
-                  &                      - e3w_n(ji  ,jj  ,jk) * ( (rhd(ji,jj  ,jk  ) + znad)                           &
-                  &                                              + (rhd(ji,jj  ,jk-1) + znad) ) * tmask(ji,jj  ,jk-1)   )
-               ! s-coordinate pressure gradient correction
-               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)
-               zvap = - zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )                    &
-                  &            * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk-1)
+               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
+                  &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj)
+               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
+                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj)
                ! add to the general momentum trend
-               va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk)
+               ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
-
-!==================================================================================     
-!===== Compute bottom cell contribution (partial cell) ============================ 
-!==================================================================================
-
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
-            iku = mbku(ji,jj)
-            ikv = mbkv(ji,jj)
-
-            IF (iku .GT. 1) THEN
-               ! remove old value (interior case)
-               zuap            = -zcoef0 * ( rhd    (ji+1,jj  ,iku) + rhd    (ji,jj,iku) + 2._wp * znad )   &
-                     &                   * ( gde3w_n(ji+1,jj  ,iku) - gde3w_n(ji,jj,iku) ) * r1_e1u(ji,jj)
-               ua(ji,jj,iku)   = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap
-               ! put new value
-               ! -zpshpi to avoid double contribution of the partial step in the top layer 
-               zuap            = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj)  * r1_e1u(ji,jj)
-               zhpi(ji,jj,iku) =  zhpi(ji,jj,iku-1) + zcoef0 * r1_e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 
-               ua(ji,jj,iku)   =  ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap
-            END IF
-            ! v direction
-            IF (ikv .GT. 1) THEN
-               ! remove old value (interior case)
-               zvap            = -zcoef0 * ( rhd    (ji  ,jj+1,ikv) + rhd    (ji,jj,ikv) + 2._wp * znad )   &
-                     &                   * ( gde3w_n(ji  ,jj+1,ikv) - gde3w_n(ji,jj,ikv) ) * r1_e2v(ji,jj)
-               va(ji,jj,ikv)   = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap
-               ! put new value
-               ! -zpshpj to avoid double contribution of the partial step in the top layer 
-               zvap            = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) * r1_e2v(ji,jj)
-               zhpj(ji,jj,ikv) =  zhpj(ji,jj,ikv-1) + zcoef0 * r1_e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj)   
-               va(ji,jj,ikv)   =  va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap
-            END IF
-         END DO
-      END DO
-     
-      ! set back to original density value into the ice shelf cell (maybe useless because it is masked)
-      rhd = zrhd
-      !
-      CALL wrk_dealloc( jpi,jpj,2, ztstop)
-      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd)
-      CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj)
+     !
+      CALL wrk_dealloc( jpi,jpj,2  , ztstop)
+      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj)
+      CALL wrk_dealloc( jpi,jpj    , zrhdtop_oce )
       !
    END SUBROUTINE hpg_isf
