Context Navigation

dynvor.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

Property svn:keywords set to Id

File size: 52.3 KB

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

Note: See TracBrowser for help on using the repository browser.

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynvor.F90 @ 12340

Download in other formats: