source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/penzAM98/MY_SRC/dynvor.F90 @ 13562

Last change on this file since 13562 was 13562, checked in by gm, 5 months ago

zgr_pse created only for NS coast

File size: 71.9 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   !!            4.x  ! 2020-03  (G. Madec, A. Nasser)  make ln_dynvor_msk truly efficient on relative vorticity
24   !!            4.x  ! 2020-03  (G. Madec, A. Nasser)  alternate direction computation of vorticity tendancy
25   !!                 !                                 for ENS, ENE
26   !!----------------------------------------------------------------------
27
28   !!----------------------------------------------------------------------
29   !!   dyn_vor       : Update the momentum trend with the vorticity trend
30   !!       vor_ens   : enstrophy conserving scheme       (ln_dynvor_ens=T)
31   !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T)
32   !!       vor_een   : energy and enstrophy conserving   (ln_dynvor_een=T)
33   !!   dyn_vor_init  : set and control of the different vorticity option
34   !!----------------------------------------------------------------------
35   USE oce            ! ocean dynamics and tracers
36   USE dom_oce        ! ocean space and time domain
37   USE dommsk         ! ocean mask
38   USE dynadv         ! momentum advection
39   USE trd_oce        ! trends: ocean variables
40   USE trddyn         ! trend manager: dynamics
41   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
42   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force
43   !
44   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
45   USE prtctl         ! Print control
46   USE in_out_manager ! I/O manager
47   USE lib_mpp        ! MPP library
48   USE timing         ! Timing
49!!anAD only
50   USE dynkeg, ONLY : dyn_kegAD
51   USE dynadv, ONLY : nn_dynkeg
52
53   IMPLICIT NONE
54   PRIVATE
55
56   PUBLIC   dyn_vor        ! routine called by step.F90
57   PUBLIC   dyn_vor_init   ! routine called by nemogcm.F90
58
59   !                                   !!* Namelist namdyn_vor: vorticity term
60   LOGICAL, PUBLIC ::   ln_dynvor_ens      !: enstrophy conserving scheme          (ENS)
61   LOGICAL, PUBLIC ::   ln_dynvor_ens_adVO   = .FALSE.   !: AD enstrophy conserving scheme       (ENS_ad)
62   LOGICAL, PUBLIC ::   ln_dynvor_ens_adKE   = .FALSE.   !: AD enstrophy conserving scheme       (ENS_ad)
63   LOGICAL, PUBLIC ::   ln_dynvor_ens_adKEVO = .FALSE.   !: AD enstrophy conserving scheme       (ENS_ad)
64   LOGICAL, PUBLIC ::   ln_dynvor_ene      !: f-point energy conserving scheme     (ENE)
65   LOGICAL, PUBLIC ::   ln_dynvor_ene_adVO   = .FALSE.   !: f-point AD energy conserving scheme  (ENE_ad)
66   LOGICAL, PUBLIC ::   ln_dynvor_ene_adKE   = .FALSE.   !: f-point AD energy conserving scheme  (ENE_ad)
67   LOGICAL, PUBLIC ::   ln_dynvor_ene_adKEVO = .FALSE.   !: f-point AD energy conserving scheme  (ENE_ad)
68   LOGICAL, PUBLIC ::   ln_dynvor_enT      !: t-point energy conserving scheme     (ENT)
69   LOGICAL, PUBLIC ::   ln_dynvor_eeT      !: t-point energy conserving scheme     (EET)
70   LOGICAL, PUBLIC ::   ln_dynvor_een      !: energy & enstrophy conserving scheme (EEN)
71   INTEGER, PUBLIC ::      nn_een_e3f         !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
72   LOGICAL, PUBLIC ::   ln_dynvor_mix      !: mixed scheme                         (MIX)
73   LOGICAL, PUBLIC ::   ln_dynvor_msk      !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
74
75   INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme
76   !                                       ! associated indices:
77   INTEGER, PUBLIC, PARAMETER ::   np_ENS = 0   ! ENS scheme
78   INTEGER, PUBLIC, PARAMETER ::   np_ENE = 1   ! ENE scheme
79   INTEGER, PUBLIC, PARAMETER ::   np_ENT = 2   ! ENT scheme (t-point vorticity)
80   INTEGER, PUBLIC, PARAMETER ::   np_EET = 3   ! EET scheme (EEN using e3t)
81   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme
82   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme
83!!an
84   INTEGER, PUBLIC, PARAMETER ::   np_ENS_adKE   = 11   ! ENS scheme - AD scheme (KE only)
85   INTEGER, PUBLIC, PARAMETER ::   np_ENS_adVO   = 12   ! ENS scheme - AD scheme (VOR only)
86   INTEGER, PUBLIC, PARAMETER ::   np_ENS_adKEVO = 13   ! ENS scheme - AD scheme (KE+VOR)
87   INTEGER, PUBLIC, PARAMETER ::   np_ENE_adKE   = 21   ! ENE scheme - AD scheme (KE only)
88   INTEGER, PUBLIC, PARAMETER ::   np_ENE_adVO   = 22   ! ENE scheme - AD scheme (VOR only)
89   INTEGER, PUBLIC, PARAMETER ::   np_ENE_adKEVO = 23   ! ENE scheme - AD scheme (KE+VOR)
90!!an
91
92!!an ds step on pourra spécifier la valeur de ntot = np_COR ou np_COR + np_RVO
93   INTEGER, PUBLIC ::   ncor, nrvm, ntot   ! choice of calculated vorticity
94   !                               ! associated indices:
95   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
96   INTEGER, PUBLIC, PARAMETER ::   np_RVO = 2         ! relative vorticity
97   INTEGER, PUBLIC, PARAMETER ::   np_MET = 3         ! metric term
98   INTEGER, PUBLIC, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity)
99   INTEGER, PUBLIC, PARAMETER ::   np_CME = 5         ! Coriolis + metric term
100
101   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation
102   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -
103   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation
104   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       -
105
106   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4
107   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8
108   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12
109
110   !! * Substitutions
111#  include "do_loop_substitute.h90"
112!!st23
113#  include "domzgr_substitute.h90"
114   !!----------------------------------------------------------------------
115   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
116   !! $Id: dynvor.F90 13154 2020-06-24 13:31:32Z gm $
117   !! Software governed by the CeCILL license (see ./LICENSE)
118   !!----------------------------------------------------------------------
119CONTAINS
120
121   SUBROUTINE dyn_vor( kt, Kbb, Kmm, puu, pvv, Krhs)
122      !!----------------------------------------------------------------------
123      !!
124      !! ** Purpose :   compute the lateral ocean tracer physics.
125      !!
126      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
127      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
128      !!               and planetary vorticity trends) and send them to trd_dyn
129      !!               for futher diagnostics (l_trddyn=T)
130      !!----------------------------------------------------------------------
131      INTEGER ::   ji, jj, jk   ! dummy loop indice
132      INTEGER                             , INTENT( in  ) ::   kt               ! ocean time-step index
133      INTEGER                             , INTENT( in  ) ::   Kmm, Krhs, Kbb   ! ocean time level indices
134      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv         ! ocean velocity field and RHS of momentum equation
135      !
136      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv
137      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zuu, zvv
138      !!----------------------------------------------------------------------
139      !
140      IF( ln_timing )   CALL timing_start('dyn_vor')
141      !
142      IF( l_trddyn ) THEN     !==  trend diagnostics case : split the added trend in two parts  ==!
143         !
144         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
145         !
146         ztrdu(:,:,:) = puu(:,:,:,Krhs)            !* planetary vorticity trend (including Stokes-Coriolis force)
147         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
148         SELECT CASE( nvor_scheme )
149         CASE( np_ENS )           ;   CALL vor_ens( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! enstrophy conserving scheme
150            IF( ln_stcor )            CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
151         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme
152            IF( ln_stcor )            CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
153         CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (T-pts)
154            IF( ln_stcor )            CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
155         CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (een with e3t)
156            IF( ln_stcor )            CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
157         CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy & enstrophy scheme
158            IF( ln_stcor )            CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
159         END SELECT
160         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
161         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
162         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt, Kmm )
163         !
164         IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case)
165            ztrdu(:,:,:) = puu(:,:,:,Krhs)
166            ztrdv(:,:,:) = pvv(:,:,:,Krhs)
167            SELECT CASE( nvor_scheme )
168            CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (T-pts)
169            CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (een with e3t)
170            CASE( np_ENE )           ;   CALL vor_ene( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme
171            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! enstrophy conserving scheme
172            CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy & enstrophy scheme
173            END SELECT
174            ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
175            ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
176            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt, Kmm )
177         ENDIF
178         !
179         DEALLOCATE( ztrdu, ztrdv )
180         !
181      ELSE              !==  total vorticity trend added to the general trend  ==!
182         !
183         SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==!
184         CASE( np_ENT )                        !* energy conserving scheme  (T-pts)
185                             CALL vor_enT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
186            IF( ln_stcor )   CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
187         CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t)
188                             CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
189            IF( ln_stcor )   CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
190         CASE( np_ENE )                        !* energy conserving scheme
191                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )
192                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
193            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
194
195         CASE( np_ENE_adVO )                     !* energy conserving scheme with alternating direction
196            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components
197
198                              !==  Alternative direction - VOR only  ==!
199
200                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
201
202                             ALLOCATE( zuu(jpi,jpj,jpk) )
203                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
204                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
205                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
206                             CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
207                             DEALLOCATE( zuu )
208            ELSE
209                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
210
211                             ALLOCATE( zvv(jpi,jpj,jpk) )
212                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
213                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
214                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
215                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
216                             DEALLOCATE( zvv )
217            ENDIF
218            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
219         CASE( np_ENE_adKE )                     !* energy conserving scheme with alternating direction
220            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components
221
222                                 !==  Alternative direction - KEG only  ==!
223
224                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
225
226                             ALLOCATE( zuu(jpi,jpj,jpk) )
227                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
228                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
229                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
230                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
231                             DEALLOCATE( zuu )
232            ELSE
233
234                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
235
236                             ALLOCATE( zvv(jpi,jpj,jpk) )
237                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
238                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
239                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
240                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm)  , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
241                             DEALLOCATE( zvv )
242            ENDIF
243            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
244
245         CASE( np_ENE_adKEVO )                     !* energy conserving scheme with alternating direction
246         ! equivalent ADE pas sur coriolis
247         ! appel vor_ene(ncor)
248            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components
249
250                              !==  Alternative direction - KE + VOR  ==!
251
252                             ALLOCATE( zuu(jpi,jpj,jpk) )
253                             !puis call que sur nrvm
254                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
255                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   !
256                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
257                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
258                             !puis call que sur nrvm
259                             CALL vor_ene(   kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
260                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )
261                             DEALLOCATE( zuu )
262            ELSE
263
264                             ALLOCATE( zvv(jpi,jpj,jpk) )
265                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
266                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )
267                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
268                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
269                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
270                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )
271                             DEALLOCATE( zvv )
272            ENDIF
273         CASE( np_ENS )                        !* enstrophy conserving scheme
274                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )
275                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend
276            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
277         CASE( np_ENS_adVO )                     !* enstrophy conserving scheme with alternating direction
278            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components
279
280                                 !==  Alternative direction - VOR only  ==!
281
282                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )
283
284                             ALLOCATE( zuu(jpi,jpj,jpk) )
285                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
286                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
287                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
288                             CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
289                             DEALLOCATE( zuu )
290            ELSE
291                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )
292
293                             ALLOCATE( zvv(jpi,jpj,jpk) )
294                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
295                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
296                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
297                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , zvv(:,:,:), pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
298                             DEALLOCATE( zvv )
299            ENDIF
300            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
301         CASE( np_ENS_adKE )                     !* enstrophy conserving scheme with alternating direction
302            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components
303
304                                !==  Alternative direction - KEG only  ==!
305
306                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend
307
308                             ALLOCATE( zuu(jpi,jpj,jpk) )
309                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
310                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
311                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
312                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
313                             DEALLOCATE( zuu )
314            ELSE
315
316                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend
317
318                             ALLOCATE( zvv(jpi,jpj,jpk) )
319                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
320                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
321                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
322                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:)  , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
323                             DEALLOCATE( zvv )
324            ENDIF
325         CASE( np_ENS_adKEVO )                     !* enstrophy conserving scheme with alternating direction
326            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components
327
328                              !==  Alternative direction - KE + VOR  ==!
329
330                             ALLOCATE( zuu(jpi,jpj,jpk) )
331                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
332                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   !
333                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
334                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
335                             CALL vor_ens(   kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
336                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )
337                             DEALLOCATE( zuu )
338            ELSE
339
340                             ALLOCATE( zvv(jpi,jpj,jpk) )
341                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
342                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )
343                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
344                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
345                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
346                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )
347                             DEALLOCATE( zvv )
348            ENDIF
349            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
350         CASE( np_MIX )                        !* mixed ene-ens scheme
351                             CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! relative vorticity or metric trend (ens)
352                             CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! planetary vorticity trend (ene)
353            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
354         CASE( np_EEN )                        !* energy and enstrophy conserving scheme
355                             CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
356            IF( ln_stcor )   CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
357         END SELECT
358         !
359      ENDIF
360      !
361      !                       ! print sum trends (used for debugging)
362      IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor  - Ua: ', mask1=umask,               &
363         &                                tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
364      !
365      IF( ln_timing )   CALL timing_stop('dyn_vor')
366      !
367   END SUBROUTINE dyn_vor
368
369
370   SUBROUTINE vor_enT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
371      !!----------------------------------------------------------------------
372      !!                  ***  ROUTINE vor_enT  ***
373      !!
374      !! ** Purpose :   Compute the now total vorticity trend and add it to
375      !!      the general trend of the momentum equation.
376      !!
377      !! ** Method  :   Trend evaluated using now fields (centered in time)
378      !!       and t-point evaluation of vorticity (planetary and relative).
379      !!       conserves the horizontal kinetic energy.
380      !!         The general trend of momentum is increased due to the vorticity
381      !!       term which is given by:
382      !!          voru = 1/bu  mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t  mj[vn] ]
383      !!          vorv = 1/bv  mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f  mj[un] ]
384      !!       where rvor is the relative vorticity at f-point
385      !!
386      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
387      !!----------------------------------------------------------------------
388      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
389      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
390      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric
391      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities
392      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend
393      !
394      INTEGER  ::   ji, jj, jk           ! dummy loop indices
395      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
396      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx, zwy, zwt   ! 2D workspace
397      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz             ! 3D workspace
398      !!----------------------------------------------------------------------
399      !
400      IF( kt == nit000 ) THEN
401         IF(lwp) WRITE(numout,*)
402         IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
403         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
404      ENDIF
405      !
406      !
407      SELECT CASE( kvor )                 !== relative vorticity considered  ==!
408      CASE ( np_RVO , np_CRV )                  !* relative vorticity at f-point is used
409         DO jk = 1, jpkm1                                ! Horizontal slab
410            DO_2D_10_10
411               zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
412                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
413            END_2D
414            IF( ln_dynvor_msk ) THEN                     ! mask relative vorticity
415               DO_2D_10_10
416                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
417               END_2D
418            ENDIF
419         END DO
420         CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
421         !
422      END SELECT
423
424      !                                                ! ===============
425      DO jk = 1, jpkm1                                 ! Horizontal slab
426         !                                             ! ===============
427         !
428         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
429         !
430         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
431            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm)
432         CASE ( np_RVO )                           !* relative vorticity
433            DO_2D_01_01
434               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   &
435                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) &
436                  &                * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
437            END_2D
438         CASE ( np_MET )                           !* metric term
439            DO_2D_01_01
440               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)     &
441                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) &
442                  &         * e3t(ji,jj,jk,Kmm)
443            END_2D
444         CASE ( np_CRV )                           !* Coriolis + relative vorticity
445            DO_2D_01_01
446               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)      &
447                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) &
448                  &         * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
449            END_2D
450         CASE ( np_CME )                           !* Coriolis + metric
451            DO_2D_01_01
452               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                             &
453                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)    &
454                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) &
455                    &       * e3t(ji,jj,jk,Kmm)
456            END_2D
457         CASE DEFAULT                                             ! error
458            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
459         END SELECT
460         !
461         !                                   !==  compute and add the vorticity term trend  =!
462         DO_2D_00_00
463            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                    &
464               &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   &
465               &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   )
466               !
467            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    &
468               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   &
469               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   )
470         END_2D
471         !                                             ! ===============
472      END DO                                           !   End of slab
473      !                                                ! ===============
474   END SUBROUTINE vor_enT
475
476
477   SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
478      !!----------------------------------------------------------------------
479      !!                  ***  ROUTINE vor_ene  ***
480      !!
481      !! ** Purpose :   Compute the now total vorticity trend and add it to
482      !!      the general trend of the momentum equation.
483      !!
484      !! ** Method  :   Trend evaluated using now fields (centered in time)
485      !!       and the Sadourny (1975) flux form formulation : conserves the
486      !!       horizontal kinetic energy.
487      !!         The general trend of momentum is increased due to the vorticity
488      !!       term which is given by:
489      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v pvv(:,:,:,Kmm)) ]
490      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u puu(:,:,:,Kmm)) ]
491      !!       where rvor is the relative vorticity
492      !!
493      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
494      !!
495      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
496      !!----------------------------------------------------------------------
497      INTEGER                                  , INTENT(in   ) ::   kt          ! ocean time-step index
498      INTEGER                                  , INTENT(in   ) ::   Kmm              ! ocean time level index
499      INTEGER                                  , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
500      REAL(wp), DIMENSION(jpi,jpj,jpk)         , INTENT(inout) ::   pu, pv    ! now velocities
501      REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
502      !
503      INTEGER  ::   ji, jj, jk           ! dummy loop indices
504      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
505      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace
506      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t        ! unpenalised scale factors
507      !!----------------------------------------------------------------------
508      !
509      IF( kt == nit000 ) THEN
510         IF(lwp) WRITE(numout,*)
511         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
512         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
513      ENDIF
514      !
515      !                                                ! ===============
516      DO jk = 1, jpkm1                                 ! Horizontal slab
517         !                                             ! ===============
518         !
519         SELECT CASE( kvor )                 !==  vorticity considered  ==!
520         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
521            zwz(:,:) = ff_f(:,:)
522         CASE ( np_RVO )                           !* relative vorticity
523            DO_2D_10_10
524               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
525                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
526            END_2D
527            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
528               DO_2D_10_10
529                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
530               END_2D
531            ENDIF
532         CASE ( np_MET )                           !* metric term
533            DO_2D_10_10
534               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
535                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
536            END_2D
537         CASE ( np_CRV )                           !* Coriolis + relative vorticity
538            DO_2D_10_10
539               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
540                  &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
541            END_2D
542            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term)
543               DO_2D_10_10
544                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
545               END_2D
546            ENDIF
547         CASE ( np_CME )                           !* Coriolis + metric
548            DO_2D_10_10
549               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
550                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
551            END_2D
552         CASE DEFAULT                                             ! error
553            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
554         END SELECT
555         !
556         IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN     !***  NO alternating direction  ***!
557            !
558            !                                   !==  horizontal fluxes and potential vorticity ==!
559            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
560            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
561            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
562            !
563            !                                   !==  compute and add the vorticity term trend  =!
564            DO_2D_00_00
565               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
566               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
567               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
568               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
569               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 )
570               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 )
571            END_2D
572            !
573            !
574         ELSEIF(       PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : i-component  ***!
575            !
576            !
577            !                                   !==  horizontal fluxes and potential vorticity ==!
578!!anAD remplacer pu par uu(Nnn) et on fait correctement la direction alternée
579!      on transportera (en AD) de la bonne manière
580            !zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
581            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * vv(:,:,jk,Kmm)
582            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
583            !
584            !                                   !==  compute and add the vorticity term trend  =!
585            DO_2D_00_00
586               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
587               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
588               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 )
589            END_2D
590            !
591         ELSEIF( .NOT. PRESENT( pu_rhs ) .AND.       PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : j-component  ***!
592            !
593            !
594            !                                   !==  horizontal fluxes and potential vorticity ==!
595!!anAD remplacer pu par uu(Nnn) et on fait correctement la direction alternée
596            !zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
597            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * uu(:,:,jk,Kmm)
598            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
599            !
600            !                                   !==  compute and add the vorticity term trend  =!
601            DO_2D_00_00
602               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
603               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
604               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 )
605            END_2D
606            !
607         ENDIF
608         !                                             ! ===============
609      END DO                                           !   End of slab
610      !                                                ! ===============
611   END SUBROUTINE vor_ene
612
613
614   SUBROUTINE vor_ens( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
615      !!----------------------------------------------------------------------
616      !!                ***  ROUTINE vor_ens  ***
617      !!
618      !! ** Purpose :   Compute the now total vorticity trend and add it to
619      !!      the general trend of the momentum equation.
620      !!
621      !! ** Method  :   Trend evaluated using now fields (centered in time)
622      !!      and the Sadourny (1975) flux FORM formulation : conserves the
623      !!      potential enstrophy of a horizontally non-divergent flow. the
624      !!      trend of the vorticity term is given by:
625      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v pvv(:,:,:,Kmm)) ]
626      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u puu(:,:,:,Kmm)) ]
627      !!      Add this trend to the general momentum trend:
628      !!          (u(rhs),v(Krhs)) = (u(rhs),v(Krhs)) + ( voru , vorv )
629      !!
630      !! ** Action : - Update (pu_rhs,pv_rhs)) arrays with the now vorticity term trend
631      !!
632      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
633      !!----------------------------------------------------------------------
634      INTEGER                                  , INTENT(in   ) ::   kt          ! ocean time-step index
635      INTEGER                                  , INTENT(in   ) ::   Kmm              ! ocean time level index
636      INTEGER                                  , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
637      REAL(wp), DIMENSION(jpi,jpj,jpk)         , INTENT(inout) ::   pu, pv    ! now velocities
638      REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
639      !
640      INTEGER  ::   ji, jj, jk   ! dummy loop indices
641      REAL(wp) ::   zuav, zvau   ! local scalars
642      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! 2D workspace
643      REAL(wp), DIMENSION(jpi,jpj) ::   zr1_e1e2f            ! 2D workspace penalised
644      !!----------------------------------------------------------------------
645      !
646      IF( kt == nit000 ) THEN
647         IF(lwp) WRITE(numout,*)
648         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
649         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
650      ENDIF
651      !
652      zr1_e1e2f(:,:) = r1_e1e2f(:,:)
653# if defined key_bvp
654      zr1_e1e2f(:,:) = r1_rpof(:,:,1) * zr1_e1e2f(:,:)
655# endif
656      !                                                ! ===============
657      DO jk = 1, jpkm1                                 ! Horizontal slab
658         !                                             ! ===============
659         !
660         SELECT CASE( kvor )                 !==  vorticity considered  ==!
661         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
662            zwz(:,:) = ff_f(:,:)
663         CASE ( np_RVO )                           !* relative vorticity
664            DO_2D_10_10
665               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
666                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * zr1_e1e2f(ji,jj)
667            END_2D
668            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
669               DO_2D_10_10
670                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
671               END_2D
672            ENDIF
673         CASE ( np_MET )                           !* metric term
674            DO_2D_10_10
675               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
676                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
677            END_2D
678         CASE ( np_CRV )                           !* Coriolis + relative vorticity
679            DO_2D_10_10
680               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
681                  &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * zr1_e1e2f(ji,jj)
682            END_2D
683            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term)
684               DO_2D_10_10
685                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
686               END_2D
687            ENDIF
688         CASE ( np_CME )                           !* Coriolis + metric
689            DO_2D_10_10
690               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
691                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
692            END_2D
693         CASE DEFAULT                                             ! error
694            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
695         END SELECT
696         !
697         !
698!!an wut ? v et u
699         IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN     !***  NO alternating direction  ***!
700            !
701            !                                   !==  horizontal fluxes and potential vorticity ==!
702            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
703            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
704            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
705            !
706            !                                   !==  compute and add the vorticity term trend  =!
707            DO_2D_00_00
708               zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
709                  &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
710               zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
711                  &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
712               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
713               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
714            END_2D
715            !
716         ELSEIF(       PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : i-component  ***!
717            !
718            !
719            !                                   !==  horizontal fluxes and potential vorticity ==!
720            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
721            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
722            !
723            !                                   !==  compute and add the vorticity term trend  =!
724            DO_2D_00_00
725               zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
726                  &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
727               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
728            END_2D
729            !
730         ELSEIF( .NOT. PRESENT( pu_rhs ) .AND.       PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : j-component  ***!
731            !
732            !
733            !                                   !==  horizontal fluxes and potential vorticity ==!
734            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
735            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
736            !
737            !                                   !==  compute and add the vorticity term trend  =!
738            DO_2D_00_00
739               zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
740                  &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
741               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
742            END_2D
743            !
744         ENDIF
745         !                                             ! ===============
746      END DO                                           !   End of slab
747      !                                                ! ===============
748   END SUBROUTINE vor_ens
749
750
751   SUBROUTINE vor_een( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
752      !!----------------------------------------------------------------------
753      !!                ***  ROUTINE vor_een  ***
754      !!
755      !! ** Purpose :   Compute the now total vorticity trend and add it to
756      !!      the general trend of the momentum equation.
757      !!
758      !! ** Method  :   Trend evaluated using now fields (centered in time)
759      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
760      !!      both the horizontal kinetic energy and the potential enstrophy
761      !!      when horizontal divergence is zero (see the NEMO documentation)
762      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
763      !!
764      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
765      !!
766      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
767      !!----------------------------------------------------------------------
768      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
769      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
770      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
771      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
772      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
773      !
774      INTEGER  ::   ji, jj, jk   ! dummy loop indices
775      INTEGER  ::   ierr         ! local integer
776      REAL(wp) ::   zua, zva     ! local scalars
777      REAL(wp) ::   zmsk, ze3f   ! local scalars
778      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy , z1_e3f
779      REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse
780      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz
781      !!----------------------------------------------------------------------
782      !
783      IF( kt == nit000 ) THEN
784         IF(lwp) WRITE(numout,*)
785         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
786         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
787      ENDIF
788      !
789      !                                                ! ===============
790      DO jk = 1, jpkm1                                 ! Horizontal slab
791         !                                             ! ===============
792         !
793         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
794         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
795            DO_2D_10_10
796               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
797                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
798                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
799                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
800               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
801               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
802               ENDIF
803            END_2D
804         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
805            DO_2D_10_10
806               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
807                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
808                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
809                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
810               zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
811                  &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
812               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3f
813               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
814               ENDIF
815            END_2D
816         END SELECT
817         !
818         SELECT CASE( kvor )                 !==  vorticity considered  ==!
819         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
820            DO_2D_10_10
821               zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
822            END_2D
823         CASE ( np_RVO )                           !* relative vorticity
824            DO_2D_10_10
825               zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
826                  &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
827            END_2D
828         CASE ( np_MET )                           !* metric term
829            DO_2D_10_10
830               zwz(ji,jj,jk) = (   ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
831                  &              - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
832            END_2D
833         CASE ( np_CRV )                           !* Coriolis + relative vorticity
834            DO_2D_10_10
835               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
836                  &                              - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  )   &
837                  &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
838            END_2D
839         CASE ( np_CME )                           !* Coriolis + metric
840            DO_2D_10_10
841               zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
842                  &                            - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
843            END_2D
844         CASE DEFAULT                                             ! error
845            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
846         END SELECT
847         !
848         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
849            DO_2D_10_10
850               zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
851            END_2D
852         ENDIF
853      END DO                                           !   End of slab
854         !
855      CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
856
857      DO jk = 1, jpkm1                                 ! Horizontal slab
858         !
859         !                                   !==  horizontal fluxes  ==!
860         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
861         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
862
863         !                                   !==  compute and add the vorticity term trend  =!
864         jj = 2
865         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
866         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
867               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
868               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
869               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
870               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
871         END DO
872         DO jj = 3, jpj
873            DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3
874               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
875               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
876               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
877               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
878            END DO
879         END DO
880         DO_2D_00_00
881            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
882               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
883            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
884               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
885            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
886            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
887         END_2D
888         !                                             ! ===============
889      END DO                                           !   End of slab
890      !                                                ! ===============
891   END SUBROUTINE vor_een
892
893
894
895   SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
896      !!----------------------------------------------------------------------
897      !!                ***  ROUTINE vor_eeT  ***
898      !!
899      !! ** Purpose :   Compute the now total vorticity trend and add it to
900      !!      the general trend of the momentum equation.
901      !!
902      !! ** Method  :   Trend evaluated using now fields (centered in time)
903      !!      and the Arakawa and Lamb (1980) vector form formulation using
904      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
905      !!      The change consists in
906      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
907      !!
908      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
909      !!
910      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
911      !!----------------------------------------------------------------------
912      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
913      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
914      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
915      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
916      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
917      !
918      INTEGER  ::   ji, jj, jk     ! dummy loop indices
919      INTEGER  ::   ierr           ! local integer
920      REAL(wp) ::   zua, zva       ! local scalars
921      REAL(wp) ::   zmsk, z1_e3t   ! local scalars
922      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy
923      REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse
924      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz
925      !!----------------------------------------------------------------------
926      !
927      IF( kt == nit000 ) THEN
928         IF(lwp) WRITE(numout,*)
929         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
930         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
931      ENDIF
932      !
933      !                                                ! ===============
934      DO jk = 1, jpkm1                                 ! Horizontal slab
935         !                                             ! ===============
936         !
937         !
938         SELECT CASE( kvor )                 !==  vorticity considered  ==!
939         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
940            DO_2D_10_10
941               zwz(ji,jj,jk) = ff_f(ji,jj)
942            END_2D
943         CASE ( np_RVO )                           !* relative vorticity
944            DO_2D_10_10
945               zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
946                  &             - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
947                  &          * r1_e1e2f(ji,jj)
948            END_2D
949         CASE ( np_MET )                           !* metric term
950            DO_2D_10_10
951               zwz(ji,jj,jk) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
952                  &          - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
953            END_2D
954         CASE ( np_CRV )                           !* Coriolis + relative vorticity
955            DO_2D_10_10
956               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
957                  &                              - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
958                  &                         * r1_e1e2f(ji,jj)    )
959            END_2D
960         CASE ( np_CME )                           !* Coriolis + metric
961            DO_2D_10_10
962               zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
963                  &                        - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
964            END_2D
965         CASE DEFAULT                                             ! error
966            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
967         END SELECT
968         !
969         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
970            DO_2D_10_10
971               zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
972            END_2D
973         ENDIF
974      END DO
975      !
976      CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
977      !
978      DO jk = 1, jpkm1                                 ! Horizontal slab
979         !
980         !                                   !==  horizontal fluxes  ==!
981         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
982         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
983
984         !                                   !==  compute and add the vorticity term trend  =!
985         jj = 2
986         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
987         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
988               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
989               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
990               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
991               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
992               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
993         END DO
994         DO jj = 3, jpj
995            DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3
996               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
997               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
998               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
999               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
1000               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
1001            END DO
1002         END DO
1003         DO_2D_00_00
1004            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
1005               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
1006            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
1007               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
1008            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
1009            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
1010         END_2D
1011         !                                             ! ===============
1012      END DO                                           !   End of slab
1013      !                                                ! ===============
1014   END SUBROUTINE vor_eeT
1015
1016
1017   SUBROUTINE dyn_vor_init
1018      !!---------------------------------------------------------------------
1019      !!                  ***  ROUTINE dyn_vor_init  ***
1020      !!
1021      !! ** Purpose :   Control the consistency between cpp options for
1022      !!              tracer advection schemes
1023      !!----------------------------------------------------------------------
1024      INTEGER ::   ji, jj, jk    ! dummy loop indices
1025      INTEGER ::   ioptio, ios   ! local integer
1026      !!
1027      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   &
1028         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk,   &
1029         &                 ln_dynvor_ens_adVO, ln_dynvor_ens_adKE, ln_dynvor_ens_adKEVO, &   ! Alternative direction parameters
1030         &                 ln_dynvor_ene_adVO, ln_dynvor_ene_adKE, ln_dynvor_ene_adKEVO
1031      !!----------------------------------------------------------------------
1032      !
1033      IF(lwp) THEN
1034         WRITE(numout,*)
1035         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
1036         WRITE(numout,*) '~~~~~~~~~~~~'
1037      ENDIF
1038      !
1039      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
1040901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' )
1041      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
1042902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' )
1043      IF(lwm) WRITE ( numond, namdyn_vor )
1044      !
1045      IF(lwp) THEN                    ! Namelist print
1046         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
1047         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
1048         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene
1049         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT
1050         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT
1051         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
1052         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
1053         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
1054         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
1055      ENDIF
1056
1057!!an      IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available')
1058
1059!!gm  this should be removed when choosing a unique strategy for fmask at the coast
1060      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
1061      ! at angles with three ocean points and one land point
1062      IF(lwp) WRITE(numout,*)
1063      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
1064      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
1065         DO_3D_10_10( 1, jpk )
1066            IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
1067               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
1068         END_3D
1069         !
1070         CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
1071         !
1072      ENDIF
1073!!gm end
1074
1075      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
1076      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF
1077      IF( ln_dynvor_ens_adVO ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENS_adVO   ;   ENDIF
1078      IF( ln_dynvor_ens_adKE ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENS_adKE   ;   ENDIF
1079      IF( ln_dynvor_ens_adKEVO ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS_adKEVO   ;   ENDIF
1080      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
1081      IF( ln_dynvor_ene_adVO ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENE_adVO   ;   ENDIF
1082      IF( ln_dynvor_ene_adKE ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENE_adKE   ;   ENDIF
1083      IF( ln_dynvor_ene_adKEVO ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE_adKEVO   ;   ENDIF
1084      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF
1085      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF
1086      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF
1087      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF
1088      !
1089      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
1090      !
1091      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
1092      ncor = np_COR                       ! planetary vorticity
1093      SELECT CASE( n_dynadv )
1094      CASE( np_LIN_dyn )
1095         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
1096         nrvm = np_COR        ! planetary vorticity
1097         ntot = np_COR        !     -         -
1098      CASE( np_VEC_c2  )
1099         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity'
1100         nrvm = np_RVO        ! relative vorticity
1101         ntot = np_CRV        ! relative + planetary vorticity
1102      CASE( np_FLX_c2 , np_FLX_ubs  )
1103         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
1104         nrvm = np_MET        ! metric term
1105         ntot = np_CME        ! Coriolis + metric term
1106         !
1107         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term:
1108         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
1109            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
1110            DO_2D_00_00
1111               di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
1112               dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
1113            END_2D
1114            CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1. , dj_e1v_2, 'T', -1. )   ! Lateral boundary conditions
1115            !
1116         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
1117            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
1118            DO_2D_10_10
1119               di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
1120               dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
1121            END_2D
1122            CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1. , dj_e1u_2e1e2f, 'F', -1. )   ! Lateral boundary conditions
1123         END SELECT
1124         !
1125      END SELECT
1126
1127      IF(lwp) THEN                   ! Print the choice
1128         WRITE(numout,*)
1129         SELECT CASE( nvor_scheme )
1130
1131         CASE( np_ENS    )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
1132         CASE( np_ENS_adVO ) ;   WRITE(numout,*) '   ==>>>   AD enstrophy conserving scheme (ENS_adVO) on vorticity only'
1133         CASE( np_ENS_adKE ) ;   WRITE(numout,*) '   ==>>>   AD enstrophy conserving scheme (ENS_adKE) on kinetic energy only'
1134         CASE( np_ENS_adKEVO ) ;   WRITE(numout,*) '   ==>>>   AD enstrophy conserving scheme (ENS_adKEVO) on kinetic energy and vorticity'
1135
1136         CASE( np_ENE    )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
1137         CASE( np_ENE_adVO ) ;   WRITE(numout,*) '   ==>>>   AD energy conserving scheme (Coriolis at F-points) (ENE_adVO) on vorticity only'
1138         CASE( np_ENE_adKE ) ;   WRITE(numout,*) '   ==>>>   AD energy conserving scheme (Coriolis at F-points) (ENE_adKE) on kinetic energy only'
1139         CASE( np_ENE_adKEVO ) ;   WRITE(numout,*) '   ==>>>   AD energy conserving scheme (Coriolis at F-points) (ENE_adKEVO) on kinetic energy and vorticity'
1140
1141         CASE( np_ENT    )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
1142         CASE( np_EET    )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
1143         CASE( np_EEN    )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
1144         CASE( np_MIX    )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
1145         END SELECT
1146      ENDIF
1147      !
1148   END SUBROUTINE dyn_vor_init
1149
1150   !!==============================================================================
1151END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.