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 10928 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynhpg.F90 – NEMO

Ignore:
Timestamp:
2019-05-03T17:44:56+02:00 (5 years ago)
Author:
davestorkey
Message:

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Finish converting DYN module, including some updates to previously processed modules, but excluding dynnxt.F90 (which needs to be completely rewritten) and wet_dry.F90 - I need to talk to Enda about this.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynhpg.F90

    r10491 r10928  
    8181CONTAINS 
    8282 
    83    SUBROUTINE dyn_hpg( kt ) 
     83   SUBROUTINE dyn_hpg( kt, Kmm, puu, pvv, Krhs ) 
    8484      !!--------------------------------------------------------------------- 
    8585      !!                  ***  ROUTINE dyn_hpg  *** 
     
    8888      !!              using the scheme defined in the namelist 
    8989      !! 
    90       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
     90      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 
    9191      !!             - send trends to trd_dyn for futher diagnostics (l_trddyn=T) 
    9292      !!---------------------------------------------------------------------- 
    93       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     93      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index 
     94      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     95      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
     96      ! 
    9497      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
    9598      !!---------------------------------------------------------------------- 
     
    97100      IF( ln_timing )   CALL timing_start('dyn_hpg') 
    98101      ! 
    99       IF( l_trddyn ) THEN                    ! Temporary saving of ua and va trends (l_trddyn) 
     102      IF( l_trddyn ) THEN                    ! Temporary saving of puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends (l_trddyn) 
    100103         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
    101          ztrdu(:,:,:) = ua(:,:,:) 
    102          ztrdv(:,:,:) = va(:,:,:) 
     104         ztrdu(:,:,:) = puu(:,:,:,Krhs) 
     105         ztrdv(:,:,:) = pvv(:,:,:,Krhs) 
    103106      ENDIF 
    104107      ! 
    105108      SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation 
    106       CASE ( np_zco )   ;   CALL hpg_zco    ( kt )      ! z-coordinate 
    107       CASE ( np_zps )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation) 
    108       CASE ( np_sco )   ;   CALL hpg_sco    ( kt )      ! s-coordinate (standard jacobian formulation) 
    109       CASE ( np_djc )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial) 
    110       CASE ( np_prj )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme) 
    111       CASE ( np_isf )   ;   CALL hpg_isf    ( kt )      ! s-coordinate similar to sco modify for ice shelf 
     109      CASE ( np_zco )   ;   CALL hpg_zco    ( kt, Kmm, puu, pvv, Krhs )  ! z-coordinate 
     110      CASE ( np_zps )   ;   CALL hpg_zps    ( kt, Kmm, puu, pvv, Krhs )  ! z-coordinate plus partial steps (interpolation) 
     111      CASE ( np_sco )   ;   CALL hpg_sco    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate (standard jacobian formulation) 
     112      CASE ( np_djc )   ;   CALL hpg_djc    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate (Density Jacobian with Cubic polynomial) 
     113      CASE ( np_prj )   ;   CALL hpg_prj    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate (Pressure Jacobian scheme) 
     114      CASE ( np_isf )   ;   CALL hpg_isf    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate similar to sco modify for ice shelf 
    112115      END SELECT 
    113116      ! 
    114117      IF( l_trddyn ) THEN      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics 
    115          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    116          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     118         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 
     119         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 
    117120         CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt ) 
    118121         DEALLOCATE( ztrdu , ztrdv ) 
    119122      ENDIF 
    120123      ! 
    121       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask,   & 
    122          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     124      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' hpg  - Ua: ', mask1=umask,   & 
     125         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    123126      ! 
    124127      IF( ln_timing )   CALL timing_stop('dyn_hpg') 
     
    127130 
    128131 
    129    SUBROUTINE dyn_hpg_init 
     132   SUBROUTINE dyn_hpg_init( Kmm ) 
    130133      !!---------------------------------------------------------------------- 
    131134      !!                 ***  ROUTINE dyn_hpg_init  *** 
     
    137140      !!      with the type of vertical coordinate used (zco, zps, sco) 
    138141      !!---------------------------------------------------------------------- 
     142      INTEGER, INTENT( in ) :: Kmm   ! ocean time level index 
     143      ! 
    139144      INTEGER ::   ioptio = 0      ! temporary integer 
    140145      INTEGER ::   ios             ! Local integer output status for namelist read 
     
    228233         ! 
    229234         DO jk = 1, jpk                   !- compute density of the water displaced by the ice shelf  
    230             CALL eos( zts_top(:,:,:), gdept_n(:,:,jk), zrhd(:,:,jk) ) 
     235            CALL eos( zts_top(:,:,:), gdept(:,:,jk,Kmm), zrhd(:,:,jk) ) 
    231236         END DO 
    232237         ! 
     
    239244            DO ji = 1, jpi                      ! divided by 2 later 
    240245               ikt = mikt(ji,jj) 
    241                ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 
     246               ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w(ji,jj,1,Kmm) * (1._wp - tmask(ji,jj,1)) 
    242247               DO jk = 2, ikt-1 
    243                   ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk) & 
     248                  ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w(ji,jj,jk,Kmm) & 
    244249                     &                              * (1._wp - tmask(ji,jj,jk)) 
    245250               END DO 
    246251               IF (ikt  >=  2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 
    247                   &                                              * ( risfdep(ji,jj) - gdept_n(ji,jj,ikt-1) ) 
     252                  &                                              * ( risfdep(ji,jj) - gdept(ji,jj,ikt-1,Kmm) ) 
    248253            END DO 
    249254         END DO 
     
    256261 
    257262 
    258    SUBROUTINE hpg_zco( kt ) 
     263   SUBROUTINE hpg_zco( kt, Kmm, puu, pvv, Krhs ) 
    259264      !!--------------------------------------------------------------------- 
    260265      !!                  ***  ROUTINE hpg_zco  *** 
     
    266271      !!      level:    zhpi = grav ..... 
    267272      !!                zhpj = grav ..... 
    268       !!      add it to the general momentum trend (ua,va). 
    269       !!            ua = ua - 1/e1u * zhpi 
    270       !!            va = va - 1/e2v * zhpj 
    271       !! 
    272       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
    273       !!---------------------------------------------------------------------- 
    274       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     273      !!      add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 
     274      !!            puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 
     275      !!            pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 
     276      !! 
     277      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 
     278      !!---------------------------------------------------------------------- 
     279      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index 
     280      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     281      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    275282      ! 
    276283      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
     
    290297      DO jj = 2, jpjm1 
    291298         DO ji = fs_2, fs_jpim1   ! vector opt. 
    292             zcoef1 = zcoef0 * e3w_n(ji,jj,1) 
     299            zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 
    293300            ! hydrostatic pressure gradient 
    294301            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 
    295302            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 
    296303            ! add to the general momentum trend 
    297             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
    298             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 
     304            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 
     305            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 
    299306         END DO 
    300307      END DO 
     
    305312         DO jj = 2, jpjm1 
    306313            DO ji = fs_2, fs_jpim1   ! vector opt. 
    307                zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 
     314               zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 
    308315               ! hydrostatic pressure gradient 
    309316               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   & 
     
    315322                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj) 
    316323               ! add to the general momentum trend 
    317                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
    318                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 
     324               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 
     325               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 
    319326            END DO 
    320327         END DO 
     
    324331 
    325332 
    326    SUBROUTINE hpg_zps( kt ) 
     333   SUBROUTINE hpg_zps( kt, Kmm, puu, pvv, Krhs ) 
    327334      !!--------------------------------------------------------------------- 
    328335      !!                 ***  ROUTINE hpg_zps  *** 
     
    330337      !! ** Method  :   z-coordinate plus partial steps case.  blahblah... 
    331338      !! 
    332       !! ** Action  : - Update (ua,va) with the now hydrastatic pressure trend 
    333       !!---------------------------------------------------------------------- 
    334       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     339      !! ** Action  : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 
     340      !!---------------------------------------------------------------------- 
     341      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index 
     342      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     343      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    335344      !! 
    336345      INTEGER  ::   ji, jj, jk                       ! dummy loop indices 
     
    347356 
    348357      ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level 
    349 !jc      CALL zps_hde    ( kt, jpts, tsn, gtsu, gtsv, rhd, gru , grv ) 
     358!jc      CALL zps_hde    ( kt, jpts, ts(:,:,:,:,Kmm), gtsu, gtsv, rhd, gru , grv ) 
    350359 
    351360      ! Local constant initialization 
     
    355364      DO jj = 2, jpjm1 
    356365         DO ji = fs_2, fs_jpim1   ! vector opt. 
    357             zcoef1 = zcoef0 * e3w_n(ji,jj,1) 
     366            zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 
    358367            ! hydrostatic pressure gradient 
    359368            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 
    360369            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 
    361370            ! add to the general momentum trend 
    362             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
    363             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 
     371            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 
     372            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 
    364373         END DO 
    365374      END DO 
     
    369378         DO jj = 2, jpjm1 
    370379            DO ji = fs_2, fs_jpim1   ! vector opt. 
    371                zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 
     380               zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 
    372381               ! hydrostatic pressure gradient 
    373382               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   & 
     
    379388                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj) 
    380389               ! add to the general momentum trend 
    381                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
    382                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 
     390               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 
     391               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 
    383392            END DO 
    384393         END DO 
     
    390399            iku = mbku(ji,jj) 
    391400            ikv = mbkv(ji,jj) 
    392             zcoef2 = zcoef0 * MIN( e3w_n(ji,jj,iku), e3w_n(ji+1,jj  ,iku) ) 
    393             zcoef3 = zcoef0 * MIN( e3w_n(ji,jj,ikv), e3w_n(ji  ,jj+1,ikv) ) 
     401            zcoef2 = zcoef0 * MIN( e3w(ji,jj,iku,Kmm), e3w(ji+1,jj  ,iku,Kmm) ) 
     402            zcoef3 = zcoef0 * MIN( e3w(ji,jj,ikv,Kmm), e3w(ji  ,jj+1,ikv,Kmm) ) 
    394403            IF( iku > 1 ) THEN            ! on i-direction (level 2 or more) 
    395                ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value 
     404               puu  (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) - zhpi(ji,jj,iku)         ! subtract old value 
    396405               zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one 
    397406                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) 
    398                ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend 
     407               puu  (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend 
    399408            ENDIF 
    400409            IF( ikv > 1 ) THEN            ! on j-direction (level 2 or more) 
    401                va  (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)         ! subtract old value 
     410               pvv  (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) - zhpj(ji,jj,ikv)         ! subtract old value 
    402411               zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one 
    403412                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) 
    404                va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend 
     413               pvv  (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend 
    405414            ENDIF 
    406415         END DO 
     
    410419 
    411420 
    412    SUBROUTINE hpg_sco( kt ) 
     421   SUBROUTINE hpg_sco( kt, Kmm, puu, pvv, Krhs ) 
    413422      !!--------------------------------------------------------------------- 
    414423      !!                  ***  ROUTINE hpg_sco  *** 
     
    422431      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ] 
    423432      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ] 
    424       !!      add it to the general momentum trend (ua,va). 
    425       !!         ua = ua - 1/e1u * zhpi 
    426       !!         va = va - 1/e2v * zhpj 
    427       !! 
    428       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
    429       !!---------------------------------------------------------------------- 
    430       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     433      !!      add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 
     434      !!         puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 
     435      !!         pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 
     436      !! 
     437      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 
     438      !!---------------------------------------------------------------------- 
     439      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index 
     440      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     441      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    431442      !! 
    432443      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices 
     
    453464        DO jj = 2, jpjm1 
    454465           DO ji = 2, jpim1  
    455              ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                & 
     466             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)               ,  ssh(ji+1,jj,Kmm) ) >                & 
    456467                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
    457                   &    MAX(  sshn(ji,jj) +  ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     468                  &    MAX(  ssh(ji,jj,Kmm) +  ht_0(ji,jj),  ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) )  & 
    458469                  &                                                       > rn_wdmin1 + rn_wdmin2 
    459              ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (       & 
    460                   &    MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
     470             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)              -  ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. (       & 
     471                  &    MAX(   ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                & 
    461472                  &    MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    462473 
     
    464475               zcpx(ji,jj) = 1.0_wp 
    465476             ELSE IF(ll_tmp2) THEN 
    466                ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    467                zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
    468                            &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) ) 
     477               ! no worries about  ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm) = 0, it won't happen ! here 
     478               zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
     479                           &    / (ssh(ji+1,jj,Kmm) - ssh(ji  ,jj,Kmm)) ) 
    469480             ELSE 
    470481               zcpx(ji,jj) = 0._wp 
    471482             END IF 
    472483       
    473              ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
     484             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                & 
    474485                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            & 
    475                   &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
     486                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) )  & 
    476487                  &                                                      > rn_wdmin1 + rn_wdmin2 
    477              ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (        & 
    478                   &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                & 
     488             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. (        & 
     489                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji,jj+1,Kmm) ) >                & 
    479490                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    480491 
     
    482493               zcpy(ji,jj) = 1.0_wp 
    483494             ELSE IF(ll_tmp2) THEN 
    484                ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    485                zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
    486                            &    / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
     495               ! no worries about  ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm) = 0, it won't happen ! here 
     496               zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
     497                           &    / (ssh(ji,jj+1,Kmm) - ssh(ji,jj  ,Kmm)) ) 
    487498             ELSE 
    488499               zcpy(ji,jj) = 0._wp 
     
    497508         DO ji = fs_2, fs_jpim1   ! vector opt. 
    498509            ! hydrostatic pressure gradient along s-surfaces 
    499             zhpi(ji,jj,1) = zcoef0 * (  e3w_n(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )    & 
    500                &                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e1u(ji,jj) 
    501             zhpj(ji,jj,1) = zcoef0 * (  e3w_n(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )    & 
    502                &                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e2v(ji,jj) 
     510            zhpi(ji,jj,1) = zcoef0 * (  e3w(ji+1,jj  ,1,Kmm) * ( znad + rhd(ji+1,jj  ,1) )    & 
     511               &                      - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e1u(ji,jj) 
     512            zhpj(ji,jj,1) = zcoef0 * (  e3w(ji  ,jj+1,1,Kmm) * ( znad + rhd(ji  ,jj+1,1) )    & 
     513               &                      - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e2v(ji,jj) 
    503514            ! s-coordinate pressure gradient correction 
    504515            zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
    505                &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 
     516               &           * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 
    506517            zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
    507                &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 
     518               &           * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 
    508519            ! 
    509520            IF( ln_wd_il ) THEN 
     
    515526            ! 
    516527            ! add to the general momentum trend 
    517             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
    518             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
     528            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + zuap 
     529            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + zvap 
    519530         END DO 
    520531      END DO 
     
    526537               ! hydrostatic pressure gradient along s-surfaces 
    527538               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   & 
    528                   &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
    529                   &              - e3w_n(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  ) 
     539                  &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
     540                  &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  ) 
    530541               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   & 
    531                   &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
    532                   &              - e3w_n(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
     542                  &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
     543                  &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
    533544               ! s-coordinate pressure gradient correction 
    534545               zuap = -zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   & 
    535                   &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) 
     546                  &           * ( gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk) ) * r1_e1u(ji,jj) 
    536547               zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   & 
    537                   &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 
     548                  &           * ( gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk) ) * r1_e2v(ji,jj) 
    538549               ! 
    539550               IF( ln_wd_il ) THEN 
     
    545556               ! 
    546557               ! add to the general momentum trend 
    547                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
    548                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
     558               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) + zuap 
     559               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) + zvap 
    549560            END DO 
    550561         END DO 
     
    556567 
    557568 
    558    SUBROUTINE hpg_isf( kt ) 
     569   SUBROUTINE hpg_isf( kt, Kmm, puu, pvv, Krhs ) 
    559570      !!--------------------------------------------------------------------- 
    560571      !!                  ***  ROUTINE hpg_isf  *** 
     
    568579      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ] 
    569580      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ] 
    570       !!      add it to the general momentum trend (ua,va). 
    571       !!         ua = ua - 1/e1u * zhpi 
    572       !!         va = va - 1/e2v * zhpj 
     581      !!      add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 
     582      !!         puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 
     583      !!         pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 
    573584      !!      iceload is added and partial cell case are added to the top and bottom 
    574585      !!       
    575       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
    576       !!---------------------------------------------------------------------- 
    577       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     586      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 
     587      !!---------------------------------------------------------------------- 
     588      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index 
     589      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     590      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    578591      !! 
    579592      INTEGER  ::   ji, jj, jk, ikt, iktp1i, iktp1j   ! dummy loop indices 
     
    596609        DO jj = 1, jpj 
    597610          ikt = mikt(ji,jj) 
    598           zts_top(ji,jj,1) = tsn(ji,jj,ikt,1) 
    599           zts_top(ji,jj,2) = tsn(ji,jj,ikt,2) 
     611          zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm) 
     612          zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm) 
    600613        END DO 
    601614      END DO 
     
    612625            ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 
    613626            ! we assume ISF is in isostatic equilibrium 
    614             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj,iktp1i)                                    & 
     627            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w(ji+1,jj,iktp1i,Kmm)                                    & 
    615628               &                                    * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) )   & 
    616                &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         & 
     629               &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         & 
    617630               &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          & 
    618631               &                                  + ( riceload(ji+1,jj) - riceload(ji,jj))                            )  
    619             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji,jj+1,iktp1j)                                    & 
     632            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm)                                    & 
    620633               &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   & 
    621                &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         &  
     634               &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         &  
    622635               &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          & 
    623636               &                                  + ( riceload(ji,jj+1) - riceload(ji,jj))                            )  
    624637            ! s-coordinate pressure gradient correction (=0 if z coordinate) 
    625638            zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
    626                &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 
     639               &           * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 
    627640            zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
    628                &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 
     641               &           * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 
    629642            ! add to the general momentum trend 
    630             ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 
    631             va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 
     643            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 
     644            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 
    632645         END DO 
    633646      END DO 
     
    641654               ! hydrostatic pressure gradient along s-surfaces 
    642655               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
    643                   &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   & 
    644                   &              - e3w_n(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   ) 
     656                  &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   & 
     657                  &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   ) 
    645658               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
    646                   &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   & 
    647                   &              - e3w_n(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   ) 
     659                  &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   & 
     660                  &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   ) 
    648661               ! s-coordinate pressure gradient correction 
    649662               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    650                   &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj) 
     663                  &           * ( gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk) ) / e1u(ji,jj) 
    651664               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    652                   &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj) 
     665                  &           * ( gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk) ) / e2v(ji,jj) 
    653666               ! add to the general momentum trend 
    654                ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 
    655                va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 
     667               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 
     668               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 
    656669            END DO 
    657670         END DO 
     
    661674 
    662675 
    663    SUBROUTINE hpg_djc( kt ) 
     676   SUBROUTINE hpg_djc( kt, Kmm, puu, pvv, Krhs ) 
    664677      !!--------------------------------------------------------------------- 
    665678      !!                  ***  ROUTINE hpg_djc  *** 
     
    669682      !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003 
    670683      !!---------------------------------------------------------------------- 
    671       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     684      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index 
     685      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     686      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    672687      !! 
    673688      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     
    687702        DO jj = 2, jpjm1 
    688703           DO ji = 2, jpim1  
    689              ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
     704             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                & 
    690705                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            & 
    691                   &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     706                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) )  & 
    692707                  &                                                      > rn_wdmin1 + rn_wdmin2 
    693              ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (        & 
    694                   &    MAX(   sshn(ji,jj)             ,  sshn(ji+1,jj) ) >                & 
     708             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. (        & 
     709                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji+1,jj,Kmm) ) >                & 
    695710                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    696711             IF(ll_tmp1) THEN 
    697712               zcpx(ji,jj) = 1.0_wp 
    698713             ELSE IF(ll_tmp2) THEN 
    699                ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    700                zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
    701                            &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) ) 
     714               ! no worries about  ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm) = 0, it won't happen ! here 
     715               zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
     716                           &    / (ssh(ji+1,jj,Kmm) - ssh(ji  ,jj,Kmm)) ) 
    702717             ELSE 
    703718               zcpx(ji,jj) = 0._wp 
    704719             END IF 
    705720       
    706              ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
     721             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                & 
    707722                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            & 
    708                   &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
     723                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) )  & 
    709724                  &                                                      > rn_wdmin1 + rn_wdmin2 
    710              ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (        & 
    711                   &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                & 
     725             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. (        & 
     726                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji,jj+1,Kmm) ) >                & 
    712727                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    713728 
     
    715730               zcpy(ji,jj) = 1.0_wp 
    716731             ELSE IF(ll_tmp2) THEN 
    717                ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    718                zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
    719                            &    / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
     732               ! no worries about  ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm) = 0, it won't happen ! here 
     733               zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
     734                           &    / (ssh(ji,jj+1,Kmm) - ssh(ji,jj  ,Kmm)) ) 
    720735             ELSE 
    721736               zcpy(ji,jj) = 0._wp 
     
    747762            DO ji = fs_2, fs_jpim1   ! vector opt. 
    748763               drhoz(ji,jj,jk) = rhd    (ji  ,jj  ,jk) - rhd    (ji,jj,jk-1) 
    749                dzz  (ji,jj,jk) = gde3w_n(ji  ,jj  ,jk) - gde3w_n(ji,jj,jk-1) 
     764               dzz  (ji,jj,jk) = gde3w(ji  ,jj  ,jk) - gde3w(ji,jj,jk-1) 
    750765               drhox(ji,jj,jk) = rhd    (ji+1,jj  ,jk) - rhd    (ji,jj,jk  ) 
    751                dzx  (ji,jj,jk) = gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk  ) 
     766               dzx  (ji,jj,jk) = gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk  ) 
    752767               drhoy(ji,jj,jk) = rhd    (ji  ,jj+1,jk) - rhd    (ji,jj,jk  ) 
    753                dzy  (ji,jj,jk) = gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk  ) 
     768               dzy  (ji,jj,jk) = gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk  ) 
    754769            END DO 
    755770         END DO 
     
    838853      DO jj = 2, jpjm1 
    839854         DO ji = fs_2, fs_jpim1   ! vector opt. 
    840             rho_k(ji,jj,1) = -grav * ( e3w_n(ji,jj,1) - gde3w_n(ji,jj,1) )               & 
     855            rho_k(ji,jj,1) = -grav * ( e3w(ji,jj,1,Kmm) - gde3w(ji,jj,1) )               & 
    841856               &                   * (  rhd(ji,jj,1)                                     & 
    842857               &                     + 0.5_wp * ( rhd    (ji,jj,2) - rhd    (ji,jj,1) )  & 
    843                &                              * ( e3w_n  (ji,jj,1) - gde3w_n(ji,jj,1) )  & 
    844                &                              / ( gde3w_n(ji,jj,2) - gde3w_n(ji,jj,1) )  ) 
     858               &                              * ( e3w  (ji,jj,1,Kmm) - gde3w(ji,jj,1) )  & 
     859               &                              / ( gde3w(ji,jj,2) - gde3w(ji,jj,1) )  ) 
    845860         END DO 
    846861      END DO 
     
    854869 
    855870               rho_k(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj,jk) + rhd    (ji,jj,jk-1) )                                   & 
    856                   &                     * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) )                                   & 
     871                  &                     * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) )                                   & 
    857872                  &            - grav * z1_10 * (                                                                   & 
    858873                  &     ( drhow  (ji,jj,jk) - drhow  (ji,jj,jk-1) )                                                     & 
    859                   &   * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   & 
     874                  &   * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   & 
    860875                  &   - ( dzw    (ji,jj,jk) - dzw    (ji,jj,jk-1) )                                                     & 
    861876                  &   * ( rhd    (ji,jj,jk) - rhd    (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) )   & 
     
    863878 
    864879               rho_i(ji,jj,jk) = zcoef0 * ( rhd    (ji+1,jj,jk) + rhd    (ji,jj,jk) )                                   & 
    865                   &                     * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) )                                   & 
     880                  &                     * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )                                   & 
    866881                  &            - grav* z1_10 * (                                                                    & 
    867882                  &     ( drhou  (ji+1,jj,jk) - drhou  (ji,jj,jk) )                                                     & 
    868                   &   * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   & 
     883                  &   * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   & 
    869884                  &   - ( dzu    (ji+1,jj,jk) - dzu    (ji,jj,jk) )                                                     & 
    870885                  &   * ( rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) )   & 
     
    872887 
    873888               rho_j(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj+1,jk) + rhd    (ji,jj,jk) )                                 & 
    874                   &                     * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) )                                   & 
     889                  &                     * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) )                                   & 
    875890                  &            - grav* z1_10 * (                                                                    & 
    876891                  &     ( drhov  (ji,jj+1,jk) - drhov  (ji,jj,jk) )                                                     & 
    877                   &   * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   & 
     892                  &   * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   & 
    878893                  &   - ( dzv    (ji,jj+1,jk) - dzv    (ji,jj,jk) )                                                     & 
    879894                  &   * ( rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) )   & 
     
    897912            ENDIF 
    898913            ! add to the general momentum trend 
    899             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
    900             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 
     914            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 
     915            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 
    901916         END DO 
    902917      END DO 
     
    920935               ENDIF 
    921936               ! add to the general momentum trend 
    922                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
    923                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 
     937               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 
     938               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 
    924939            END DO 
    925940         END DO 
     
    931946 
    932947 
    933    SUBROUTINE hpg_prj( kt ) 
     948   SUBROUTINE hpg_prj( kt, Kmm, puu, pvv, Krhs ) 
    934949      !!--------------------------------------------------------------------- 
    935950      !!                  ***  ROUTINE hpg_prj  *** 
     
    940955      !!      all vertical coordinate systems 
    941956      !! 
    942       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
     957      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 
    943958      !!---------------------------------------------------------------------- 
    944959      INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear 
    945       INTEGER, INTENT(in) ::   kt                   ! ocean time-step index 
     960      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index 
     961      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     962      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    946963      !! 
    947964      INTEGER  ::   ji, jj, jk, jkk                 ! dummy loop indices 
     
    975992         DO jj = 2, jpjm1 
    976993           DO ji = 2, jpim1  
    977              ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
     994             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                & 
    978995                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            & 
    979                   &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     996                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) )  & 
    980997                  &                                                      > rn_wdmin1 + rn_wdmin2 
    981              ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (         & 
    982                   &    MAX(   sshn(ji,jj)             ,  sshn(ji+1,jj) ) >                & 
     998             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. (         & 
     999                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji+1,jj,Kmm) ) >                & 
    9831000                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    9841001 
     
    9861003               zcpx(ji,jj) = 1.0_wp 
    9871004             ELSE IF(ll_tmp2) THEN 
    988                ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    989                zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
    990                            &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
     1005               ! no worries about  ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm) = 0, it won't happen ! here 
     1006               zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
     1007                           &    / (ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm)) ) 
    9911008               
    9921009                zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 
     
    9951012             END IF 
    9961013       
    997              ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
     1014             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                & 
    9981015                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            & 
    999                   &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
     1016                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) )  & 
    10001017                  &                                                      > rn_wdmin1 + rn_wdmin2 
    1001              ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (      & 
    1002                   &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                & 
     1018             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. (      & 
     1019                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji,jj+1,Kmm) ) >                & 
    10031020                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    10041021 
     
    10061023               zcpy(ji,jj) = 1.0_wp 
    10071024             ELSE IF(ll_tmp2) THEN 
    1008                ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    1009                zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
    1010                            &    / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
     1025               ! no worries about  ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm) = 0, it won't happen ! here 
     1026               zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
     1027                           &    / (ssh(ji,jj+1,Kmm) - ssh(ji,jj  ,Kmm)) ) 
    10111028                zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 
    10121029 
     
    10311048          ELSEIF( jk < jpkm1 ) THEN 
    10321049             DO jkk = jk+1, jpk 
    1033                 zrhh(ji,jj,jkk) = interp1(gde3w_n(ji,jj,jkk  ), gde3w_n(ji,jj,jkk-1),   & 
    1034                    &                      gde3w_n(ji,jj,jkk-2), rhd    (ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 
     1050                zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk  ), gde3w(ji,jj,jkk-1),   & 
     1051                   &                      gde3w(ji,jj,jkk-2), rhd    (ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 
    10351052             END DO 
    10361053          ENDIF 
     
    10411058      DO jj = 1, jpj 
    10421059         DO ji = 1, jpi 
    1043             zdept(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) - sshn(ji,jj) * znad 
     1060            zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad 
    10441061         END DO 
    10451062      END DO 
     
    10481065         DO jj = 1, jpj 
    10491066            DO ji = 1, jpi 
    1050                zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w_n(ji,jj,jk) 
     1067               zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm) 
    10511068            END DO 
    10521069         END DO 
     
    10651082        DO ji = 2, jpi 
    10661083          zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1),  & 
    1067              &                                              csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1) 
     1084             &                                              csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm) 
    10681085 
    10691086          ! assuming linear profile across the top half surface layer 
    1070           zhpi(ji,jj,1) =  0.5_wp * e3w_n(ji,jj,1) * zrhdt1 
     1087          zhpi(ji,jj,1) =  0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt1 
    10711088        END DO 
    10721089      END DO 
     
    10901107        DO ji = 2, jpim1 
    10911108!!gm BUG ?    if it is ssh at u- & v-point then it should be: 
    1092 !          zsshu_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * & 
     1109!          zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * & 
    10931110!                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp  
    1094 !          zsshv_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * & 
     1111!          zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * & 
    10951112!                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp  
    10961113!!gm not this: 
    1097           zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * & 
     1114          zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 
    10981115                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp  
    1099           zsshv_n(ji,jj) = (e1e2v(ji,jj) * sshn(ji,jj) + e1e2v(ji+1, jj) * sshn(ji,jj+1)) * & 
     1116          zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 
    11001117                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp  
    11011118        END DO 
     
    11061123      DO jj = 2, jpjm1 
    11071124        DO ji = 2, jpim1 
    1108           zu(ji,jj,1) = - ( e3u_n(ji,jj,1) - zsshu_n(ji,jj) * znad)  
    1109           zv(ji,jj,1) = - ( e3v_n(ji,jj,1) - zsshv_n(ji,jj) * znad) 
     1125          zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad)  
     1126          zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad) 
    11101127        END DO 
    11111128      END DO 
     
    11141131        DO jj = 2, jpjm1 
    11151132          DO ji = 2, jpim1 
    1116             zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk) 
    1117             zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk) 
     1133            zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm) 
     1134            zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm) 
    11181135          END DO 
    11191136        END DO 
     
    11231140        DO jj = 2, jpjm1 
    11241141          DO ji = 2, jpim1 
    1125             zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u_n(ji,jj,jk) 
    1126             zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v_n(ji,jj,jk) 
     1142            zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm) 
     1143            zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm) 
    11271144          END DO 
    11281145        END DO 
     
    11761193               DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) 
    11771194                 IF( jk1 == 1 ) THEN 
    1178                    zdeps = zdept(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad) 
     1195                   zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad) 
    11791196                   zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), & 
    11801197                                                     bsp(jid,jj,1),   csp(jid,jj,1), & 
     
    11961213               IF( .NOT.ln_linssh ) THEN 
    11971214                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 
    1198                     &    ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 
     1215                    &    ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) ) 
    11991216                ELSE 
    12001217                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 
     
    12041221                  zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 
    12051222               ENDIF 
    1206                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
     1223               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
    12071224            ENDIF 
    12081225 
     
    12341251               DO WHILE ( -zdept(ji,jjd,jk1) < zvijk ) 
    12351252                 IF( jk1 == 1 ) THEN 
    1236                    zdeps = zdept(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad) 
     1253                   zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad) 
    12371254                   zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), & 
    12381255                                                     bsp(ji,jjd,1),   csp(ji,jjd,1), & 
     
    12551272               IF( .NOT.ln_linssh ) THEN 
    12561273                  zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 
    1257                           ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 
     1274                          ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) ) 
    12581275               ELSE 
    12591276                  zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 
     
    12641281               ENDIF 
    12651282 
    1266                va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 
     1283               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 
    12671284            ENDIF 
    12681285               ! 
Note: See TracChangeset for help on using the changeset viewer.