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

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynvor.F90 @ 14618

Last change on this file since 14618 was 14547, checked in by techene, 3 years ago

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