Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 5722)
+++ branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 5758)
@@ -93,16 +93,15 @@
       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad of prd at u- & v-pts (bottom)
       !
-      INTEGER  ::   ji, jj, jn      ! Dummy loop indices
-      INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points
-      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos
-      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             ! 
+      INTEGER  ::   ji, jj, jn                  ! Dummy loop indices
+      INTEGER  ::   iku, ikv, ikum1, ikvm1      ! partial step level (ocean bottom level) at u- and v-points
+      REAL(wp) ::   ze3wu, ze3wv, zmaxu, zmaxv  ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj)      ::   zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos
+      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::   zti, ztj             ! 
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde')
-      !
-      pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ;
-      zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ;
-      zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ;
+      IF( nn_timing == 1 )   CALL timing_start( 'zps_hde')
+      !
+      pgtu(:,:,:)=0._wp   ;   zti (:,:,:)=0._wp   ;   zhi (:,:  )=0._wp
+      pgtv(:,:,:)=0._wp   ;   ztj (:,:,:)=0._wp   ;   zhj (:,:  )=0._wp
       !
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
@@ -149,8 +148,8 @@
          !
       END DO
-
-      ! horizontal derivative of density anomalies (rd)
-      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
-         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ; 
+      !                
+      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
+         pgru(:,:) = 0._wp
+         pgrv(:,:) = 0._wp                ! depth of the partial step level
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
@@ -167,12 +166,9 @@
             END DO
          END DO
-
-         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
-         ! step and store it in  zri, zrj for each  case
-         CALL eos( zti, zhi, zri )  
-         CALL eos( ztj, zhj, zrj )
-
-         ! Gradient of density at the last level 
-         DO jj = 1, jpjm1
+         !
+         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj 
+         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj 
+         !
+         DO jj = 1, jpjm1                 ! Gradient of density at the last level 
             DO ji = 1, jpim1
                iku = mbku(ji,jj)
@@ -192,11 +188,12 @@
       END IF
       !
-      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde')
+      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde')
       !
    END SUBROUTINE zps_hde
-   !
-   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv,   &
-      &                          prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv,  &
-      &                   pgtui, pgtvi, pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi )
+
+
+   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu , pgtv , pgtui, pgtvi,                                   &
+      &                              prd, pgru , pgrv , pmru , pmrv , pgzu , pgzv , pge3ru , pge3rv ,   &
+      &                                   pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE zps_hde  ***
@@ -245,18 +242,20 @@
       !!              - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points 
       !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index
-      INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields
-      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts 
-      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui, pgtvi  ! hor. grad. of stra at u- & v-pts (ISF)
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv      ! hor. grad of prd at u- & v-pts (bottom)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru, pmrv      ! hor. sum  of prd at u- & v-pts (bottom)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu, pgzv      ! hor. grad of z   at u- & v-pts (bottom)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru, pge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi      ! hor. grad of prd at u- & v-pts (top)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui, pmrvi      ! hor. sum  of prd at u- & v-pts (top)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui, pgzvi      ! hor. grad of z   at u- & v-pts (top)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui, pge3rvi  ! hor. grad of prd weighted by local e3w at u- & v-pts (top)
+      INTEGER                              , INTENT(in   )           ::  kt                ! ocean time-step index
+      INTEGER                              , INTENT(in   )           ::  kjpt              ! number of tracers
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta               ! 4D tracers fields
+      !                                                              !!  u-point ! v-point !
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu    , pgtv    ! bottom GRADh( ptra )  
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui   , pgtvi   ! top    GRADh( ptra )
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd               ! 3D density anomaly fields
+      !                                                              !!  u-point ! v-point !
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru    , pgrv    ! bottom GRADh( prd  )
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru    , pmrv    ! bottom SUM  ( prd  )
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu    , pgzv    ! bottom GRADh( z    ) 
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru  , pge3rv  ! bottom GRADh( prd  ) weighted by e3w
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui   , pgrvi   ! top    GRADh( prd  ) 
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui   , pmrvi   ! top    SUM  ( prd  ) 
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui   , pgzvi   ! top    GRADh( z    ) 
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui , pge3rvi ! top    GRADh( prd  ) weighted by e3w
       !
       INTEGER  ::   ji, jj, jn      ! Dummy loop indices
@@ -269,8 +268,8 @@
       IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_isf')
       !
-      pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ;
-      pgtui(:,:,:)=0.0_wp ; pgtvi(:,:,:)=0.0_wp ;
-      zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ;
-      zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ;
+      pgtu (:,:,:) = 0._wp   ;   pgtv (:,:,:) =0._wp
+      pgtui(:,:,:) = 0._wp   ;   pgtvi(:,:,:) =0._wp
+      zti  (:,:,:) = 0._wp   ;   ztj  (:,:,:) =0._wp
+      zhi  (:,:  ) = 0._wp   ;   zhj  (:,:  ) =0._wp
       !
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
@@ -322,11 +321,12 @@
       END DO
 
-      ! horizontal derivative of density anomalies (rd)
-      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
-         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ; 
-         pgzu(:,:)=0.0_wp   ; pgzv(:,:)=0.0_wp ;
-         pmru(:,:)=0.0_wp   ; pmru(:,:)=0.0_wp ;
-         pge3ru(:,:)=0.0_wp ; pge3rv(:,:)=0.0_wp ;
-         DO jj = 1, jpjm1
+      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
+         !
+         pgru  (:,:)=0._wp   ;   pgrv  (:,:) = 0._wp
+         pgzu  (:,:)=0._wp   ;   pgzv  (:,:) = 0._wp 
+         pmru  (:,:)=0._wp   ;   pmru  (:,:) = 0._wp 
+         pge3ru(:,:)=0._wp   ;   pge3rv(:,:) = 0._wp 
+         !
+         DO jj = 1, jpjm1                 ! depth of the partial step level
             DO ji = 1, jpim1
                iku = mbku(ji,jj)
@@ -334,5 +334,5 @@
                ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
-
+               !
                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu     ! i-direction: case 1
                ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) + ze3wu    ! -     -      case 2
@@ -343,12 +343,9 @@
             END DO
          END DO
-         
-         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
-         ! step and store it in  zri, zrj for each  case
-         CALL eos( zti, zhi, zri )  
-         CALL eos( ztj, zhj, zrj )
-
-         ! Gradient of density at the last level 
-         DO jj = 1, jpjm1
+         !
+         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj 
+         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj 
+
+         DO jj = 1, jpjm1                 ! Gradient of density at the last level 
             DO ji = 1, jpim1
                iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
@@ -394,5 +391,7 @@
          !
       END IF
-         ! (ISH)  compute grui and gruvi
+      !
+      !     !==  (ISH)  compute grui and gruvi  ==!
+      !
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            !
          DO jj = 1, jpjm1
@@ -442,12 +441,12 @@
       END DO
 
-      ! horizontal derivative of density anomalies (rd)
-      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
+      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
+         !
          pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ;
          pgzui(:,:)  =0.0_wp ; pgzvi(:,:)  =0.0_wp ;
          pmrui(:,:)  =0.0_wp ; pmrui(:,:)  =0.0_wp ;
          pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ;
-
-         DO jj = 1, jpjm1
+         !
+         DO jj = 1, jpjm1        ! depth of the partial step level
             DO ji = 1, jpim1
                iku = miku(ji,jj)
@@ -455,5 +454,5 @@
                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))
                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( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1
                ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2
@@ -464,12 +463,9 @@
             END DO
          END DO
-
-         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
-         ! step and store it in  zri, zrj for each  case
-         CALL eos( zti, zhi, zri )  
-         CALL eos( ztj, zhj, zrj )
-
-         ! Gradient of density at the last level 
-         DO jj = 1, jpjm1
+         !
+         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj 
+         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj 
+         !
+         DO jj = 1, jpjm1                 ! Gradient of density at the last level 
             DO ji = 1, jpim1
                iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1
@@ -482,6 +478,6 @@
                  pmrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1 
                  pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                  &
-                                * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   &
-                                   - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1
+                    &           * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   &
+                    &              - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1
                ELSE
                  pgzui  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu)
@@ -489,6 +485,6 @@
                  pmrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2
                  pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                   &
-                                * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  &
-                                   -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2
+                    &           * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  &
+                    &              -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2
                ENDIF
                IF( ze3wv >= 0._wp ) THEN
@@ -497,6 +493,6 @@
                  pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1
                  pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                  & 
-                                * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  &
-                                   - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1
+                     &           * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  &
+                     &                - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1
                                   ! + 2 due to the formulation in density and not in anomalie in hpg sco
                ELSE
@@ -505,6 +501,6 @@
                  pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2
                  pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                   &
-                                * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) &
-                                   -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2
+                    &           * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) &
+                    &              -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2
                ENDIF
             END DO
@@ -517,5 +513,5 @@
       END IF  
       !
-      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde_isf')
+      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde_isf')
       !
    END SUBROUTINE zps_hde_isf
