New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 6808 for branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90 – NEMO

Ignore:
Timestamp:
2016-07-19T10:38:35+02:00 (8 years ago)
Author:
jamesharle
Message:

merge with trunk@6232 for consistency with SSB code

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r5120 r6808  
    3232 
    3333   !! * Substitutions 
    34 #  include "domzgr_substitute.h90" 
    3534#  include "vectopt_loop_substitute.h90" 
    3635   !!---------------------------------------------------------------------- 
     
    9392      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad of prd at u- & v-pts (bottom) 
    9493      ! 
    95       INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    96       INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points 
    97       REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars 
    98       REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos 
    99       REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             !  
    100       !!---------------------------------------------------------------------- 
    101       ! 
    102       IF( nn_timing == 1 )  CALL timing_start( 'zps_hde') 
    103       ! 
    104       pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 
    105       zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 
    106       zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ; 
     94      INTEGER  ::   ji, jj, jn                  ! Dummy loop indices 
     95      INTEGER  ::   iku, ikv, ikum1, ikvm1      ! partial step level (ocean bottom level) at u- and v-points 
     96      REAL(wp) ::   ze3wu, ze3wv, zmaxu, zmaxv  ! local scalars 
     97      REAL(wp), DIMENSION(jpi,jpj)      ::   zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos 
     98      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::   zti, ztj             !  
     99      !!---------------------------------------------------------------------- 
     100      ! 
     101      IF( nn_timing == 1 )   CALL timing_start( 'zps_hde') 
     102      ! 
     103      pgtu(:,:,:)=0._wp   ;   zti (:,:,:)=0._wp   ;   zhi (:,:  )=0._wp 
     104      pgtv(:,:,:)=0._wp   ;   ztj (:,:,:)=0._wp   ;   zhj (:,:  )=0._wp 
    107105      ! 
    108106      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
     
    112110               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    113111               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
    114                ze3wu = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    115                ze3wv = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     112!!gm BUG ? when applied to before fields, e3w_b should be used.... 
     113               ze3wu = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku) 
     114               ze3wv = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv) 
    116115               ! 
    117116               ! i- direction 
    118117               IF( ze3wu >= 0._wp ) THEN      ! case 1 
    119                   zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
     118                  zmaxu =  ze3wu / e3w_n(ji+1,jj,iku) 
    120119                  ! interpolated values of tracers 
    121120                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     
    123122                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    124123               ELSE                           ! case 2 
    125                   zmaxu = -ze3wu / fse3w(ji,jj,iku) 
     124                  zmaxu = -ze3wu / e3w_n(ji,jj,iku) 
    126125                  ! interpolated values of tracers 
    127126                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
     
    132131               ! j- direction 
    133132               IF( ze3wv >= 0._wp ) THEN      ! case 1 
    134                   zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
     133                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikv) 
    135134                  ! interpolated values of tracers 
    136135                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
     
    138137                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    139138               ELSE                           ! case 2 
    140                   zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
     139                  zmaxv =  -ze3wv / e3w_n(ji,jj,ikv) 
    141140                  ! interpolated values of tracers 
    142141                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
     
    149148         ! 
    150149      END DO 
    151  
    152       ! horizontal derivative of density anomalies (rd) 
    153       IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    154          pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ;  
     150      !                 
     151      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part) 
     152         pgru(:,:) = 0._wp 
     153         pgrv(:,:) = 0._wp                ! depth of the partial step level 
    155154         DO jj = 1, jpjm1 
    156155            DO ji = 1, jpim1 
    157156               iku = mbku(ji,jj) 
    158157               ikv = mbkv(ji,jj) 
    159                ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    160                ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
    161                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1 
    162                ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2 
    163                ENDIF 
    164                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1 
    165                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2 
    166                ENDIF 
    167             END DO 
    168          END DO 
    169  
    170          ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
    171          ! step and store it in  zri, zrj for each  case 
    172          CALL eos( zti, zhi, zri )   
    173          CALL eos( ztj, zhj, zrj ) 
    174  
    175          ! Gradient of density at the last level  
    176          DO jj = 1, jpjm1 
     158               ze3wu  = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku) 
     159               ze3wv  = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv) 
     160               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)     ! i-direction: case 1 
     161               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)     ! -     -      case 2 
     162               ENDIF 
     163               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)     ! j-direction: case 1 
     164               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)     ! -     -      case 2 
     165               ENDIF 
     166            END DO 
     167         END DO 
     168         ! 
     169         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj  
     170         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj  
     171         ! 
     172         DO jj = 1, jpjm1                 ! Gradient of density at the last level  
    177173            DO ji = 1, jpim1 
    178174               iku = mbku(ji,jj) 
    179175               ikv = mbkv(ji,jj) 
    180                ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    181                ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     176               ze3wu  = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku) 
     177               ze3wv  = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv) 
    182178               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1 
    183179               ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2 
     
    192188      END IF 
    193189      ! 
    194       IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde') 
     190      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde') 
    195191      ! 
    196192   END SUBROUTINE zps_hde 
    197193   ! 
    198    SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv,   & 
    199       &                          prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv,  & 
    200       &                   pgtui, pgtvi, pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi ) 
    201       !!---------------------------------------------------------------------- 
    202       !!                     ***  ROUTINE zps_hde  *** 
     194   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  & 
     195      &                          prd, pgru, pgrv, pgrui, pgrvi ) 
     196      !!---------------------------------------------------------------------- 
     197      !!                     ***  ROUTINE zps_hde_isf  *** 
    203198      !!                     
    204199      !! ** Purpose :   Compute the horizontal derivative of T, S and rho 
    205200      !!      at u- and v-points with a linear interpolation for z-coordinate 
    206       !!      with partial steps. 
     201      !!      with partial steps for top (ice shelf) and bottom. 
    207202      !! 
    208203      !! ** Method  :   In z-coord with partial steps, scale factors on last  
    209204      !!      levels are different for each grid point, so that T, S and rd  
    210205      !!      points are not at the same depth as in z-coord. To have horizontal 
    211       !!      gradients again, we interpolate T and S at the good depth :  
     206      !!      gradients again, we interpolate T and S at the good depth : 
     207      !!      For the bottom case: 
    212208      !!      Linear interpolation of T, S    
    213209      !!         Computation of di(tb) and dj(tb) by vertical interpolation: 
     
    238234      !!          di(rho) = rd~ - rd(i,j,k)   or   rd(i+1,j,k) - rd~ 
    239235      !! 
     236      !!      For the top case (ice shelf): As for the bottom case but upside down 
     237      !! 
    240238      !! ** Action  : compute for top and bottom interfaces 
    241239      !!              - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points 
    242240      !!              - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points 
    243       !!              - pmru, pmrv, pmrui, pmrvi: horizontal sum of rho at u- & v- point (used in dynhpg with vvl) 
    244       !!              - pgzu, pgzv, pgzui, pgzvi: horizontal gradient of z at u- and v- point (used in dynhpg with vvl) 
    245       !!              - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points  
    246       !!---------------------------------------------------------------------- 
    247       INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index 
    248       INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers 
    249       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields 
    250       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts  
    251       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui, pgtvi  ! hor. grad. of stra at u- & v-pts (ISF) 
    252       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
    253       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv      ! hor. grad of prd at u- & v-pts (bottom) 
    254       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru, pmrv      ! hor. sum  of prd at u- & v-pts (bottom) 
    255       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu, pgzv      ! hor. grad of z   at u- & v-pts (bottom) 
    256       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru, pge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom) 
    257       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi      ! hor. grad of prd at u- & v-pts (top) 
    258       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui, pmrvi      ! hor. sum  of prd at u- & v-pts (top) 
    259       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui, pgzvi      ! hor. grad of z   at u- & v-pts (top) 
    260       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui, pge3rvi  ! hor. grad of prd weighted by local e3w at u- & v-pts (top) 
     241      !!---------------------------------------------------------------------- 
     242      INTEGER                              , INTENT(in   )           ::  kt           ! ocean time-step index 
     243      INTEGER                              , INTENT(in   )           ::  kjpt         ! number of tracers 
     244      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta          ! 4D tracers fields 
     245      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv   ! hor. grad. of ptra at u- & v-pts  
     246      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 
     247      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd          ! 3D density anomaly fields 
     248      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv   ! hor. grad of prd at u- & v-pts (bottom) 
     249      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) 
    261250      ! 
    262251      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    263252      INTEGER  ::   iku, ikv, ikum1, ikvm1,ikup1, ikvp1   ! partial step level (ocean bottom level) at u- and v-points 
    264       REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv, zdzwu, zdzwv, zdzwuip1, zdzwvjp1  ! temporary scalars 
     253      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv             ! temporary scalars 
    265254      REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos 
    266255      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             !  
     
    269258      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_isf') 
    270259      ! 
    271       pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 
    272       pgtui(:,:,:)=0.0_wp ; pgtvi(:,:,:)=0.0_wp ; 
    273       zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 
    274       zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ; 
     260      pgtu (:,:,:) = 0._wp   ;   pgtv (:,:,:) =0._wp 
     261      pgtui(:,:,:) = 0._wp   ;   pgtvi(:,:,:) =0._wp 
     262      zti  (:,:,:) = 0._wp   ;   ztj  (:,:,:) =0._wp 
     263      zhi  (:,:  ) = 0._wp   ;   zhj  (:,:  ) =0._wp 
    275264      ! 
    276265      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
     
    278267         DO jj = 1, jpjm1 
    279268            DO ji = 1, jpim1 
    280                iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    281                ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
     269 
     270               iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
     271               ikv = mbkv(ji,jj); ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
     272               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku) 
     273               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv) 
     274               ! 
     275               ! i- direction 
     276               IF( ze3wu >= 0._wp ) THEN      ! case 1 
     277                  zmaxu =  ze3wu / e3w_n(ji+1,jj,iku) 
     278                  ! interpolated values of tracers 
     279                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     280                  ! gradient of  tracers 
     281                  pgtu(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     282               ELSE                           ! case 2 
     283                  zmaxu = -ze3wu / e3w_n(ji,jj,iku) 
     284                  ! interpolated values of tracers 
     285                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
     286                  ! gradient of tracers 
     287                  pgtu(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     288               ENDIF 
     289               ! 
     290               ! j- direction 
     291               IF( ze3wv >= 0._wp ) THEN      ! case 1 
     292                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikv) 
     293                  ! interpolated values of tracers 
     294                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
     295                  ! gradient of tracers 
     296                  pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     297               ELSE                           ! case 2 
     298                  zmaxv =  -ze3wv / e3w_n(ji,jj,ikv) 
     299                  ! interpolated values of tracers 
     300                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
     301                  ! gradient of tracers 
     302                  pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     303               ENDIF 
     304 
     305            END DO 
     306         END DO 
     307         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     308         ! 
     309      END DO 
     310 
     311      ! horizontal derivative of density anomalies (rd) 
     312      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
     313         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ;  
     314         ! 
     315         DO jj = 1, jpjm1 
     316            DO ji = 1, jpim1 
     317 
     318               iku = mbku(ji,jj) 
     319               ikv = mbkv(ji,jj) 
     320               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku) 
     321               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv) 
     322               ! 
     323               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)    ! i-direction: case 1 
     324               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)    ! -     -      case 2 
     325               ENDIF 
     326               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)    ! j-direction: case 1 
     327               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)    ! -     -      case 2 
     328               ENDIF 
     329 
     330            END DO 
     331         END DO 
     332 
     333         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
     334         ! step and store it in  zri, zrj for each  case 
     335         CALL eos( zti, zhi, zri ) 
     336         CALL eos( ztj, zhj, zrj ) 
     337 
     338         DO jj = 1, jpjm1                 ! Gradient of density at the last level  
     339            DO ji = 1, jpim1 
     340               iku = mbku(ji,jj) 
     341               ikv = mbkv(ji,jj) 
     342               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku) 
     343               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv) 
     344 
     345               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1 
     346               ELSE                        ;   pgru(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2 
     347               ENDIF 
     348               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = ssvmask(ji,jj) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1 
     349               ELSE                        ;   pgrv(ji,jj) = ssvmask(ji,jj) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2 
     350               ENDIF 
     351 
     352            END DO 
     353         END DO 
     354 
     355         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions 
     356         ! 
     357      END IF 
     358      ! 
     359      !     !==  (ISH)  compute grui and gruvi  ==! 
     360      ! 
     361      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            ! 
     362         DO jj = 1, jpjm1 
     363            DO ji = 1, jpim1 
     364               iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1 
     365               ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1 
     366               ! 
    282367               ! (ISF) case partial step top and bottom in adjacent cell in vertical 
    283368               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 
    284369               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 
    285370               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 
    286                ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
    287                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    288                ! 
     371               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku) 
     372               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv)  
     373 
    289374               ! i- direction 
    290375               IF( ze3wu >= 0._wp ) THEN      ! case 1 
    291                   zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
    292                   ! interpolated values of tracers 
    293                   zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     376                  zmaxu = ze3wu / e3w_n(ji+1,jj,ikup1) 
     377                  ! interpolated values of tracers 
     378                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikup1,jn) - pta(ji+1,jj,iku,jn) ) 
     379                  ! gradient of tracers 
     380                  pgtui(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     381               ELSE                           ! case 2 
     382                  zmaxu = - ze3wu / e3w_n(ji,jj,ikup1) 
     383                  ! interpolated values of tracers 
     384                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikup1,jn) - pta(ji,jj,iku,jn) ) 
    294385                  ! gradient of  tracers 
    295                   pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    296                ELSE                           ! case 2 
    297                   zmaxu = -ze3wu / fse3w(ji,jj,iku) 
    298                   ! interpolated values of tracers 
    299                   zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
    300                   ! gradient of tracers 
    301                   pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     386                  pgtui(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
    302387               ENDIF 
    303388               ! 
    304389               ! j- direction 
    305390               IF( ze3wv >= 0._wp ) THEN      ! case 1 
    306                   zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
    307                   ! interpolated values of tracers 
    308                   ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
    309                   ! gradient of tracers 
    310                   pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    311                ELSE                           ! case 2 
    312                   zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
    313                   ! interpolated values of tracers 
    314                   ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
    315                   ! gradient of tracers 
    316                   pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
    317                ENDIF 
    318             END DO 
    319          END DO 
    320          CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     391                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikvp1) 
     392                  ! interpolated values of tracers 
     393                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvp1,jn) - pta(ji,jj+1,ikv,jn) ) 
     394                  ! gradient of tracers 
     395                  pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     396               ELSE                           ! case 2 
     397                  zmaxv =  - ze3wv / e3w_n(ji,jj,ikvp1) 
     398                  ! interpolated values of tracers 
     399                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvp1,jn) - pta(ji,jj,ikv,jn) ) 
     400                  ! gradient of tracers 
     401                  pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     402               ENDIF 
     403 
     404            END DO 
     405         END DO 
     406         CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. ); CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
    321407         ! 
    322408      END DO 
    323409 
    324       ! horizontal derivative of density anomalies (rd) 
    325       IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    326          pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ;  
    327          pgzu(:,:)=0.0_wp   ; pgzv(:,:)=0.0_wp ; 
    328          pmru(:,:)=0.0_wp   ; pmru(:,:)=0.0_wp ; 
    329          pge3ru(:,:)=0.0_wp ; pge3rv(:,:)=0.0_wp ; 
    330          DO jj = 1, jpjm1 
    331             DO ji = 1, jpim1 
    332                iku = mbku(ji,jj) 
    333                ikv = mbkv(ji,jj) 
    334                ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
    335                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    336  
    337                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu     ! i-direction: case 1 
    338                ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) + ze3wu    ! -     -      case 2 
    339                ENDIF 
    340                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) - ze3wv    ! j-direction: case 1 
    341                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) + ze3wv    ! -     -      case 2 
    342                ENDIF 
    343             END DO 
    344          END DO 
    345           
    346          ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
    347          ! step and store it in  zri, zrj for each  case 
    348          CALL eos( zti, zhi, zri )   
    349          CALL eos( ztj, zhj, zrj ) 
    350  
    351          ! Gradient of density at the last level  
    352          DO jj = 1, jpjm1 
    353             DO ji = 1, jpim1 
    354                iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    355                ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    356                ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
    357                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    358                IF( ze3wu >= 0._wp ) THEN  
    359                   pgzu(ji,jj) = (fsde3w(ji+1,jj,iku) - ze3wu) - fsde3w(ji,jj,iku) 
    360                   pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1 
    361                   pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) + prd(ji,jj,iku) )   ! i: 1  
    362                   pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  & 
    363                                 * ( (fse3w(ji+1,jj,iku) - ze3wu )* ( zri(ji  ,jj    ) + prd(ji+1,jj,ikum1) + 2._wp) & 
    364                                    - fse3w(ji  ,jj,iku)          * ( prd(ji  ,jj,iku) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2 
    365                ELSE   
    366                   pgzu(ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) + ze3wu) 
    367                   pgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2 
    368                   pmru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) )   ! i: 2 
    369                   pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  & 
    370                                 * (  fse3w(ji+1,jj,iku)          * ( prd(ji+1,jj,iku) + prd(ji+1,jj,ikum1) + 2._wp) & 
    371                                    -(fse3w(ji  ,jj,iku) + ze3wu) * ( zri(ji  ,jj    ) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2 
    372                ENDIF 
    373                IF( ze3wv >= 0._wp ) THEN 
    374                   pgzv(ji,jj) = (fsde3w(ji,jj+1,ikv) - ze3wv) - fsde3w(ji,jj,ikv)  
    375                   pgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1 
    376                   pmrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )   ! j: 1 
    377                   pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  & 
    378                                 * ( (fse3w(ji,jj+1,ikv) - ze3wv )* ( zrj(ji,jj      ) + prd(ji,jj+1,ikvm1) + 2._wp) & 
    379                                    - fse3w(ji,jj  ,ikv)          * ( prd(ji,jj  ,ikv) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2 
    380                ELSE  
    381                   pgzv(ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) + ze3wv) 
    382                   pgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2 
    383                   pmrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )   ! j: 2 
    384                   pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  & 
    385                                 * (  fse3w(ji,jj+1,ikv)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikvm1) + 2._wp) & 
    386                                    -(fse3w(ji,jj  ,ikv) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2 
    387                ENDIF 
    388             END DO 
    389          END DO 
    390          CALL lbc_lnk( pgru   , 'U', -1. )   ;   CALL lbc_lnk( pgrv   , 'V', -1. )   ! Lateral boundary conditions 
    391          CALL lbc_lnk( pmru   , 'U',  1. )   ;   CALL lbc_lnk( pmrv   , 'V',  1. )   ! Lateral boundary conditions 
    392          CALL lbc_lnk( pgzu   , 'U', -1. )   ;   CALL lbc_lnk( pgzv   , 'V', -1. )   ! Lateral boundary conditions 
    393          CALL lbc_lnk( pge3ru , 'U', -1. )   ;   CALL lbc_lnk( pge3rv , 'V', -1. )   ! Lateral boundary conditions 
    394          ! 
    395       END IF 
    396          ! (ISH)  compute grui and gruvi 
    397       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            ! 
    398          DO jj = 1, jpjm1 
    399             DO ji = 1, jpim1 
    400                iku = miku(ji,jj)   ;  ikup1 = miku(ji,jj) + 1 
    401                ikv = mikv(ji,jj)   ;  ikvp1 = mikv(ji,jj) + 1 
    402                ! 
    403                ! (ISF) case partial step top and bottom in adjacent cell in vertical 
    404                ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 
    405                ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 
    406                ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 
    407                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))  
    408                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)) 
    409                ! i- direction 
    410                IF( ze3wu >= 0._wp ) THEN      ! case 1 
    411                   zmaxu = ze3wu / fse3w(ji+1,jj,iku+1) 
    412                   ! interpolated values of tracers 
    413                   zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 
    414                   ! gradient of tracers 
    415                   pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    416                ELSE                           ! case 2 
    417                   zmaxu = - ze3wu / fse3w(ji,jj,iku+1) 
    418                   ! interpolated values of tracers 
    419                   zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 
    420                   ! gradient of  tracers 
    421                   pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
    422                ENDIF 
    423                ! 
    424                ! j- direction 
    425                IF( ze3wv >= 0._wp ) THEN      ! case 1 
    426                   zmaxv =  ze3wv / fse3w(ji,jj+1,ikv+1) 
    427                   ! interpolated values of tracers 
    428                   ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 
    429                   ! gradient of tracers 
    430                   pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    431                ELSE                           ! case 2 
    432                   zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1) 
    433                   ! interpolated values of tracers 
    434                   ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 
    435                   ! gradient of tracers 
    436                   pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
    437                ENDIF 
    438             END DO!! 
    439          END DO!! 
    440          CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
    441          ! 
    442       END DO 
    443  
    444       ! horizontal derivative of density anomalies (rd) 
    445       IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    446          pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ; 
    447          pgzui(:,:)  =0.0_wp ; pgzvi(:,:)  =0.0_wp ; 
    448          pmrui(:,:)  =0.0_wp ; pmrui(:,:)  =0.0_wp ; 
    449          pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ; 
    450  
    451          DO jj = 1, jpjm1 
    452             DO ji = 1, jpim1 
     410      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part) 
     411         ! 
     412         pgrui(:,:)  =0.0_wp; pgrvi(:,:)  =0.0_wp; 
     413         DO jj = 1, jpjm1 
     414            DO ji = 1, jpim1 
     415 
    453416               iku = miku(ji,jj) 
    454417               ikv = mikv(ji,jj) 
    455                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)) 
    456                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)) 
    457  
    458                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1 
    459                ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2 
    460                ENDIF 
    461                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1 
    462                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2 
    463                ENDIF 
    464             END DO 
    465          END DO 
    466  
    467          ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
    468          ! step and store it in  zri, zrj for each  case 
    469          CALL eos( zti, zhi, zri )   
    470          CALL eos( ztj, zhj, zrj ) 
    471  
    472          ! Gradient of density at the last level  
    473          DO jj = 1, jpjm1 
    474             DO ji = 1, jpim1 
    475                iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1 
    476                ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1 
    477                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)) 
    478                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)) 
    479                IF( ze3wu >= 0._wp ) THEN 
    480                  pgzui  (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 
    481                  pgrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1 
    482                  pmrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1  
    483                  pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                  & 
    484                                 * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   & 
    485                                    - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1 
    486                ELSE 
    487                  pgzui  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 
    488                  pgrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2 
    489                  pmrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2 
    490                  pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                   & 
    491                                 * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  & 
    492                                    -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2 
    493                ENDIF 
    494                IF( ze3wv >= 0._wp ) THEN 
    495                  pgzvi  (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv)  
    496                  pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1 
    497                  pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1 
    498                  pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                  &  
    499                                 * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  & 
    500                                    - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1 
    501                                   ! + 2 due to the formulation in density and not in anomalie in hpg sco 
    502                ELSE 
    503                  pgzvi  (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 
    504                  pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2 
    505                  pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2 
    506                  pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                   & 
    507                                 * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 
    508                                    -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2 
    509                ENDIF 
    510             END DO 
    511          END DO 
    512          CALL lbc_lnk( pgrui   , 'U', -1. )   ;   CALL lbc_lnk( pgrvi   , 'V', -1. )   ! Lateral boundary conditions 
    513          CALL lbc_lnk( pmrui   , 'U',  1. )   ;   CALL lbc_lnk( pmrvi   , 'V',  1. )   ! Lateral boundary conditions 
    514          CALL lbc_lnk( pgzui   , 'U', -1. )   ;   CALL lbc_lnk( pgzvi   , 'V', -1. )   ! Lateral boundary conditions 
    515          CALL lbc_lnk( pge3rui , 'U', -1. )   ;   CALL lbc_lnk( pge3rvi , 'V', -1. )   ! Lateral boundary conditions 
     418               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku) 
     419               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv)  
     420               ! 
     421               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)    ! i-direction: case 1 
     422               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)    ! -     -      case 2 
     423               ENDIF 
     424 
     425               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)    ! j-direction: case 1 
     426               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)    ! -     -      case 2 
     427               ENDIF 
     428 
     429            END DO 
     430         END DO 
     431         ! 
     432         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj  
     433         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj  
     434         ! 
     435         DO jj = 1, jpjm1                 ! Gradient of density at the last level  
     436            DO ji = 1, jpim1 
     437               iku = miku(ji,jj)  
     438               ikv = mikv(ji,jj)  
     439               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku) 
     440               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv)  
     441 
     442               IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) ) ! i: 1 
     443               ELSE                      ; pgrui(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj  ,iku) - zri(ji,jj    ) ) ! i: 2 
     444               ENDIF 
     445               IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( zrj(ji  ,jj      ) - prd(ji,jj,ikv) ) ! j: 1 
     446               ELSE                      ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( prd(ji  ,jj+1,ikv) - zrj(ji,jj    ) ) ! j: 2 
     447               ENDIF 
     448 
     449            END DO 
     450         END DO 
     451         CALL lbc_lnk( pgrui   , 'U', -1. ); CALL lbc_lnk( pgrvi   , 'V', -1. )   ! Lateral boundary conditions 
    516452         ! 
    517453      END IF   
    518454      ! 
    519       IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde_isf') 
     455      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde_isf') 
    520456      ! 
    521457   END SUBROUTINE zps_hde_isf 
Note: See TracChangeset for help on using the changeset viewer.