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 10877 – NEMO

Changeset 10877


Ignore:
Timestamp:
2019-04-16T11:50:59+02:00 (6 years ago)
Author:
davestorkey
Message:

branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : momentum advection and dynvor modules.

Location:
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv.F90

    r10874 r10877  
    5353CONTAINS 
    5454 
    55    SUBROUTINE dyn_adv( kt ) 
     55   SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs ) 
    5656      !!--------------------------------------------------------------------- 
    5757      !!                  ***  ROUTINE dyn_adv  *** 
    … …  
    5959      !! ** Purpose :   compute the ocean momentum advection trend. 
    6060      !! 
    61       !! ** Method  : - Update (ua,va) with the advection term following n_dynadv 
     61      !! ** Method  : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the advection term following n_dynadv 
    6262      !! 
    6363      !!      NB: in flux form advection (ln_dynadv_cen2 or ln_dynadv_ubs=T)  
    … …  
    6666      !!      (see dynvor module). 
    6767      !!---------------------------------------------------------------------- 
    68       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     68      INTEGER                             , INTENT( in )  ::  kt               ! ocean time-step index 
     69      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs   ! ocean time level indices 
     70      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation 
    6971      !!---------------------------------------------------------------------- 
    7072      ! 
    … …  
    7375      SELECT CASE( n_dynadv )    !==  compute advection trend and add it to general trend  ==! 
    7476      CASE( np_VEC_c2  )      
    75          CALL dyn_keg     ( kt, nn_dynkeg )    ! vector form : horizontal gradient of kinetic energy 
    76          CALL dyn_zad     ( kt )               ! vector form : vertical advection 
     77         CALL dyn_keg     ( kt, nn_dynkeg, Kmm, puu, pvv, Krhs )    ! vector form : horizontal gradient of kinetic energy 
     78         CALL dyn_zad     ( kt,            Kmm, puu, pvv, Krhs )    ! vector form : vertical advection 
    7779      CASE( np_FLX_c2  )  
    78          CALL dyn_adv_cen2( kt )               ! 2nd order centered scheme 
     80         CALL dyn_adv_cen2( kt,      Kmm, puu, pvv, Krhs )          ! 2nd order centered scheme 
    7981      CASE( np_FLX_ubs )    
    80          CALL dyn_adv_ubs ( kt )               ! 3rd order UBS      scheme (UP3) 
     82         CALL dyn_adv_ubs ( kt, Kbb, Kmm, puu, pvv, Krhs )          ! 3rd order UBS      scheme (UP3) 
    8183      END SELECT 
    8284      ! 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_cen2.F90

    r10874 r10877  
    3535CONTAINS 
    3636 
    37    SUBROUTINE dyn_adv_cen2( kt ) 
     37   SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs ) 
    3838      !!---------------------------------------------------------------------- 
    3939      !!                  ***  ROUTINE dyn_adv_cen2  *** 
    … …  
    4444      !! ** Method  :   Trend evaluated using now fields (centered in time)  
    4545      !! 
    46       !! ** Action  :   (ua,va) updated with the now vorticity term trend 
     46      !! ** Action  :   (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the now vorticity term trend 
    4747      !!---------------------------------------------------------------------- 
    48       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     48      INTEGER                             , INTENT( in )  ::  kt           ! ocean time-step index 
     49      INTEGER                             , INTENT( in )  ::  Kmm, Krhs    ! ocean time level indices 
     50      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv     ! ocean velocities and RHS of momentum equation 
    4951      ! 
    5052      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    … …  
    6062      ! 
    6163      IF( l_trddyn ) THEN           ! trends: store the input trends 
    62          zfu_uw(:,:,:) = ua(:,:,:) 
    63          zfv_vw(:,:,:) = va(:,:,:) 
     64         zfu_uw(:,:,:) = puu(:,:,:,Krhs) 
     65         zfv_vw(:,:,:) = pvv(:,:,:,Krhs) 
    6466      ENDIF 
    6567      ! 
    … …  
    6769      ! 
    6870      DO jk = 1, jpkm1                    ! horizontal transport 
    69          zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
    70          zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
     71         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) 
     72         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) 
    7173         DO jj = 1, jpjm1                 ! horizontal momentum fluxes (at T- and F-point) 
    7274            DO ji = 1, fs_jpim1   ! vector opt. 
    73                zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji+1,jj  ,jk) ) 
    74                zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) ) 
    75                zfu_f(ji  ,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) ) 
    76                zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) ) 
     75               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) ) 
     76               zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) ) 
     77               zfu_f(ji  ,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) ) 
     78               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) ) 
    7779            END DO 
    7880         END DO 
    7981         DO jj = 2, jpjm1                 ! divergence of horizontal momentum fluxes 
    8082            DO ji = fs_2, fs_jpim1   ! vector opt. 
    81                ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    & 
    82                   &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
    83                va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    & 
    84                   &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
     83               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    & 
     84                  &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
     85               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    & 
     86                  &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
    8587            END DO 
    8688         END DO 
    … …  
    8890      ! 
    8991      IF( l_trddyn ) THEN           ! trends: send trend to trddyn for diagnostic 
    90          zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 
    91          zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 
     92         zfu_uw(:,:,:) = puu(:,:,:,Krhs) - zfu_uw(:,:,:) 
     93         zfv_vw(:,:,:) = pvv(:,:,:,Krhs) - zfv_vw(:,:,:) 
    9294         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt ) 
    93          zfu_t(:,:,:) = ua(:,:,:) 
    94          zfv_t(:,:,:) = va(:,:,:) 
     95         zfu_t(:,:,:) = puu(:,:,:,Krhs) 
     96         zfv_t(:,:,:) = pvv(:,:,:,Krhs) 
    9597      ENDIF 
    9698      ! 
    … …  
    106108         DO jj = 2, jpjm1 
    107109            DO ji = fs_2, fs_jpim1 
    108                zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1) 
    109                zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1) 
     110               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm) 
     111               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm) 
    110112            END DO 
    111113         END DO 
    … …  
    114116         DO jj = 2, jpj                       ! 1/4 * Vertical transport 
    115117            DO ji = 2, jpi 
    116                zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
     118               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk) 
    117119            END DO 
    118120         END DO 
    119121         DO jj = 2, jpjm1 
    120122            DO ji = fs_2, fs_jpim1   ! vector opt. 
    121                zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 
    122                zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 
     123               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj  ,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) ) 
     124               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji  ,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) ) 
    123125            END DO 
    124126         END DO 
    … …  
    127129         DO jj = 2, jpjm1  
    128130            DO ji = fs_2, fs_jpim1   ! vector opt. 
    129                ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
    130                va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
     131               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
     132               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
    131133            END DO 
    132134         END DO 
    … …  
    134136      ! 
    135137      IF( l_trddyn ) THEN                 ! trends: send trend to trddyn for diagnostic 
    136          zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 
    137          zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 
     138         zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:) 
     139         zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:) 
    138140         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 
    139141      ENDIF 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90

    r10874 r10877  
    4141CONTAINS 
    4242 
    43    SUBROUTINE dyn_adv_ubs( kt ) 
     43   SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs ) 
    4444      !!---------------------------------------------------------------------- 
    4545      !!                  ***  ROUTINE dyn_adv_ubs  *** 
    … …  
    6464      !!      gamma1=1/3 and gamma2=1/32. 
    6565      !! 
    66       !! ** Action : - (ua,va) updated with the 3D advective momentum trends 
     66      !! ** Action : - (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the 3D advective momentum trends 
    6767      !! 
    6868      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.  
    6969      !!---------------------------------------------------------------------- 
    70       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     70      INTEGER                             , INTENT( in )  ::  kt              ! ocean time-step index 
     71      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs  ! ocean time level indices 
     72      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv        ! ocean velocities and RHS of momentum equation 
    7173      ! 
    7274      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    … …  
    9597      ! 
    9698      IF( l_trddyn ) THEN           ! trends: store the input trends 
    97          zfu_uw(:,:,:) = ua(:,:,:) 
    98          zfv_vw(:,:,:) = va(:,:,:) 
     99         zfu_uw(:,:,:) = puu(:,:,:,Krhs) 
     100         zfv_vw(:,:,:) = pvv(:,:,:,Krhs) 
    99101      ENDIF 
    100102      !                                      ! =========================== ! 
    … …  
    102104         !                                   ! =========================== ! 
    103105         !                                         ! horizontal volume fluxes 
    104          zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
    105          zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
     106         zfu(:,:,jk) = e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) 
     107         zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) 
    106108         !             
    107109         DO jj = 2, jpjm1                          ! laplacian 
    108110            DO ji = fs_2, fs_jpim1   ! vector opt. 
    109                zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj  ,jk) - 2.*ub (ji,jj,jk) + ub (ji-1,jj  ,jk) ) * umask(ji,jj,jk) 
    110                zlv_vv(ji,jj,jk,1) = ( vb (ji  ,jj+1,jk) - 2.*vb (ji,jj,jk) + vb (ji  ,jj-1,jk) ) * vmask(ji,jj,jk) 
    111                zlu_uv(ji,jj,jk,1) = ( ub (ji  ,jj+1,jk) - ub (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   & 
    112                   &               - ( ub (ji  ,jj  ,jk) - ub (ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk) 
    113                zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj  ,jk) - vb (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   & 
    114                   &               - ( vb (ji  ,jj  ,jk) - vb (ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk) 
     111               zlu_uu(ji,jj,jk,1) = ( puu (ji+1,jj  ,jk,Kbb) - 2.*puu (ji,jj,jk,Kbb) + puu (ji-1,jj  ,jk,Kbb) ) * umask(ji,jj,jk) 
     112               zlv_vv(ji,jj,jk,1) = ( pvv (ji  ,jj+1,jk,Kbb) - 2.*pvv (ji,jj,jk,Kbb) + pvv (ji  ,jj-1,jk,Kbb) ) * vmask(ji,jj,jk) 
     113               zlu_uv(ji,jj,jk,1) = ( puu (ji  ,jj+1,jk,Kbb) - puu (ji  ,jj  ,jk,Kbb) ) * fmask(ji  ,jj  ,jk)   & 
     114                  &               - ( puu (ji  ,jj  ,jk,Kbb) - puu (ji  ,jj-1,jk,Kbb) ) * fmask(ji  ,jj-1,jk) 
     115               zlv_vu(ji,jj,jk,1) = ( pvv (ji+1,jj  ,jk,Kbb) - pvv (ji  ,jj  ,jk,Kbb) ) * fmask(ji  ,jj  ,jk)   & 
     116                  &               - ( pvv (ji  ,jj  ,jk,Kbb) - pvv (ji-1,jj  ,jk,Kbb) ) * fmask(ji-1,jj  ,jk) 
    115117               ! 
    116118               zlu_uu(ji,jj,jk,2) = ( zfu(ji+1,jj  ,jk) - 2.*zfu(ji,jj,jk) + zfu(ji-1,jj  ,jk) ) * umask(ji,jj,jk) 
    … …  
    132134      DO jk = 1, jpkm1                       ! ====================== ! 
    133135         !                                         ! horizontal volume fluxes 
    134          zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
    135          zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
     136         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) 
     137         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) 
    136138         ! 
    137139         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point 
    138140            DO ji = 1, fs_jpim1   ! vector opt. 
    139                zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) ) 
    140                zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) ) 
     141               zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) ) 
     142               zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) ) 
    141143               ! 
    142144               IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1) 
    … …  
    164166               ! 
    165167               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   & 
    166                   &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u ) 
     168                  &                * ( puu(ji,jj,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) - gamma1 * zl_u ) 
    167169               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   & 
    168                   &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v ) 
     170                  &                * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) - gamma1 * zl_v ) 
    169171            END DO 
    170172         END DO 
    171173         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes 
    172174            DO ji = fs_2, fs_jpim1   ! vector opt. 
    173                ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    & 
    174                   &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
    175                va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    & 
    176                   &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
     175               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    & 
     176                  &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
     177               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    & 
     178                  &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
    177179            END DO 
    178180         END DO 
    179181      END DO 
    180182      IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic 
    181          zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 
    182          zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 
     183         zfu_uw(:,:,:) = puu(:,:,:,Krhs) - zfu_uw(:,:,:) 
     184         zfv_vw(:,:,:) = pvv(:,:,:,Krhs) - zfv_vw(:,:,:) 
    183185         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt ) 
    184          zfu_t(:,:,:) = ua(:,:,:) 
    185          zfv_t(:,:,:) = va(:,:,:) 
     186         zfu_t(:,:,:) = puu(:,:,:,Krhs) 
     187         zfv_t(:,:,:) = pvv(:,:,:,Krhs) 
    186188      ENDIF 
    187189      !                                      ! ==================== ! 
    … …  
    199201         DO jj = 2, jpjm1 
    200202            DO ji = fs_2, fs_jpim1 
    201                zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1) 
    202                zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1) 
     203               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm) 
     204               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm) 
    203205            END DO 
    204206         END DO 
    … …  
    207209         DO jj = 2, jpj 
    208210            DO ji = 2, jpi 
    209                zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
     211               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk) 
    210212            END DO 
    211213         END DO 
    212214         DO jj = 2, jpjm1 
    213215            DO ji = fs_2, fs_jpim1   ! vector opt. 
    214                zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 
    215                zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 
     216               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) ) 
     217               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) ) 
    216218            END DO 
    217219         END DO 
    … …  
    220222         DO jj = 2, jpjm1 
    221223            DO ji = fs_2, fs_jpim1   ! vector opt. 
    222                ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
    223                va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
     224               puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
     225               pvv(ji,jj,jk,Krhs) =  pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
    224226            END DO 
    225227         END DO 
    … …  
    227229      ! 
    228230      IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic 
    229          zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 
    230          zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 
     231         zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:) 
     232         zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:) 
    231233         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 
    232234      ENDIF 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynkeg.F90

    r10874 r10877  
    4444CONTAINS 
    4545 
    46    SUBROUTINE dyn_keg( kt, kscheme ) 
     46   SUBROUTINE dyn_keg( kt, kscheme, Kmm, puu, pvv, Krhs ) 
    4747      !!---------------------------------------------------------------------- 
    4848      !!                  ***  ROUTINE dyn_keg  *** 
    … …  
    5757      !!              * kscheme = nkeg_HW : Hollingsworth correction following 
    5858      !!      Arakawa (2001). The now horizontal kinetic energy is given by: 
    59       !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((un(j+1)+un(j-1))/2)^2  ) 
    60       !!                    + mj-1(  2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2  ) ] 
     59      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((puu(j+1,Kmm)+puu(j-1,Kmm))/2)^2  ) 
     60      !!                    + mj-1(  2 * vn^2 + ((pvv(i+1,Kmm)+pvv(i-1,Kmm))/2)^2  ) ] 
    6161      !!       
    6262      !!      Take its horizontal gradient and add it to the general momentum 
    63       !!      trend (ua,va). 
    64       !!         ua = ua - 1/e1u di[ zhke ] 
    65       !!         va = va - 1/e2v dj[ zhke ] 
     63      !!      trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 
     64      !!         puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u di[ zhke ] 
     65      !!         pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v dj[ zhke ] 
    6666      !! 
    67       !! ** Action : - Update the (ua, va) with the hor. ke gradient trend 
     67      !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend 
    6868      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing 
    6969      !! 
    … …  
    7171      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983. 
    7272      !!---------------------------------------------------------------------- 
    73       INTEGER, INTENT( in ) ::   kt        ! ocean time-step index 
    74       INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme  
     73      INTEGER                             , INTENT( in )  ::  kt               ! ocean time-step index 
     74      INTEGER                             , INTENT( in )  ::  kscheme          ! =0/1   type of KEG scheme  
     75      INTEGER                             , INTENT( in )  ::  Kmm, Krhs        ! ocean time level indices 
     76      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation 
    7577      ! 
    7678      INTEGER  ::   ji, jj, jk, jb    ! dummy loop indices 
    … …  
    9294      IF( l_trddyn ) THEN           ! Save the input trends 
    9395         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
    94          ztrdu(:,:,:) = ua(:,:,:)  
    95          ztrdv(:,:,:) = va(:,:,:)  
     96         ztrdu(:,:,:) = puu(:,:,:,Krhs)  
     97         ztrdv(:,:,:) = pvv(:,:,:,Krhs)  
    9698      ENDIF 
    9799       
    … …  
    109111                     ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
    110112                     ifu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) ) 
    111                      un(ii-ifu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk) 
     113                     puu(ii-ifu,ij,jk,Kmm) = puu(ii,ij,jk,Kmm) * umask(ii,ij,jk) 
    112114                  END DO 
    113115               END DO 
    … …  
    119121                     ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
    120122                     ifv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) ) 
    121                      vn(ii,ij-ifv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk) 
     123                     pvv(ii,ij-ifv,jk,Kmm) = pvv(ii,ij,jk,Kmm) * vmask(ii,ij,jk) 
    122124                  END DO 
    123125               END DO 
    … …  
    132134            DO jj = 2, jpj 
    133135               DO ji = fs_2, jpi   ! vector opt. 
    134                   zu =    un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   & 
    135                      &  + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) 
    136                   zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   & 
    137                      &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) 
     136                  zu =    puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)   & 
     137                     &  + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm) 
     138                  zv =    pvv(ji  ,jj-1,jk,Kmm) * pvv(ji  ,jj-1,jk,Kmm)   & 
     139                     &  + pvv(ji  ,jj  ,jk,Kmm) * pvv(ji  ,jj  ,jk,Kmm) 
    138140                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 
    139141               END DO   
    … …  
    145147            DO jj = 2, jpjm1        
    146148               DO ji = fs_2, jpim1   ! vector opt. 
    147                   zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    & 
    148                      &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  & 
    149                      &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   & 
    150                      &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) 
     149                  zu = 8._wp * ( puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)    & 
     150                     &         + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm) )  & 
     151                     &   +     ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) * ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) )   & 
     152                     &   +     ( puu(ji  ,jj-1,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) ) * ( puu(ji  ,jj-1,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) ) 
    151153                     ! 
    152                   zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    & 
    153                      &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  & 
    154                      &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   & 
    155                      &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) 
     154                  zv = 8._wp * ( pvv(ji  ,jj-1,jk,Kmm) * pvv(ji  ,jj-1,jk,Kmm)    & 
     155                     &         + pvv(ji  ,jj  ,jk,Kmm) * pvv(ji  ,jj  ,jk,Kmm) )  & 
     156                     &  +      ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) * ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) )   & 
     157                     &  +      ( pvv(ji-1,jj  ,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) ) * ( pvv(ji-1,jj  ,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) ) 
    156158                  zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 
    157159               END DO   
    … …  
    164166      IF (ln_bdy) THEN 
    165167         ! restore velocity masks at points outside boundary 
    166          un(:,:,:) = un(:,:,:) * umask(:,:,:) 
    167          vn(:,:,:) = vn(:,:,:) * vmask(:,:,:) 
     168         puu(:,:,:,Kmm) = puu(:,:,:,Kmm) * umask(:,:,:) 
     169         pvv(:,:,:,Kmm) = pvv(:,:,:,Kmm) * vmask(:,:,:) 
    168170      ENDIF       
    169171 
    … …  
    172174         DO jj = 2, jpjm1 
    173175            DO ji = fs_2, fs_jpim1   ! vector opt. 
    174                ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 
    175                va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 
     176               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 
     177               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 
    176178            END DO  
    177179         END DO 
    … …  
    179181      ! 
    180182      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic 
    181          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    182          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     183         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 
     184         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 
    183185         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt ) 
    184186         DEALLOCATE( ztrdu , ztrdv ) 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90

    r10874 r10877  
    9696CONTAINS 
    9797 
    98    SUBROUTINE dyn_vor( kt ) 
     98   SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs ) 
    9999      !!---------------------------------------------------------------------- 
    100100      !! 
    101101      !! ** Purpose :   compute the lateral ocean tracer physics. 
    102102      !! 
    103       !! ** Action : - Update (ua,va) with the now vorticity term trend 
     103      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend 
    104104      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative 
    105105      !!               and planetary vorticity trends) and send them to trd_dyn  
    106106      !!               for futher diagnostics (l_trddyn=T) 
    107107      !!---------------------------------------------------------------------- 
    108       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     108      INTEGER                             , INTENT( in  ) ::   kt          ! ocean time-step index 
     109      INTEGER                             , INTENT( in  ) ::   Kmm, Krhs   ! ocean time level indices 
     110      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv    ! ocean velocity field and RHS of momentum equation 
    109111      ! 
    110112      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
    … …  
    117119         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
    118120         ! 
    119          ztrdu(:,:,:) = ua(:,:,:)            !* planetary vorticity trend (including Stokes-Coriolis force) 
    120          ztrdv(:,:,:) = va(:,:,:) 
     121         ztrdu(:,:,:) = puu(:,:,:,Krhs)            !* planetary vorticity trend (including Stokes-Coriolis force) 
     122         ztrdv(:,:,:) = pvv(:,:,:,Krhs) 
    121123         SELECT CASE( nvor_scheme ) 
    122          CASE( np_ENS )           ;   CALL vor_ens( kt, ncor, un , vn , ua, va )   ! enstrophy conserving scheme 
    123             IF( ln_stcor )            CALL vor_ens( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    124          CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, ncor, un , vn , ua, va )   ! energy conserving scheme 
    125             IF( ln_stcor )            CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    126          CASE( np_ENT )           ;   CALL vor_enT( kt, ncor, un , vn , ua, va )   ! energy conserving scheme (T-pts) 
    127             IF( ln_stcor )            CALL vor_enT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    128          CASE( np_EET )           ;   CALL vor_eeT( kt, ncor, un , vn , ua, va )   ! energy conserving scheme (een with e3t) 
    129             IF( ln_stcor )            CALL vor_eeT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    130          CASE( np_EEN )           ;   CALL vor_een( kt, ncor, un , vn , ua, va )   ! energy & enstrophy scheme 
    131             IF( ln_stcor )            CALL vor_een( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
     124         CASE( np_ENS )           ;   CALL vor_ens( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! enstrophy conserving scheme 
     125            IF( ln_stcor )            CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
     126         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme 
     127            IF( ln_stcor )            CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
     128         CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (T-pts) 
     129            IF( ln_stcor )            CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
     130         CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (een with e3t) 
     131            IF( ln_stcor )            CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
     132         CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy & enstrophy scheme 
     133            IF( ln_stcor )            CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
    132134         END SELECT 
    133          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    134          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     135         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 
     136         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 
    135137         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    136138         ! 
    137139         IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case) 
    138             ztrdu(:,:,:) = ua(:,:,:) 
    139             ztrdv(:,:,:) = va(:,:,:) 
     140            ztrdu(:,:,:) = puu(:,:,:,Krhs) 
     141            ztrdv(:,:,:) = pvv(:,:,:,Krhs) 
    140142            SELECT CASE( nvor_scheme ) 
    141             CASE( np_ENT )           ;   CALL vor_enT( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (T-pts) 
    142             CASE( np_EET )           ;   CALL vor_eeT( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (een with e3t) 
    143             CASE( np_ENE )           ;   CALL vor_ene( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme 
    144             CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, nrvm, un , vn , ua, va )  ! enstrophy conserving scheme 
    145             CASE( np_EEN )           ;   CALL vor_een( kt, nrvm, un , vn , ua, va )  ! energy & enstrophy scheme 
     143            CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (T-pts) 
     144            CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (een with e3t) 
     145            CASE( np_ENE )           ;   CALL vor_ene( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme 
     146            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! enstrophy conserving scheme 
     147            CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy & enstrophy scheme 
    146148            END SELECT 
    147             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    148             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     149            ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 
     150            ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 
    149151            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 
    150152         ENDIF 
    … …  
    156158         SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==! 
    157159         CASE( np_ENT )                        !* energy conserving scheme  (T-pts) 
    158                              CALL vor_enT( kt, ntot, un , vn , ua, va )   ! total vorticity trend 
    159             IF( ln_stcor )   CALL vor_enT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
     160                             CALL vor_enT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend 
     161            IF( ln_stcor )   CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
    160162         CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t) 
    161                              CALL vor_eeT( kt, ntot, un , vn , ua, va )   ! total vorticity trend 
    162             IF( ln_stcor )   CALL vor_eeT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
     163                             CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend 
     164            IF( ln_stcor )   CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
    163165         CASE( np_ENE )                        !* energy conserving scheme 
    164                              CALL vor_ene( kt, ntot, un , vn , ua, va )   ! total vorticity trend 
    165             IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
     166                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend 
     167            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
    166168         CASE( np_ENS )                        !* enstrophy conserving scheme 
    167                              CALL vor_ens( kt, ntot, un , vn , ua, va )  ! total vorticity trend 
    168             IF( ln_stcor )   CALL vor_ens( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend 
     169                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend 
     170            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend 
    169171         CASE( np_MIX )                        !* mixed ene-ens scheme 
    170                              CALL vor_ens( kt, nrvm, un , vn , ua, va )   ! relative vorticity or metric trend (ens) 
    171                              CALL vor_ene( kt, ncor, un , vn , ua, va )   ! planetary vorticity trend (ene) 
    172             IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
     172                             CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! relative vorticity or metric trend (ens) 
     173                             CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! planetary vorticity trend (ene) 
     174            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
    173175         CASE( np_EEN )                        !* energy and enstrophy conserving scheme 
    174                              CALL vor_een( kt, ntot, un , vn , ua, va )   ! total vorticity trend 
    175             IF( ln_stcor )   CALL vor_een( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
     176                             CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend 
     177            IF( ln_stcor )   CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
    176178         END SELECT 
    177179         ! 
    … …  
    187189 
    188190 
    189    SUBROUTINE vor_enT( kt, kvor, pu, pv, pu_rhs, pv_rhs ) 
     191   SUBROUTINE vor_enT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs ) 
    190192      !!---------------------------------------------------------------------- 
    191193      !!                  ***  ROUTINE vor_enT  *** 
    … …  
    203205      !!       where rvor is the relative vorticity at f-point 
    204206      !! 
    205       !! ** Action : - Update (ua,va) with the now vorticity term trend 
     207      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend 
    206208      !!---------------------------------------------------------------------- 
    207209      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index 
     210      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index 
    208211      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric 
    209212      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities 
    … …  
    270273         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==! 
    271274         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
    272             zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t_n(:,:,jk) 
     275            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm) 
    273276         CASE ( np_RVO )                           !* relative vorticity 
    274277            DO jj = 2, jpj 
    275278               DO ji = 2, jpi   ! vector opt. 
    276279                  zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   & 
    277                      &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
     280                     &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
    278281               END DO 
    279282            END DO 
    … …  
    282285               DO ji = 2, jpi 
    283286                  zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   & 
    284                      &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t_n(ji,jj,jk) 
     287                     &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t(ji,jj,jk,Kmm) 
    285288               END DO 
    286289            END DO 
    … …  
    289292               DO ji = 2, jpi   ! vector opt. 
    290293                  zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    & 
    291                      &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
     294                     &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
    292295               END DO 
    293296            END DO 
    … …  
    297300                  zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           & 
    298301                       &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  & 
    299                        &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t_n(ji,jj,jk) 
     302                       &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t(ji,jj,jk,Kmm) 
    300303               END DO 
    301304            END DO 
    … …  
    307310         DO jj = 2, jpjm1 
    308311            DO ji = 2, jpim1   ! vector opt. 
    309                pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)                    & 
     312               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                    & 
    310313                  &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   & 
    311314                  &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   ) 
    312315                  ! 
    313                pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)                    & 
     316               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    & 
    314317                  &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   &  
    315318                  &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   )  
    … …  
    322325 
    323326 
    324    SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva ) 
     327   SUBROUTINE vor_ene( kt, Kmm, kvor, pu_mm, pv_mm, pu_rhs, pv_rhs ) 
    325328      !!---------------------------------------------------------------------- 
    326329      !!                  ***  ROUTINE vor_ene  *** 
    … …  
    334337      !!         The general trend of momentum is increased due to the vorticity  
    335338      !!       term which is given by: 
    336       !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v vn) ] 
    337       !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u un) ] 
     339      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v pvv(:,:,:,Kmm)) ] 
     340      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u puu(:,:,:,Kmm)) ] 
    338341      !!       where rvor is the relative vorticity 
    339342      !! 
    340       !! ** Action : - Update (ua,va) with the now vorticity term trend 
     343      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend 
    341344      !! 
    342345      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    343346      !!---------------------------------------------------------------------- 
    344347      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index 
     348      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index 
    345349      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
    346       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities 
    347       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend 
     350      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_mm, pv_mm    ! now velocities 
     351      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend 
    348352      ! 
    349353      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    … …  
    368372            DO jj = 1, jpjm1 
    369373               DO ji = 1, fs_jpim1   ! vector opt. 
    370                   zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
    371                      &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     374                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv_mm(ji+1,jj  ,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)    & 
     375                     &          - e1u(ji  ,jj+1) * pu_mm(ji  ,jj+1,jk) + e1u(ji,jj) * pu_mm(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    372376               END DO 
    373377            END DO 
    … …  
    375379            DO jj = 1, jpjm1 
    376380               DO ji = 1, fs_jpim1   ! vector opt. 
    377                   zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
    378                      &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
     381                  zwz(ji,jj) = ( pv_mm(ji+1,jj  ,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     382                     &       - ( pu_mm(ji  ,jj+1,jk) + pu_mm(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    379383               END DO 
    380384            END DO 
    … …  
    382386            DO jj = 1, jpjm1 
    383387               DO ji = 1, fs_jpim1   ! vector opt. 
    384                   zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      & 
    385                      &                        - e1u(ji,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     388                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pv_mm(ji+1,jj,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)      & 
     389                     &                        - e1u(ji,jj+1) * pu_mm(ji,jj+1,jk) + e1u(ji,jj) * pu_mm(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    386390               END DO 
    387391            END DO 
    … …  
    389393            DO jj = 1, jpjm1 
    390394               DO ji = 1, fs_jpim1   ! vector opt. 
    391                   zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
    392                      &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
     395                  zwz(ji,jj) = ff_f(ji,jj) + ( pv_mm(ji+1,jj  ,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     396                     &                     - ( pu_mm(ji  ,jj+1,jk) + pu_mm(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    393397               END DO 
    394398            END DO 
    … …  
    406410 
    407411         IF( ln_sco ) THEN 
    408             zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
    409             zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
    410             zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
     412            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     413            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu_mm(:,:,jk) 
     414            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv_mm(:,:,jk) 
    411415         ELSE 
    412             zwx(:,:) = e2u(:,:) * pun(:,:,jk) 
    413             zwy(:,:) = e1v(:,:) * pvn(:,:,jk) 
     416            zwx(:,:) = e2u(:,:) * pu_mm(:,:,jk) 
     417            zwy(:,:) = e1v(:,:) * pv_mm(:,:,jk) 
    414418         ENDIF 
    415419         !                                   !==  compute and add the vorticity term trend  =! 
    … …  
    420424               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 
    421425               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
    422                pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    423                pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
     426               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     427               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
    424428            END DO   
    425429         END DO   
    … …  
    430434 
    431435 
    432    SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva ) 
     436   SUBROUTINE vor_ens( kt, Kmm, kvor, pu_mm, pv_mm, pu_rhs, pv_rhs ) 
    433437      !!---------------------------------------------------------------------- 
    434438      !!                ***  ROUTINE vor_ens  *** 
    … …  
    441445      !!      potential enstrophy of a horizontally non-divergent flow. the 
    442446      !!      trend of the vorticity term is given by: 
    443       !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ] 
    444       !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u un) ] 
    445       !!      Add this trend to the general momentum trend (ua,va): 
    446       !!          (ua,va) = (ua,va) + ( voru , vorv ) 
    447       !! 
    448       !! ** Action : - Update (ua,va) arrays with the now vorticity term trend 
     447      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v pvv(:,:,:,Kmm)) ] 
     448      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u puu(:,:,:,Kmm)) ] 
     449      !!      Add this trend to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)): 
     450      !!          (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) = (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) + ( voru , vorv ) 
     451      !! 
     452      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) arrays with the now vorticity term trend 
    449453      !! 
    450454      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    451455      !!---------------------------------------------------------------------- 
    452456      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index 
     457      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index 
    453458      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
    454       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities 
    455       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend 
     459      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_mm, pv_mm    ! now velocities 
     460      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend 
    456461      ! 
    457462      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    … …  
    475480            DO jj = 1, jpjm1 
    476481               DO ji = 1, fs_jpim1   ! vector opt. 
    477                   zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
    478                      &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     482                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv_mm(ji+1,jj  ,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)    & 
     483                     &          - e1u(ji  ,jj+1) * pu_mm(ji  ,jj+1,jk) + e1u(ji,jj) * pu_mm(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    479484               END DO 
    480485            END DO 
    … …  
    482487            DO jj = 1, jpjm1 
    483488               DO ji = 1, fs_jpim1   ! vector opt. 
    484                   zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
    485                      &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
     489                  zwz(ji,jj) = ( pv_mm(ji+1,jj  ,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     490                     &       - ( pu_mm(ji  ,jj+1,jk) + pu_mm(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    486491               END DO 
    487492            END DO 
    … …  
    489494            DO jj = 1, jpjm1 
    490495               DO ji = 1, fs_jpim1   ! vector opt. 
    491                   zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)  & 
    492                      &                        - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     496                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv_mm(ji+1,jj  ,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)  & 
     497                     &                        - e1u(ji  ,jj+1) * pu_mm(ji  ,jj+1,jk) + e1u(ji,jj) * pu_mm(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    493498               END DO 
    494499            END DO 
    … …  
    496501            DO jj = 1, jpjm1 
    497502               DO ji = 1, fs_jpim1   ! vector opt. 
    498                   zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
    499                      &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
     503                  zwz(ji,jj) = ff_f(ji,jj) + ( pv_mm(ji+1,jj  ,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     504                     &                     - ( pu_mm(ji  ,jj+1,jk) + pu_mm(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    500505               END DO 
    501506            END DO 
    … …  
    513518         ! 
    514519         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==! 
    515             zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
    516             zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
    517             zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
     520            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     521            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu_mm(:,:,jk) 
     522            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv_mm(:,:,jk) 
    518523         ELSE 
    519             zwx(:,:) = e2u(:,:) * pun(:,:,jk) 
    520             zwy(:,:) = e1v(:,:) * pvn(:,:,jk) 
     524            zwx(:,:) = e2u(:,:) * pu_mm(:,:,jk) 
     525            zwy(:,:) = e1v(:,:) * pv_mm(:,:,jk) 
    521526         ENDIF 
    522527         !                                   !==  compute and add the vorticity term trend  =! 
    … …  
    527532               zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  & 
    528533                  &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  ) 
    529                pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
    530                pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
     534               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
     535               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
    531536            END DO   
    532537         END DO   
    … …  
    537542 
    538543 
    539    SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva ) 
     544   SUBROUTINE vor_een( kt, Kmm, kvor, pu_mm, pv_mm, pu_rhs, pv_rhs ) 
    540545      !!---------------------------------------------------------------------- 
    541546      !!                ***  ROUTINE vor_een  *** 
    … …  
    548553      !!      both the horizontal kinetic energy and the potential enstrophy 
    549554      !!      when horizontal divergence is zero (see the NEMO documentation) 
    550       !!      Add this trend to the general momentum trend (ua,va). 
    551       !! 
    552       !! ** Action : - Update (ua,va) with the now vorticity term trend 
     555      !!      Add this trend to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 
     556      !! 
     557      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend 
    553558      !! 
    554559      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 
    555560      !!---------------------------------------------------------------------- 
    556561      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index 
     562      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index 
    557563      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
    558       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities 
    559       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend 
     564      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_mm, pv_mm    ! now velocities 
     565      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend 
    560566      ! 
    561567      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    … …  
    582588            DO jj = 1, jpjm1 
    583589               DO ji = 1, fs_jpim1   ! vector opt. 
    584                   ze3f = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
    585                      &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  ) 
     590                  ze3f = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     591                     &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
    586592                  IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f 
    587593                  ELSE                       ;   z1_e3f(ji,jj) = 0._wp 
    … …  
    592598            DO jj = 1, jpjm1 
    593599               DO ji = 1, fs_jpim1   ! vector opt. 
    594                   ze3f = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
    595                      &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  ) 
     600                  ze3f = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     601                     &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
    596602                  zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   & 
    597603                     &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  ) 
    … …  
    613619            DO jj = 1, jpjm1 
    614620               DO ji = 1, fs_jpim1   ! vector opt. 
    615                   zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)  & 
    616                      &            - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj) 
     621                  zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pv_mm(ji+1,jj,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)  & 
     622                     &            - e1u(ji  ,jj+1) * pu_mm(ji,jj+1,jk) + e1u(ji,jj) * pu_mm(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj) 
    617623               END DO 
    618624            END DO 
    … …  
    620626            DO jj = 1, jpjm1 
    621627               DO ji = 1, fs_jpim1   ! vector opt. 
    622                   zwz(ji,jj,jk) = (   ( pvn(ji+1,jj,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
    623                      &              - ( pun(ji,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
     628                  zwz(ji,jj,jk) = (   ( pv_mm(ji+1,jj,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     629                     &              - ( pu_mm(ji,jj+1,jk) + pu_mm(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
    624630               END DO 
    625631            END DO 
    … …  
    627633            DO jj = 1, jpjm1 
    628634               DO ji = 1, fs_jpim1   ! vector opt. 
    629                   zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      & 
    630                      &                              - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  )   & 
     635                  zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv_mm(ji+1,jj,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)      & 
     636                     &                              - e1u(ji  ,jj+1) * pu_mm(ji,jj+1,jk) + e1u(ji,jj) * pu_mm(ji,jj,jk)  )   & 
    631637                     &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
    632638               END DO 
    … …  
    635641            DO jj = 1, jpjm1 
    636642               DO ji = 1, fs_jpim1   ! vector opt. 
    637                   zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
    638                      &                            - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
     643                  zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pv_mm(ji+1,jj  ,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     644                     &                            - ( pu_mm(ji  ,jj+1,jk) + pu_mm(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
    639645               END DO 
    640646            END DO 
    … …  
    657663         ! 
    658664         !                                   !==  horizontal fluxes  ==! 
    659          zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
    660          zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
     665         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu_mm(:,:,jk) 
     666         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv_mm(:,:,jk) 
    661667 
    662668         !                                   !==  compute and add the vorticity term trend  =! 
    … …  
    683689               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   & 
    684690                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    685                pua(ji,jj,jk) = pua(ji,jj,jk) + zua 
    686                pva(ji,jj,jk) = pva(ji,jj,jk) + zva 
     691               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua 
     692               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva 
    687693            END DO   
    688694         END DO   
    … …  
    694700 
    695701 
    696    SUBROUTINE vor_eeT( kt, kvor, pun, pvn, pua, pva ) 
     702   SUBROUTINE vor_eeT( kt, Kmm, kvor, pu_mm, pv_mm, pu_rhs, pv_rhs ) 
    697703      !!---------------------------------------------------------------------- 
    698704      !!                ***  ROUTINE vor_eeT  *** 
    … …  
    705711      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een). 
    706712      !!      The change consists in  
    707       !!      Add this trend to the general momentum trend (ua,va). 
    708       !! 
    709       !! ** Action : - Update (ua,va) with the now vorticity term trend 
     713      !!      Add this trend to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 
     714      !! 
     715      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend 
    710716      !! 
    711717      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 
    712718      !!---------------------------------------------------------------------- 
    713719      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index 
     720      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index 
    714721      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
    715       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities 
    716       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend 
     722      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_mm, pv_mm    ! now velocities 
     723      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend 
    717724      ! 
    718725      INTEGER  ::   ji, jj, jk     ! dummy loop indices 
    … …  
    746753            DO jj = 1, jpjm1 
    747754               DO ji = 1, fs_jpim1   ! vector opt. 
    748                   zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
    749                      &             - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
     755                  zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pv_mm(ji+1,jj  ,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)    & 
     756                     &             - e1u(ji  ,jj+1) * pu_mm(ji  ,jj+1,jk) + e1u(ji,jj) * pu_mm(ji,jj,jk)  ) & 
    750757                     &          * r1_e1e2f(ji,jj) 
    751758               END DO 
    … …  
    754761            DO jj = 1, jpjm1 
    755762               DO ji = 1, fs_jpim1   ! vector opt. 
    756                   zwz(ji,jj,jk) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
    757                      &          - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
     763                  zwz(ji,jj,jk) = ( pv_mm(ji+1,jj  ,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     764                     &          - ( pu_mm(ji  ,jj+1,jk) + pu_mm(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    758765               END DO 
    759766            END DO 
    … …  
    761768            DO jj = 1, jpjm1 
    762769               DO ji = 1, fs_jpim1   ! vector opt. 
    763                   zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
    764                      &                              - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
     770                  zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv_mm(ji+1,jj  ,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)    & 
     771                     &                              - e1u(ji  ,jj+1) * pu_mm(ji  ,jj+1,jk) + e1u(ji,jj) * pu_mm(ji,jj,jk)  ) & 
    765772                     &                         * r1_e1e2f(ji,jj)    ) 
    766773               END DO 
    … …  
    769776            DO jj = 1, jpjm1 
    770777               DO ji = 1, fs_jpim1   ! vector opt. 
    771                   zwz(ji,jj,jk) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
    772                      &                        - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
     778                  zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv_mm(ji+1,jj  ,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     779                     &                        - ( pu_mm(ji  ,jj+1,jk) + pu_mm(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    773780               END DO 
    774781            END DO 
    … …  
    791798 
    792799      !                                   !==  horizontal fluxes  ==! 
    793          zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
    794          zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
     800         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu_mm(:,:,jk) 
     801         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv_mm(:,:,jk) 
    795802 
    796803         !                                   !==  compute and add the vorticity term trend  =! 
    … …  
    798805         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0 
    799806         DO ji = 2, jpi          ! split in 2 parts due to vector opt. 
    800                z1_e3t = 1._wp / e3t_n(ji,jj,jk) 
     807               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 
    801808               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t 
    802809               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t 
    … …  
    806813         DO jj = 3, jpj 
    807814            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3 
    808                z1_e3t = 1._wp / e3t_n(ji,jj,jk) 
     815               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 
    809816               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t 
    810817               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t 
    … …  
    819826               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   & 
    820827                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    821                pua(ji,jj,jk) = pua(ji,jj,jk) + zua 
    822                pva(ji,jj,jk) = pva(ji,jj,jk) + zva 
     828               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua 
     829               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva 
    823830            END DO   
    824831         END DO   

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzad.F90

    r10874 r10877  
    3636CONTAINS 
    3737 
    38    SUBROUTINE dyn_zad ( kt ) 
     38   SUBROUTINE dyn_zad ( kt, Kmm, puu, pvv, Krhs ) 
    3939      !!---------------------------------------------------------------------- 
    4040      !!                  ***  ROUTINE dynzad  *** 
    … …  
    4444      !! 
    4545      !! ** Method  :   The now vertical advection of momentum is given by: 
    46       !!         w dz(u) = ua + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*wn) dk(un) ] 
    47       !!         w dz(v) = va + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*wn) dk(vn) ] 
    48       !!      Add this trend to the general trend (ua,va): 
    49       !!         (ua,va) = (ua,va) + w dz(u,v) 
     46      !!         w dz(u) = puu(:,:,:,Krhs) + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*ww) dk(puu(:,:,:,Kmm)) ] 
     47      !!         w dz(v) = pvv(:,:,:,Krhs) + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*ww) dk(pvv(:,:,:,Kmm)) ] 
     48      !!      Add this trend to the general trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)): 
     49      !!         (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) = (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) + w dz(u,v) 
    5050      !! 
    51       !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends 
     51      !! ** Action  : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the vert. momentum adv. trends 
    5252      !!              - Send the trends to trddyn for diagnostics (l_trddyn=T) 
    5353      !!---------------------------------------------------------------------- 
    54       INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx 
     54      INTEGER                             , INTENT( in )  ::  kt               ! ocean time-step inedx 
     55      INTEGER                             , INTENT( in )  ::  Kmm, Krhs        ! ocean time level indices 
     56      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation 
    5557      ! 
    5658      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    … …  
    6870      ENDIF 
    6971 
    70       IF( l_trddyn )   THEN         ! Save ua and va trends 
     72      IF( l_trddyn )   THEN         ! Save puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends 
    7173         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )  
    72          ztrdu(:,:,:) = ua(:,:,:)  
    73          ztrdv(:,:,:) = va(:,:,:)  
     74         ztrdu(:,:,:) = puu(:,:,:,Krhs)  
     75         ztrdv(:,:,:) = pvv(:,:,:,Krhs)  
    7476      ENDIF 
    7577       
    … …  
    7779         DO jj = 2, jpj                   ! vertical fluxes  
    7880            DO ji = fs_2, jpi             ! vector opt. 
    79                zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
     81               zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk) 
    8082            END DO 
    8183         END DO 
    8284         DO jj = 2, jpjm1                 ! vertical momentum advection at w-point 
    8385            DO ji = fs_2, fs_jpim1        ! vector opt. 
    84                zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) 
    85                zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) 
     86               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( puu(ji,jj,jk-1,Kmm) - puu(ji,jj,jk,Kmm) ) 
     87               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( pvv(ji,jj,jk-1,Kmm) - pvv(ji,jj,jk,Kmm) ) 
    8688            END DO   
    8789         END DO    
    … …  
    101103         DO jj = 2, jpjm1 
    102104            DO ji = fs_2, fs_jpim1       ! vector opt. 
    103                ua(ji,jj,jk) = ua(ji,jj,jk) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
    104                va(ji,jj,jk) = va(ji,jj,jk) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
     105               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
     106               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
    105107            END DO   
    106108         END DO   
    … …  
    108110 
    109111      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic 
    110          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    111          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     112         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 
     113         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 
    112114         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt ) 
    113115         DEALLOCATE( ztrdu, ztrdv )  

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90

    r10876 r10877  
    188188               &         CALL Agrif_Sponge_dyn        ! momentum sponge 
    189189#endif 
    190                          CALL dyn_adv       ( kstp )  ! advection (vector or flux form) 
    191                          CALL dyn_vor       ( kstp )  ! vorticity term including Coriolis 
     190                         CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)   ==> RHS 
     191                         CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS 
    192192                         CALL dyn_ldf       ( kstp )  ! lateral mixing 
    193193      IF( ln_zdfosm  )   CALL dyn_osm       ( kstp )  ! OSMOSIS non-local velocity fluxes 
Note: See TracChangeset for help on using the changeset viewer.