Ignore:
Timestamp:
2019-03-21T16:15:22+01:00 (2 years ago)
Author:
davestorkey
Message:

branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps: Convert first batch of DYN routines and "wn" → "ww".

File:
1 edited

Legend:

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

    r10425 r10789  
    9696CONTAINS 
    9797 
    98    SUBROUTINE dyn_vor( kt ) 
     98   SUBROUTINE dyn_vor( kt, ktlev, pu_rhs, pv_rhs ) 
    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 (pu_rhs,pv_rhs) 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 ) ::   ktlev   ! time level index for source terms 
     110      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends 
    109111      ! 
    110112      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
     
    120122         ztrdv(:,:,:) = va(:,:,:) 
    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, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! enstrophy conserving scheme 
     125            IF( ln_stcor )            CALL vor_ens( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend 
     126         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! energy conserving scheme 
     127            IF( ln_stcor )            CALL vor_ene( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend 
     128         CASE( np_ENT )           ;   CALL vor_enT( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! energy conserving scheme (T-pts) 
     129            IF( ln_stcor )            CALL vor_enT( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend 
     130         CASE( np_EET )           ;   CALL vor_eeT( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! energy conserving scheme (een with e3t) 
     131            IF( ln_stcor )            CALL vor_eeT( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend 
     132         CASE( np_EEN )           ;   CALL vor_een( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! energy & enstrophy scheme 
     133            IF( ln_stcor )            CALL vor_een( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend 
    132134         END SELECT 
    133135         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     
    139141            ztrdv(:,:,:) = va(:,:,:) 
    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, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! energy conserving scheme (T-pts) 
     144            CASE( np_EET )           ;   CALL vor_eeT( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! energy conserving scheme (een with e3t) 
     145            CASE( np_ENE )           ;   CALL vor_ene( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! energy conserving scheme 
     146            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! enstrophy conserving scheme 
     147            CASE( np_EEN )           ;   CALL vor_een( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! energy & enstrophy scheme 
    146148            END SELECT 
    147149            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     
    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, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! total vorticity trend 
     161            IF( ln_stcor )   CALL vor_enT( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! 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, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! total vorticity trend 
     164            IF( ln_stcor )   CALL vor_eeT( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! 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, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! total vorticity trend 
     167            IF( ln_stcor )   CALL vor_ene( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! 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, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! total vorticity trend 
     170            IF( ln_stcor )   CALL vor_ens( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )  ! 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, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! relative vorticity or metric trend (ens) 
     173                             CALL vor_ene( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! planetary vorticity trend (ene) 
     174            IF( ln_stcor )   CALL vor_ene( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! 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, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! total vorticity trend 
     177            IF( ln_stcor )   CALL vor_een( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! 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, ktlev, 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 (u_rhs,v_rhs) with the now vorticity term trend 
    206208      !!---------------------------------------------------------------------- 
    207209      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index 
     210      INTEGER                         , INTENT( in )  ::   ktlev            ! time level index for source terms 
    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,ktlev) 
    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,ktlev) 
    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,ktlev) 
    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,ktlev) 
    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,ktlev) 
    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,ktlev)                    & 
    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,ktlev)                    & 
    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, ktlev, kvor, pu, pv, pu_rhs, pva_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 vv(:,:,:,ktlev)) ] 
     340      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u uu(:,:,:,ktlev)) ] 
    338341      !!       where rvor is the relative vorticity 
    339342      !! 
    340       !! ** Action : - Update (ua,va) with the now vorticity term trend 
     343      !! ** Action : - Update (u_rhs,v_rhs) 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 )  ::   ktlev            ! time level index for source terms 
    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, pv    ! now velocities 
     351      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_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(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    & 
     375                     &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(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(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     382                     &       - ( pu(ji  ,jj+1,jk) + pu(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(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      & 
     389                     &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(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(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     396                     &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    393397               END DO 
    394398            END DO 
     
    407411         IF( ln_sco ) THEN 
    408412            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
    409             zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
    410             zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
     413            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,ktlev) * pu(:,:,jk) 
     414            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,ktlev) * pv(:,:,jk) 
    411415         ELSE 
    412             zwx(:,:) = e2u(:,:) * pun(:,:,jk) 
    413             zwy(:,:) = e1v(:,:) * pvn(:,:,jk) 
     416            zwx(:,:) = e2u(:,:) * pu(:,:,jk) 
     417            zwy(:,:) = e1v(:,:) * pv(:,:,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               pva_rhs(ji,jj,jk) = pva_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, ktlev, kvor, pu, pv, pu_rhs, pva_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 vv(:,:,:,ktlev)) ] 
     448      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u uu(:,:,:,ktlev)) ] 
     449      !!      Add this trend to the general momentum trend (u_rhs,v_rhs): 
     450      !!          (u_rhs,v_rhs) = (u_rhs,v_rhs) + ( voru , vorv ) 
     451      !! 
     452      !! ** Action : - Update (u_rhs,v_rhs) 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 )  ::   ktlev            ! time level index for source terms 
    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, pv    ! now velocities 
     460      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_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(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    & 
     483                     &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(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(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     490                     &       - ( pu(ji  ,jj+1,jk) + pu(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(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)  & 
     497                     &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(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(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     504                     &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    500505               END DO 
    501506            END DO 
     
    514519         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==! 
    515520            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
    516             zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
    517             zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
     521            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,ktlev) * pu(:,:,jk) 
     522            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,ktlev) * pv(:,:,jk) 
    518523         ELSE 
    519             zwx(:,:) = e2u(:,:) * pun(:,:,jk) 
    520             zwy(:,:) = e1v(:,:) * pvn(:,:,jk) 
     524            zwx(:,:) = e2u(:,:) * pu(:,:,jk) 
     525            zwy(:,:) = e1v(:,:) * pv(:,:,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               pva_rhs(ji,jj,jk) = pva_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, ktlev, kvor, pu, pv, pu_rhs, pva_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 (u_rhs,v_rhs). 
     556      !! 
     557      !! ** Action : - Update (u_rhs,v_rhs) 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 )  ::   ktlev            ! time level index for source terms 
    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, pv    ! now velocities 
     565      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_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,ktlev)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,ktlev)*tmask(ji+1,jj+1,jk)   & 
     591                     &    + e3t(ji,jj  ,jk,ktlev)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,ktlev)*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,ktlev)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,ktlev)*tmask(ji+1,jj+1,jk)   & 
     601                     &    + e3t(ji,jj  ,jk,ktlev)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,ktlev)*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(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  & 
     622                     &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(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(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     629                     &              - ( pu(ji,jj+1,jk) + pu(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(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      & 
     636                     &                              - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(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(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     644                     &                            - ( pu(ji  ,jj+1,jk) + pu(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,ktlev) * pu(:,:,jk) 
     666         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,ktlev) * pv(:,:,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               pva_rhs(ji,jj,jk) = pva_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, ktlev, kvor, pu, pv, pu_rhs, pva_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 (u_rhs,v_rhs). 
     714      !! 
     715      !! ** Action : - Update (u_rhs,v_rhs) 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 )  ::   ktlev            ! time level index for source terms 
    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, pv    ! now velocities 
     723      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_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(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    & 
     756                     &             - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(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(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     764                     &          - ( pu(ji  ,jj+1,jk) + pu(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(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    & 
     771                     &                              - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(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(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     779                     &                        - ( pu(ji  ,jj+1,jk) + pu(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,ktlev) * pu(:,:,jk) 
     801         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,ktlev) * pv(:,:,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,ktlev) 
    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,ktlev) 
    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               pva_rhs(ji,jj,jk) = pva_rhs(ji,jj,jk) + zva 
    823830            END DO   
    824831         END DO   
     
    866873         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f 
    867874         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix 
    868          WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk 
     875         WRITE(numout,*) '      masked (=T) or uu(:,:,:,ktlev)masked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk 
    869876      ENDIF 
    870877 
    871878      IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available') 
    872879 
    873 !!gm  this should be removed when choosing a unique strategy for fmask at the coast 
     880!!gm  this should be removed when choosing a uu(:,:,:,ktlev)ique strategy for fmask at the coast 
    874881      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks 
    875882      ! at angles with three ocean points and one land point 
Note: See TracChangeset for help on using the changeset viewer.