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/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN – NEMO

source: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynvor.F90 @ 14787

Last change on this file since 14787 was 14787, checked in by hadcv, 3 years ago

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