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.
dynvor.F90 in NEMO/branches/2020/dev_r13898_Tiling_Cleanup_MPI3/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r13898_Tiling_Cleanup_MPI3/src/OCE/DYN/dynvor.F90 @ 13906

Last change on this file since 13906 was 13906, checked in by mocavero, 3 years ago

Merge with dev_r13296_HPC-07_mocavero_mpi3
  • Property svn:keywords set to Id
File size: 54.0 KB
Line 
1MODULE dynvor
2   !!======================================================================
3   !!                       ***  MODULE  dynvor  ***
4   !! Ocean dynamics: Update the momentum trend with the relative and
5   !!                 planetary vorticity trends
6   !!======================================================================
7   !! History :  OPA  ! 1989-12  (P. Andrich)  vor_ens: Original code
8   !!            5.0  ! 1991-11  (G. Madec)  vor_ene, vor_mix: Original code
9   !!            6.0  ! 1996-01  (G. Madec)  s-coord, suppress work arrays
10   !!   NEMO     0.5  ! 2002-08  (G. Madec)  F90: Free form and module
11   !!            1.0  ! 2004-02  (G. Madec)  vor_een: Original code
12   !!             -   ! 2003-08  (G. Madec)  add vor_ctl
13   !!             -   ! 2005-11  (G. Madec)  add dyn_vor (new step architecture)
14   !!            2.0  ! 2006-11  (G. Madec)  flux form advection: add metric term
15   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
16   !!            3.3  ! 2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase
17   !!            3.7  ! 2014-04  (G. Madec)  trend simplification: suppress jpdyn_trd_dat vorticity
18   !!             -   ! 2014-06  (G. Madec)  suppression of velocity curl from in-core memory
19   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
20   !!            4.0  ! 2017-07  (G. Madec)  linear dynamics + trends diag. with Stokes-Coriolis
21   !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet)
22   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation
23   !!----------------------------------------------------------------------
24
25   !!----------------------------------------------------------------------
26   !!   dyn_vor       : Update the momentum trend with the vorticity trend
27   !!       vor_ens   : enstrophy conserving scheme       (ln_dynvor_ens=T)
28   !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T)
29   !!       vor_een   : energy and enstrophy conserving   (ln_dynvor_een=T)
30   !!   dyn_vor_init  : set and control of the different vorticity option
31   !!----------------------------------------------------------------------
32   USE oce            ! ocean dynamics and tracers
33   USE dom_oce        ! ocean space and time domain
34   USE dommsk         ! ocean mask
35   USE dynadv         ! momentum advection
36   USE trd_oce        ! trends: ocean variables
37   USE trddyn         ! trend manager: dynamics
38   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
39   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force
40   !
41   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
42   USE prtctl         ! Print control
43   USE in_out_manager ! I/O manager
44   USE lib_mpp        ! MPP library
45   USE timing         ! Timing
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC   dyn_vor        ! routine called by step.F90
51   PUBLIC   dyn_vor_init   ! routine called by nemogcm.F90
52
53   !                                   !!* Namelist namdyn_vor: vorticity term
54   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme          (ENS)
55   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: f-point energy conserving scheme     (ENE)
56   LOGICAL, PUBLIC ::   ln_dynvor_enT   !: t-point energy conserving scheme     (ENT)
57   LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET)
58   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN)
59   INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
60   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX)
61   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
62
63   INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme
64   !                                       ! associated indices:
65   INTEGER, PUBLIC, PARAMETER ::   np_ENS = 0   ! ENS scheme
66   INTEGER, PUBLIC, PARAMETER ::   np_ENE = 1   ! ENE scheme
67   INTEGER, PUBLIC, PARAMETER ::   np_ENT = 2   ! ENT scheme (t-point vorticity)
68   INTEGER, PUBLIC, PARAMETER ::   np_EET = 3   ! EET scheme (EEN using e3t)
69   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme
70   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme
71
72   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity
73   !                               ! associated indices:
74   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
75   INTEGER, PUBLIC, PARAMETER ::   np_RVO = 2         ! relative vorticity
76   INTEGER, PUBLIC, PARAMETER ::   np_MET = 3         ! metric term
77   INTEGER, PUBLIC, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity)
78   INTEGER, PUBLIC, PARAMETER ::   np_CME = 5         ! Coriolis + metric term
79
80   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation
81   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -
82   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2v)/(2*e1e2f)  used in F-point metric term calculation
83   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1u)/(2*e1e2f)   -        -      -       -
84   
85   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4
86   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8
87   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12
88   
89   !! * Substitutions
90#  include "do_loop_substitute.h90"
91#  include "domzgr_substitute.h90"
92
93   !!----------------------------------------------------------------------
94   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
95   !! $Id$
96   !! Software governed by the CeCILL license (see ./LICENSE)
97   !!----------------------------------------------------------------------
98CONTAINS
99
100   SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs )
101      !!----------------------------------------------------------------------
102      !!
103      !! ** Purpose :   compute the lateral ocean tracer physics.
104      !!
105      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
106      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
107      !!               and planetary vorticity trends) and send them to trd_dyn
108      !!               for futher diagnostics (l_trddyn=T)
109      !!----------------------------------------------------------------------
110      INTEGER                             , INTENT( in  ) ::   kt          ! ocean time-step index
111      INTEGER                             , INTENT( in  ) ::   Kmm, Krhs   ! ocean time level indices
112      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv    ! ocean velocity field and RHS of momentum equation
113      !
114      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv
115      !!----------------------------------------------------------------------
116      !
117      IF( ln_timing )   CALL timing_start('dyn_vor')
118      !
119      IF( l_trddyn ) THEN     !==  trend diagnostics case : split the added trend in two parts  ==!
120         !
121         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
122         !
123         ztrdu(:,:,:) = puu(:,:,:,Krhs)            !* planetary vorticity trend (including Stokes-Coriolis force)
124         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
125         SELECT CASE( nvor_scheme )
126         CASE( np_ENS )           ;   CALL vor_ens( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! enstrophy conserving scheme
127            IF( ln_stcor )            CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
128         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme
129            IF( ln_stcor )            CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
130         CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (T-pts)
131            IF( ln_stcor )            CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
132         CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (een with e3t)
133            IF( ln_stcor )            CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
134         CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy & enstrophy scheme
135            IF( ln_stcor )            CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
136         END SELECT
137         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
138         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
139         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt, Kmm )
140         !
141         IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case)
142            ztrdu(:,:,:) = puu(:,:,:,Krhs)
143            ztrdv(:,:,:) = pvv(:,:,:,Krhs)
144            SELECT CASE( nvor_scheme )
145            CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (T-pts)
146            CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (een with e3t)
147            CASE( np_ENE )           ;   CALL vor_ene( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme
148            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! enstrophy conserving scheme
149            CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy & enstrophy scheme
150            END SELECT
151            ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
152            ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
153            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt, Kmm )
154         ENDIF
155         !
156         DEALLOCATE( ztrdu, ztrdv )
157         !
158      ELSE              !==  total vorticity trend added to the general trend  ==!
159         !
160         SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==!
161         CASE( np_ENT )                        !* energy conserving scheme  (T-pts)
162                             CALL vor_enT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
163            IF( ln_stcor )   CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
164         CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t)
165                             CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
166            IF( ln_stcor )   CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
167         CASE( np_ENE )                        !* energy conserving scheme
168                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
169            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
170         CASE( np_ENS )                        !* enstrophy conserving scheme
171                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend
172            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
173         CASE( np_MIX )                        !* mixed ene-ens scheme
174                             CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! relative vorticity or metric trend (ens)
175                             CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! planetary vorticity trend (ene)
176            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
177         CASE( np_EEN )                        !* energy and enstrophy conserving scheme
178                             CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
179            IF( ln_stcor )   CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
180         END SELECT
181         !
182      ENDIF
183      !
184      !                       ! print sum trends (used for debugging)
185      IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor  - Ua: ', mask1=umask,               &
186         &                                tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
187      !
188      IF( ln_timing )   CALL timing_stop('dyn_vor')
189      !
190   END SUBROUTINE dyn_vor
191
192
193   SUBROUTINE vor_enT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
194      !!----------------------------------------------------------------------
195      !!                  ***  ROUTINE vor_enT  ***
196      !!
197      !! ** Purpose :   Compute the now total vorticity trend and add it to
198      !!      the general trend of the momentum equation.
199      !!
200      !! ** Method  :   Trend evaluated using now fields (centered in time)
201      !!       and t-point evaluation of vorticity (planetary and relative).
202      !!       conserves the horizontal kinetic energy.
203      !!         The general trend of momentum is increased due to the vorticity
204      !!       term which is given by:
205      !!          voru = 1/bu  mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t  mj[vn] ]
206      !!          vorv = 1/bv  mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f  mj[un] ]
207      !!       where rvor is the relative vorticity at f-point
208      !!
209      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
210      !!----------------------------------------------------------------------
211      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
212      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
213      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric
214      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities
215      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend
216      !
217      INTEGER  ::   ji, jj, jk           ! dummy loop indices
218      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
219      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx, zwy, zwt   ! 2D workspace
220      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz             ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined
221      !!----------------------------------------------------------------------
222      !
223      IF( kt == nit000 ) THEN
224         IF(lwp) WRITE(numout,*)
225         IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
226         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
227      ENDIF
228      !
229      !
230      SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
231      CASE ( np_RVO )                           !* relative vorticity
232         DO jk = 1, jpkm1                                 ! Horizontal slab
233            DO_2D( 1, 0, 1, 0 )
234               zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
235                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
236            END_2D
237            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
238               DO_2D( 1, 0, 1, 0 )
239                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
240               END_2D
241            ENDIF
242         END DO
243
244#if defined key_mpi3
245         CALL lbc_lnk_nc_multi( 'dynvor', zwz, 'F', 1.0_wp )
246#else
247         CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
248#endif
249
250      CASE ( np_CRV )                           !* Coriolis + relative vorticity
251         DO jk = 1, jpkm1                                 ! Horizontal slab
252            DO_2D( 1, 0, 1, 0 )                          ! relative vorticity
253               zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   &
254                  &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj)
255            END_2D
256            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
257               DO_2D( 1, 0, 1, 0 )
258                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
259               END_2D
260            ENDIF
261         END DO
262
263#if defined key_mpi3
264         CALL lbc_lnk_nc_multi( 'dynvor', zwz, 'F', 1.0_wp )
265#else
266         CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
267#endif
268
269      END SELECT
270
271      !                                                ! ===============
272      DO jk = 1, jpkm1                                 ! Horizontal slab
273      !                                                ! ===============
274
275         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
276         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
277            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm)
278         CASE ( np_RVO )                           !* relative vorticity
279            DO_2D( 0, 1, 0, 1 )
280               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   &
281                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) &
282                  &                  * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
283            END_2D
284         CASE ( np_MET )                           !* metric term
285            DO_2D( 0, 1, 0, 1 )
286               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   &
287                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) &
288                  &             * e3t(ji,jj,jk,Kmm)
289            END_2D
290         CASE ( np_CRV )                           !* Coriolis + relative vorticity
291            DO_2D( 0, 1, 0, 1 )
292               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    &
293                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) &
294                  &                                 * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
295            END_2D
296         CASE ( np_CME )                           !* Coriolis + metric
297            DO_2D( 0, 1, 0, 1 )
298               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           &
299                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  &
300                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) &
301                    &          * e3t(ji,jj,jk,Kmm)
302            END_2D
303         CASE DEFAULT                                             ! error
304            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
305         END SELECT
306         !
307         !                                   !==  compute and add the vorticity term trend  =!
308         DO_2D( 0, 0, 0, 0 )
309            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                    &
310               &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   &
311               &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   )
312               !
313            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    &
314               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
315               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
316         END_2D
317         !                                             ! ===============
318      END DO                                           !   End of slab
319      !                                                ! ===============
320   END SUBROUTINE vor_enT
321
322
323   SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
324      !!----------------------------------------------------------------------
325      !!                  ***  ROUTINE vor_ene  ***
326      !!
327      !! ** Purpose :   Compute the now total vorticity trend and add it to
328      !!      the general trend of the momentum equation.
329      !!
330      !! ** Method  :   Trend evaluated using now fields (centered in time)
331      !!       and the Sadourny (1975) flux form formulation : conserves the
332      !!       horizontal kinetic energy.
333      !!         The general trend of momentum is increased due to the vorticity
334      !!       term which is given by:
335      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v pvv(:,:,:,Kmm)) ]
336      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u puu(:,:,:,Kmm)) ]
337      !!       where rvor is the relative vorticity
338      !!
339      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
340      !!
341      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
342      !!----------------------------------------------------------------------
343      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
344      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
345      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
346      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
347      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
348      !
349      INTEGER  ::   ji, jj, jk           ! dummy loop indices
350      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
351      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace
352      !!----------------------------------------------------------------------
353      !
354      IF( kt == nit000 ) THEN
355         IF(lwp) WRITE(numout,*)
356         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
357         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
358      ENDIF
359      !
360      !                                                ! ===============
361      DO jk = 1, jpkm1                                 ! Horizontal slab
362         !                                             ! ===============
363         !
364         SELECT CASE( kvor )                 !==  vorticity considered  ==!
365         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
366            zwz(:,:) = ff_f(:,:) 
367         CASE ( np_RVO )                           !* relative vorticity
368            DO_2D( 1, 0, 1, 0 )
369               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
370                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
371            END_2D
372         CASE ( np_MET )                           !* metric term
373            DO_2D( 1, 0, 1, 0 )
374               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
375                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
376            END_2D
377         CASE ( np_CRV )                           !* Coriolis + relative vorticity
378            DO_2D( 1, 0, 1, 0 )
379               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
380                  &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
381            END_2D
382         CASE ( np_CME )                           !* Coriolis + metric
383            DO_2D( 1, 0, 1, 0 )
384               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
385                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
386            END_2D
387         CASE DEFAULT                                             ! error
388            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
389         END SELECT
390         !
391         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
392            DO_2D( 1, 0, 1, 0 )
393               zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
394            END_2D
395         ENDIF
396
397         IF( ln_sco ) THEN
398            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
399            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
400            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
401         ELSE
402            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
403            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
404         ENDIF
405         !                                   !==  compute and add the vorticity term trend  =!
406         DO_2D( 0, 0, 0, 0 )
407            zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
408            zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
409            zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
410            zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
411            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 )
412            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 ) 
413         END_2D
414         !                                             ! ===============
415      END DO                                           !   End of slab
416      !                                                ! ===============
417   END SUBROUTINE vor_ene
418
419
420   SUBROUTINE vor_ens( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
421      !!----------------------------------------------------------------------
422      !!                ***  ROUTINE vor_ens  ***
423      !!
424      !! ** Purpose :   Compute the now total vorticity trend and add it to
425      !!      the general trend of the momentum equation.
426      !!
427      !! ** Method  :   Trend evaluated using now fields (centered in time)
428      !!      and the Sadourny (1975) flux FORM formulation : conserves the
429      !!      potential enstrophy of a horizontally non-divergent flow. the
430      !!      trend of the vorticity term is given by:
431      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v pvv(:,:,:,Kmm)) ]
432      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u puu(:,:,:,Kmm)) ]
433      !!      Add this trend to the general momentum trend:
434      !!          (u(rhs),v(Krhs)) = (u(rhs),v(Krhs)) + ( voru , vorv )
435      !!
436      !! ** Action : - Update (pu_rhs,pv_rhs)) arrays with the now vorticity term trend
437      !!
438      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
439      !!----------------------------------------------------------------------
440      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
441      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
442      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
443      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
444      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
445      !
446      INTEGER  ::   ji, jj, jk   ! dummy loop indices
447      REAL(wp) ::   zuav, zvau   ! local scalars
448      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! 2D workspace
449      !!----------------------------------------------------------------------
450      !
451      IF( kt == nit000 ) THEN
452         IF(lwp) WRITE(numout,*)
453         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
454         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
455      ENDIF
456      !                                                ! ===============
457      DO jk = 1, jpkm1                                 ! Horizontal slab
458         !                                             ! ===============
459         !
460         SELECT CASE( kvor )                 !==  vorticity considered  ==!
461         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
462            zwz(:,:) = ff_f(:,:) 
463         CASE ( np_RVO )                           !* relative vorticity
464            DO_2D( 1, 0, 1, 0 )
465               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
466                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
467            END_2D
468         CASE ( np_MET )                           !* metric term
469            DO_2D( 1, 0, 1, 0 )
470               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
471                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
472            END_2D
473         CASE ( np_CRV )                           !* Coriolis + relative vorticity
474            DO_2D( 1, 0, 1, 0 )
475               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
476                  &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
477            END_2D
478         CASE ( np_CME )                           !* Coriolis + metric
479            DO_2D( 1, 0, 1, 0 )
480               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
481                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
482            END_2D
483         CASE DEFAULT                                             ! error
484            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
485         END SELECT
486         !
487         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==!
488            DO_2D( 1, 0, 1, 0 )
489               zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
490            END_2D
491         ENDIF
492         !
493         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
494            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
495            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
496            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
497         ELSE
498            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
499            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
500         ENDIF
501         !                                   !==  compute and add the vorticity term trend  =!
502         DO_2D( 0, 0, 0, 0 )
503            zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
504               &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
505            zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
506               &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
507            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
508            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
509         END_2D
510         !                                             ! ===============
511      END DO                                           !   End of slab
512      !                                                ! ===============
513   END SUBROUTINE vor_ens
514
515
516   SUBROUTINE vor_een( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
517      !!----------------------------------------------------------------------
518      !!                ***  ROUTINE vor_een  ***
519      !!
520      !! ** Purpose :   Compute the now total vorticity trend and add it to
521      !!      the general trend of the momentum equation.
522      !!
523      !! ** Method  :   Trend evaluated using now fields (centered in time)
524      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
525      !!      both the horizontal kinetic energy and the potential enstrophy
526      !!      when horizontal divergence is zero (see the NEMO documentation)
527      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
528      !!
529      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
530      !!
531      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
532      !!----------------------------------------------------------------------
533      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
534      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
535      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
536      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
537      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
538      !
539      INTEGER  ::   ji, jj, jk   ! dummy loop indices
540      INTEGER  ::   ierr         ! local integer
541      REAL(wp) ::   zua, zva     ! local scalars
542      REAL(wp) ::   zmsk, ze3f   ! local scalars
543      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx , zwy , z1_e3f
544      REAL(wp), DIMENSION(jpi,jpj)       ::   ztnw, ztne, ztsw, ztse
545      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz   ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
546      !!----------------------------------------------------------------------
547      !
548      IF( kt == nit000 ) THEN
549         IF(lwp) WRITE(numout,*)
550         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
551         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
552      ENDIF
553      !
554      !                                                ! ===============
555      DO jk = 1, jpkm1                                 ! Horizontal slab
556         !                                             ! ===============
557         !
558         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
559         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
560            DO_2D( 1, 0, 1, 0 )
561               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
562                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
563                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
564                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
565               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
566               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
567               ENDIF
568            END_2D
569         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
570            DO_2D( 1, 0, 1, 0 )
571               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
572                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
573                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
574                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
575               zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
576                  &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
577               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3f
578               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
579               ENDIF
580            END_2D
581         END SELECT
582         !
583         SELECT CASE( kvor )                 !==  vorticity considered  ==!
584         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
585            DO_2D( 1, 0, 1, 0 )
586               zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
587            END_2D
588         CASE ( np_RVO )                           !* relative vorticity
589            DO_2D( 1, 0, 1, 0 )
590               zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
591                  &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
592            END_2D
593         CASE ( np_MET )                           !* metric term
594            DO_2D( 1, 0, 1, 0 )
595               zwz(ji,jj,jk) = (   ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
596                  &              - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
597            END_2D
598         CASE ( np_CRV )                           !* Coriolis + relative vorticity
599            DO_2D( 1, 0, 1, 0 )
600               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
601                  &                              - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  )   &
602                  &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
603            END_2D
604         CASE ( np_CME )                           !* Coriolis + metric
605            DO_2D( 1, 0, 1, 0 )
606               zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
607                  &                            - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
608            END_2D
609         CASE DEFAULT                                             ! error
610            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
611         END SELECT
612         !
613         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
614            DO_2D( 1, 0, 1, 0 )
615               zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
616            END_2D
617         ENDIF
618      END DO                                           !   End of slab
619         !
620#if defined key_mpi3
621      CALL lbc_lnk_nc_multi( 'dynvor', zwz, 'F', 1.0_wp )
622#else
623      CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
624#endif
625
626      DO jk = 1, jpkm1                                 ! Horizontal slab
627         !
628         !                                   !==  horizontal fluxes  ==!
629         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
630         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
631
632         !                                   !==  compute and add the vorticity term trend  =!
633         jj = 2
634         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
635         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
636               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
637               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
638               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
639               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
640         END DO
641         DO jj = 3, jpj
642            DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3
643               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
644               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
645               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
646               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
647            END DO
648         END DO
649         DO_2D( 0, 0, 0, 0 )
650            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
651               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
652            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
653               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
654            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
655            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
656         END_2D
657         !                                             ! ===============
658      END DO                                           !   End of slab
659      !                                                ! ===============
660   END SUBROUTINE vor_een
661
662
663
664   SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
665      !!----------------------------------------------------------------------
666      !!                ***  ROUTINE vor_eeT  ***
667      !!
668      !! ** Purpose :   Compute the now total vorticity trend and add it to
669      !!      the general trend of the momentum equation.
670      !!
671      !! ** Method  :   Trend evaluated using now fields (centered in time)
672      !!      and the Arakawa and Lamb (1980) vector form formulation using
673      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
674      !!      The change consists in
675      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
676      !!
677      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
678      !!
679      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
680      !!----------------------------------------------------------------------
681      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
682      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
683      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
684      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
685      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
686      !
687      INTEGER  ::   ji, jj, jk     ! dummy loop indices
688      INTEGER  ::   ierr           ! local integer
689      REAL(wp) ::   zua, zva       ! local scalars
690      REAL(wp) ::   zmsk, z1_e3t   ! local scalars
691      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx , zwy 
692      REAL(wp), DIMENSION(jpi,jpj)       ::   ztnw, ztne, ztsw, ztse
693      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz   ! 3D workspace, avoid lbc_lnk on jpk that is not defined
694      !!----------------------------------------------------------------------
695      !
696      IF( kt == nit000 ) THEN
697         IF(lwp) WRITE(numout,*)
698         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
699         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
700      ENDIF
701      !
702      !                                                ! ===============
703      DO jk = 1, jpkm1                                 ! Horizontal slab
704         !                                             ! ===============
705         !
706         !
707         SELECT CASE( kvor )                 !==  vorticity considered  ==!
708         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
709            DO_2D( 1, 0, 1, 0 )
710               zwz(ji,jj,jk) = ff_f(ji,jj)
711            END_2D
712         CASE ( np_RVO )                           !* relative vorticity
713            DO_2D( 1, 0, 1, 0 )
714               zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
715                  &             - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
716                  &          * r1_e1e2f(ji,jj)
717            END_2D
718         CASE ( np_MET )                           !* metric term
719            DO_2D( 1, 0, 1, 0 )
720               zwz(ji,jj,jk) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
721                  &          - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
722            END_2D
723         CASE ( np_CRV )                           !* Coriolis + relative vorticity
724            DO_2D( 1, 0, 1, 0 )
725               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
726                  &                              - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
727                  &                         * r1_e1e2f(ji,jj)    )
728            END_2D
729         CASE ( np_CME )                           !* Coriolis + metric
730            DO_2D( 1, 0, 1, 0 )
731               zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
732                  &                        - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
733            END_2D
734         CASE DEFAULT                                             ! error
735            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
736         END SELECT
737         !
738         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
739            DO_2D( 1, 0, 1, 0 )
740               zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
741            END_2D
742         ENDIF
743      END DO
744      !
745#if defined key_mpi3
746      CALL lbc_lnk_nc_multi( 'dynvor', zwz, 'F', 1.0_wp )
747#else
748      CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
749#endif
750      !
751      DO jk = 1, jpkm1                                 ! Horizontal slab
752
753      !                                   !==  horizontal fluxes  ==!
754         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
755         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
756
757         !                                   !==  compute and add the vorticity term trend  =!
758         jj = 2
759         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
760         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
761               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
762               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
763               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
764               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
765               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
766         END DO
767         DO jj = 3, jpj
768            DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3
769               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
770               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
771               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
772               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
773               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
774            END DO
775         END DO
776         DO_2D( 0, 0, 0, 0 )
777            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
778               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
779            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
780               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
781            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
782            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
783         END_2D
784         !                                             ! ===============
785      END DO                                           !   End of slab
786      !                                                ! ===============
787   END SUBROUTINE vor_eeT
788
789
790   SUBROUTINE dyn_vor_init
791      !!---------------------------------------------------------------------
792      !!                  ***  ROUTINE dyn_vor_init  ***
793      !!
794      !! ** Purpose :   Control the consistency between cpp options for
795      !!              tracer advection schemes
796      !!----------------------------------------------------------------------
797      INTEGER ::   ji, jj, jk    ! dummy loop indices
798      INTEGER ::   ioptio, ios   ! local integer
799      !!
800      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   &
801         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk
802      !!----------------------------------------------------------------------
803      !
804      IF(lwp) THEN
805         WRITE(numout,*)
806         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
807         WRITE(numout,*) '~~~~~~~~~~~~'
808      ENDIF
809      !
810      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
811901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' )
812      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
813902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' )
814      IF(lwm) WRITE ( numond, namdyn_vor )
815      !
816      IF(lwp) THEN                    ! Namelist print
817         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
818         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
819         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene
820         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT
821         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT
822         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
823         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
824         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
825         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
826      ENDIF
827
828      IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available')
829
830!!gm  this should be removed when choosing a unique strategy for fmask at the coast
831      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
832      ! at angles with three ocean points and one land point
833      IF(lwp) WRITE(numout,*)
834      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
835      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
836         DO_3D( 1, 0, 1, 0, 1, jpk )
837            IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
838               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
839         END_3D
840         !
841#if defined key_mpi3
842         CALL lbc_lnk_nc_multi( 'dynvor', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
843#else
844         CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
845#endif
846         !
847      ENDIF
848!!gm end
849
850      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
851      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF
852      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
853      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF
854      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF
855      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF
856      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF
857      !
858      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
859      !                     
860      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
861      ncor = np_COR                       ! planetary vorticity
862      SELECT CASE( n_dynadv )
863      CASE( np_LIN_dyn )
864         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
865         nrvm = np_COR        ! planetary vorticity
866         ntot = np_COR        !     -         -
867      CASE( np_VEC_c2  )
868         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
869         nrvm = np_RVO        ! relative vorticity
870         ntot = np_CRV        ! relative + planetary vorticity         
871      CASE( np_FLX_c2 , np_FLX_ubs  )
872         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
873         nrvm = np_MET        ! metric term
874         ntot = np_CME        ! Coriolis + metric term
875         !
876         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term:
877         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
878            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
879            DO_2D( 0, 0, 0, 0 )
880               di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
881               dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
882            END_2D
883#if defined key_mpi3
884            CALL lbc_lnk_nc_multi( 'dynvor', di_e2u_2, 'T', -1.0_wp , dj_e1v_2, 'T', -1.0_wp )   ! Lateral boundary conditions
885#else
886            CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1.0_wp , dj_e1v_2, 'T', -1.0_wp )   ! Lateral boundary conditions
887#endif
888            !
889         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
890            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
891            DO_2D( 1, 0, 1, 0 )
892               di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
893               dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
894            END_2D
895#if defined key_mpi3
896            CALL lbc_lnk_nc_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1.0_wp , dj_e1u_2e1e2f, 'F', -1.0_wp )   ! Lateral boundary conditions
897#else
898            CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1.0_wp , dj_e1u_2e1e2f, 'F', -1.0_wp )   ! Lateral boundary conditions
899#endif
900         END SELECT
901         !
902      END SELECT
903     
904      IF(lwp) THEN                   ! Print the choice
905         WRITE(numout,*)
906         SELECT CASE( nvor_scheme )
907         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
908         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
909         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
910         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
911         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
912         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
913         END SELECT         
914      ENDIF
915      !
916   END SUBROUTINE dyn_vor_init
917
918   !!==============================================================================
919END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.