Index: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5234)
+++ branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5443)
@@ -16,4 +16,5 @@
    !!            3.4  !  2011-11  (H. Liu) hpg_prj: Original code for s-coordinates
    !!                 !           (A. Coward) suppression of hel, wdj and rot options
+   !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity
    !!----------------------------------------------------------------------
 
@@ -25,4 +26,5 @@
    !!       hpg_zps  : z-coordinate plus partial steps (interpolation)
    !!       hpg_sco  : s-coordinate (standard jacobian formulation)
+   !!       hpg_isf  : s-coordinate (sco formulation) adapted to ice shelf
    !!       hpg_djc  : s-coordinate (Density Jacobian with Cubic polynomial)
    !!       hpg_prj  : s-coordinate (Pressure Jacobian with Cubic polynomial)
@@ -55,4 +57,5 @@
    LOGICAL , PUBLIC ::   ln_hpg_djc      !: s-coordinate (Density Jacobian with Cubic polynomial)
    LOGICAL , PUBLIC ::   ln_hpg_prj      !: s-coordinate (Pressure Jacobian scheme)
+   LOGICAL , PUBLIC ::   ln_hpg_isf      !: s-coordinate similar to sco modify for isf
    LOGICAL , PUBLIC ::   ln_dynhpg_imp   !: semi-implicite hpg flag
 
@@ -97,4 +100,5 @@
       CASE (  3 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial)
       CASE (  4 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme)
+      CASE (  5 )   ;   CALL hpg_isf    ( kt )      ! s-coordinate similar to sco modify for ice shelf
       END SELECT
       !
@@ -128,5 +132,5 @@
       !!
       NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     &
-         &                 ln_hpg_djc, ln_hpg_prj, ln_dynhpg_imp
+         &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, ln_dynhpg_imp
       !!----------------------------------------------------------------------
       !
@@ -148,4 +152,5 @@
          WRITE(numout,*) '      z-coord. - partial steps (interpolation)          ln_hpg_zps    = ', ln_hpg_zps
          WRITE(numout,*) '      s-coord. (standard jacobian formulation)          ln_hpg_sco    = ', ln_hpg_sco
+         WRITE(numout,*) '      s-coord. (standard jacobian formulation) for isf  ln_hpg_isf    = ', ln_hpg_isf
          WRITE(numout,*) '      s-coord. (Density Jacobian: Cubic polynomial)     ln_hpg_djc    = ', ln_hpg_djc
          WRITE(numout,*) '      s-coord. (Pressure Jacobian: Cubic polynomial)    ln_hpg_prj    = ', ln_hpg_prj
@@ -158,8 +163,13 @@
                            & either  ln_hpg_sco or  ln_hpg_prj instead')
       !
-      IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) )   &
+      IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )   &
          &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:&
                            & the standard jacobian formulation hpg_sco or &
                            & the pressure jacobian formulation hpg_prj')
+
+      IF(       ln_hpg_isf .AND. .NOT. ln_isfcav )   &
+         &   CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' )
+      IF( .NOT. ln_hpg_isf .AND.       ln_isfcav )   &
+         &   CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' )
       !
       !                               ! Set nhpg from ln_hpg_... flags
@@ -169,4 +179,5 @@
       IF( ln_hpg_djc )   nhpg = 3
       IF( ln_hpg_prj )   nhpg = 4
+      IF( ln_hpg_isf )   nhpg = 5
       !
       !                               ! Consistency check
@@ -177,7 +188,9 @@
       IF( ln_hpg_djc )   ioptio = ioptio + 1
       IF( ln_hpg_prj )   ioptio = ioptio + 1
+      IF( ln_hpg_isf )   ioptio = ioptio + 1
       IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
-      IF( (ln_hpg_zco .OR. ln_hpg_zps .OR. ln_hpg_djc .OR. ln_hpg_prj ) .AND. nn_isf .NE. 0 )   &
-          &  CALL ctl_stop( 'Only hpg_sco has been corrected to work with ice shelf cavity.' )
+      ! 
+      ! initialisation of ice load
+      riceload(:,:)=0.0
       !
    END SUBROUTINE dyn_hpg_init
@@ -345,5 +358,4 @@
    END SUBROUTINE hpg_zps
 
-
    SUBROUTINE hpg_sco( kt )
       !!---------------------------------------------------------------------
@@ -366,4 +378,91 @@
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
+      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
+      REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
+      !!----------------------------------------------------------------------
+      !
+      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used'
+      ENDIF
+
+      ! Local constant initialization
+      zcoef0 = - grav * 0.5_wp
+      ! To use density and not density anomaly
+      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume
+      ELSE                 ;     znad = 0._wp         ! Fixed volume
+      ENDIF
+
+      ! Surface value
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ! hydrostatic pressure gradient along s-surfaces
+            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   &
+               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
+            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   &
+               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
+            ! s-coordinate pressure gradient correction
+            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
+               &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
+            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
+               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
+            ! add to the general momentum trend
+            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
+            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
+         END DO
+      END DO
+
+      ! 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(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
+                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
+                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
+               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
+                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
+                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
+               ! s-coordinate pressure gradient correction
+               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
+                  &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj)
+               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
+                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj)
+               ! add to the general momentum trend
+               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
+               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
+            END DO
+         END DO
+      END DO
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
+      !
+   END SUBROUTINE hpg_sco
+
+   SUBROUTINE hpg_isf( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE hpg_sco  ***
+      !!
+      !! ** Method  :   s-coordinate case. Jacobian scheme.
+      !!      The now hydrostatic pressure gradient at a given level, jk,
+      !!      is computed by taking the vertical integral of the in-situ
+      !!      density gradient along the model level from the suface to that
+      !!      level. s-coordinates (ln_sco): a corrective term is added
+      !!      to the horizontal pressure gradient :
+      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ]
+      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ]
+      !!      add it to the general momentum trend (ua,va).
+      !!         ua = ua - 1/e1u * zhpi
+      !!         va = va - 1/e2v * zhpj
+      !!      iceload is added and partial cell case are added to the top and bottom
+      !!      
+      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
+      !!----------------------------------------------------------------------
+      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
@@ -379,5 +478,5 @@
      IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
+         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
@@ -565,11 +664,6 @@
 !==================================================================================
 
-# if defined key_vectopt_loop
-         jj = 1
-         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-# else
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
-# endif
             iku = mbku(ji,jj)
             ikv = mbkv(ji,jj)
@@ -598,7 +692,5 @@
                va(ji,jj,ikv)   =  va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap
             END IF
-# if ! defined key_vectopt_loop
-         END DO
-# endif
+         END DO
       END DO
      
@@ -610,5 +702,5 @@
       CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj)
       !
-   END SUBROUTINE hpg_sco
+   END SUBROUTINE hpg_isf
 
 
@@ -864,8 +956,10 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdept, zrhh
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zsshu_n, zsshv_n
       !!----------------------------------------------------------------------
       !
       CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
       CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh )
+      CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n )
       !
       IF( kt == nit000 ) THEN
@@ -948,8 +1042,19 @@
 
       ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1)
+
+      ! Prepare zsshu_n and zsshv_n
       DO jj = 2, jpjm1
         DO ji = 2, jpim1
-          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshu_n for ztilde compilation
-          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshv_n for ztilde compilation
+          zsshu_n(ji,jj) = (e12u(ji,jj) * sshn(ji,jj) + e12u(ji+1, jj) * sshn(ji+1,jj)) * &
+                         & r1_e12u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
+          zsshv_n(ji,jj) = (e12v(ji,jj) * sshn(ji,jj) + e12v(ji+1, jj) * sshn(ji,jj+1)) * &
+                         & r1_e12v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
+        END DO
+      END DO
+
+      DO jj = 2, jpjm1
+        DO ji = 2, jpim1
+          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - zsshu_n(ji,jj) * znad) 
+          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - zsshv_n(ji,jj) * znad)
         END DO
       END DO
@@ -1113,4 +1218,5 @@
       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
       CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh )
+      CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n )
       !
    END SUBROUTINE hpg_prj
