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/UKMO/dev_r10037_dynvor_EEUV/src/OCE/DYN – NEMO

source: NEMO/branches/UKMO/dev_r10037_dynvor_EEUV/src/OCE/DYN/dynvor.F90 @ 10094

Last change on this file since 10094 was 10094, checked in by davestorkey, 5 years ago

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