Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynvor.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 4 years ago
The Dr Hook changes from my perl code.
File size: 43.4 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	!!----------------------------------------------------------------------
19
20	!!----------------------------------------------------------------------
21	!! dyn_vor : Update the momentum trend with the vorticity trend
22	!! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T)
23	!! vor_ene : energy conserving scheme (ln_dynvor_ene=T)
24	!! vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=T)
25	!! vor_een : energy and enstrophy conserving (ln_dynvor_een=T)
26	!! dyn_vor_init : set and control of the different vorticity option
27	!!----------------------------------------------------------------------
28	USE oce ! ocean dynamics and tracers
29	USE dom_oce ! ocean space and time domain
30	USE dommsk ! ocean mask
31	USE dynadv ! momentum advection (use ln_dynadv_vec value)
32	USE trd_oce ! trends: ocean variables
33	USE trddyn ! trend manager: dynamics
34	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
35	USE prtctl ! Print control
36	USE in_out_manager ! I/O manager
37	USE lib_mpp ! MPP library
38	USE wrk_nemo ! Memory Allocation
39	USE timing ! Timing
40	USE lib_fortran
41
42
43	USE yomhook, ONLY: lhook, dr_hook
44	USE parkind1, ONLY: jprb, jpim
45
46	IMPLICIT NONE
47	PRIVATE
48
49	PUBLIC dyn_vor ! routine called by step.F90
50	PUBLIC dyn_vor_init ! routine called by opa.F90
51
52	! !!* Namelist namdyn_vor: vorticity term
53	LOGICAL, PUBLIC :: ln_dynvor_ene !: energy conserving scheme
54	LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme
55	LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme
56	LOGICAL, PUBLIC :: ln_dynvor_een !: energy and enstrophy conserving scheme
57	LOGICAL, PUBLIC :: ln_dynvor_een_old !: energy and enstrophy conserving scheme (original formulation)
58
59	INTEGER :: nvor = 0 ! type of vorticity trend used
60	INTEGER :: ncor = 1 ! coriolis
61	INTEGER :: nrvm = 2 ! =2 relative vorticity ; =3 metric term
62	INTEGER :: ntot = 4 ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term
63
64	!! * Substitutions
65	# include "domzgr_substitute.h90"
66	# include "vectopt_loop_substitute.h90"
67	!!----------------------------------------------------------------------
68	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
69	!! $Id$
70	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
71	!!----------------------------------------------------------------------
72	CONTAINS
73
74	SUBROUTINE dyn_vor( kt )
75	!!----------------------------------------------------------------------
76	!!
77	!! ** Purpose : compute the lateral ocean tracer physics.
78	!!
79	!! ** Action : - Update (ua,va) with the now vorticity term trend
80	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
81	!! and planetary vorticity trends) and send them to trd_dyn
82	!! for futher diagnostics (l_trddyn=T)
83	!!----------------------------------------------------------------------
84	INTEGER, INTENT( in ) :: kt ! ocean time-step index
85	!
86	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv
87	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
88	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
89	REAL(KIND=jprb) :: zhook_handle
90
91	CHARACTER(LEN=*), PARAMETER :: RoutineName='DYN_VOR'
92
93	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
94
95	!!----------------------------------------------------------------------
96	!
97	IF( nn_timing == 1 ) CALL timing_start('dyn_vor')
98	!
99	IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
100	!
101	! ! vorticity term
102	SELECT CASE ( nvor ) ! compute the vorticity trend and add it to the general trend
103	!
104	CASE ( -1 ) ! esopa: test all possibility with control print
105	CALL vor_ene( kt, ntot, ua, va )
106	CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, &
107	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
108	CALL vor_ens( kt, ntot, ua, va )
109	CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, &
110	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
111	CALL vor_mix( kt )
112	CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, &
113	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
114	CALL vor_een( kt, ntot, ua, va )
115	CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, &
116	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
117	!
118	CASE ( 0 ) ! energy conserving scheme
119	IF( l_trddyn ) THEN
120	ztrdu(:,:,:) = ua(:,:,:)
121	ztrdv(:,:,:) = va(:,:,:)
122	CALL vor_ene( kt, nrvm, ua, va ) ! relative vorticity or metric trend
123	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
124	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
125	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
126	ztrdu(:,:,:) = ua(:,:,:)
127	ztrdv(:,:,:) = va(:,:,:)
128	CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend
129	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
130	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
131	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
132	ELSE
133	CALL vor_ene( kt, ntot, ua, va ) ! total vorticity
134	ENDIF
135	!
136	CASE ( 1 ) ! enstrophy conserving scheme
137	IF( l_trddyn ) THEN
138	ztrdu(:,:,:) = ua(:,:,:)
139	ztrdv(:,:,:) = va(:,:,:)
140	CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend
141	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
142	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
143	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
144	ztrdu(:,:,:) = ua(:,:,:)
145	ztrdv(:,:,:) = va(:,:,:)
146	CALL vor_ens( kt, ncor, ua, va ) ! planetary vorticity trend
147	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
148	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
149	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
150	ELSE
151	CALL vor_ens( kt, ntot, ua, va ) ! total vorticity
152	ENDIF
153	!
154	CASE ( 2 ) ! mixed ene-ens scheme
155	IF( l_trddyn ) THEN
156	ztrdu(:,:,:) = ua(:,:,:)
157	ztrdv(:,:,:) = va(:,:,:)
158	CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens)
159	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
160	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
161	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
162	ztrdu(:,:,:) = ua(:,:,:)
163	ztrdv(:,:,:) = va(:,:,:)
164	CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene)
165	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
166	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
167	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
168	ELSE
169	CALL vor_mix( kt ) ! total vorticity (mix=ens-ene)
170	ENDIF
171	!
172	CASE ( 3 ) ! energy and enstrophy conserving scheme
173	IF( l_trddyn ) THEN
174	ztrdu(:,:,:) = ua(:,:,:)
175	ztrdv(:,:,:) = va(:,:,:)
176	CALL vor_een( kt, nrvm, ua, va ) ! relative vorticity or metric trend
177	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
178	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
179	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
180	ztrdu(:,:,:) = ua(:,:,:)
181	ztrdv(:,:,:) = va(:,:,:)
182	CALL vor_een( kt, ncor, ua, va ) ! planetary vorticity trend
183	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
184	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
185	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
186	ELSE
187	CALL vor_een( kt, ntot, ua, va ) ! total vorticity
188	ENDIF
189	!
190	END SELECT
191	!
192	! ! print sum trends (used for debugging)
193	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, &
194	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
195	!
196	IF( l_trddyn ) CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
197	!
198	IF( nn_timing == 1 ) CALL timing_stop('dyn_vor')
199	!
200	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
201	END SUBROUTINE dyn_vor
202
203
204	SUBROUTINE vor_ene( kt, kvor, pua, pva )
205	!!----------------------------------------------------------------------
206	!! * ROUTINE vor_ene *
207	!!
208	!! ** Purpose : Compute the now total vorticity trend and add it to
209	!! the general trend of the momentum equation.
210	!!
211	!! ** Method : Trend evaluated using now fields (centered in time)
212	!! and the Sadourny (1975) flux form formulation : conserves the
213	!! horizontal kinetic energy.
214	!! The trend of the vorticity term is given by:
215	!! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
216	!! voru = 1/e1u mj-1[ (rotn+f)/e3f mi(e1v*e3v vn) ]
217	!! vorv = 1/e2v mi-1[ (rotn+f)/e3f mj(e2u*e3u un) ]
218	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
219	!! voru = 1/e1u mj-1[ (rotn+f) mi(e1v vn) ]
220	!! vorv = 1/e2v mi-1[ (rotn+f) mj(e2u un) ]
221	!! Add this trend to the general momentum trend (ua,va):
222	!! (ua,va) = (ua,va) + ( voru , vorv )
223	!!
224	!! ** Action : - Update (ua,va) with the now vorticity term trend
225	!!
226	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
227	!!----------------------------------------------------------------------
228	!
229	INTEGER , INTENT(in ) :: kt ! ocean time-step index
230	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
231	! ! =nrvm (relative vorticity or metric)
232	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
233	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
234	!
235	INTEGER :: ji, jj, jk ! dummy loop indices
236	REAL(wp) :: zx1, zy1, zfact2, zx2, zy2 ! local scalars
237	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz
238	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
239	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
240	REAL(KIND=jprb) :: zhook_handle
241
242	CHARACTER(LEN=*), PARAMETER :: RoutineName='VOR_ENE'
243
244	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
245
246	!!----------------------------------------------------------------------
247	!
248	IF( nn_timing == 1 ) CALL timing_start('vor_ene')
249	!
250	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )
251	!
252	IF( kt == nit000 ) THEN
253	IF(lwp) WRITE(numout,*)
254	IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
255	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
256	ENDIF
257
258	zfact2 = 0.5 * 0.5 ! Local constant initialization
259
260	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
261	! ! ===============
262	DO jk = 1, jpkm1 ! Horizontal slab
263	! ! ===============
264	!
265	! Potential vorticity and horizontal fluxes
266	! -----------------------------------------
267	SELECT CASE( kvor ) ! vorticity considered
268	CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis)
269	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity
270	CASE ( 3 ) ! metric term
271	DO jj = 1, jpjm1
272	DO ji = 1, fs_jpim1 ! vector opt.
273	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
274	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
275	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
276	END DO
277	END DO
278	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity)
279	CASE ( 5 ) ! total (coriolis + metric)
280	DO jj = 1, jpjm1
281	DO ji = 1, fs_jpim1 ! vector opt.
282	zwz(ji,jj) = ( ff (ji,jj) &
283	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
284	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
285	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
286	& )
287	END DO
288	END DO
289	END SELECT
290
291	IF( ln_sco ) THEN
292	zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
293	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
294	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
295	ELSE
296	zwx(:,:) = e2u(:,:) * un(:,:,jk)
297	zwy(:,:) = e1v(:,:) * vn(:,:,jk)
298	ENDIF
299
300	! Compute and add the vorticity term trend
301	! ----------------------------------------
302	DO jj = 2, jpjm1
303	DO ji = fs_2, fs_jpim1 ! vector opt.
304	zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
305	zy2 = zwy(ji,jj ) + zwy(ji+1,jj )
306	zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
307	zx2 = zwx(ji ,jj) + zwx(ji ,jj+1)
308	pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
309	pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
310	END DO
311	END DO
312	! ! ===============
313	END DO ! End of slab
314	! ! ===============
315	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )
316	!
317	IF( nn_timing == 1 ) CALL timing_stop('vor_ene')
318	!
319	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
320	END SUBROUTINE vor_ene
321
322
323	SUBROUTINE vor_mix( kt )
324	!!----------------------------------------------------------------------
325	!! * ROUTINE vor_mix *
326	!!
327	!! ** Purpose : Compute the now total vorticity trend and add it to
328	!! the general trend of the momentum equation.
329	!!
330	!! ** Method : Trend evaluated using now fields (centered in time)
331	!! Mixte formulation : conserves the potential enstrophy of a hori-
332	!! zontally non-divergent flow for (rotzu x uh), the relative vor-
333	!! ticity term and the horizontal kinetic energy for (f x uh), the
334	!! coriolis term. the now trend of the vorticity term is given by:
335	!! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
336	!! voru = 1/e1u mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]
337	!! +1/e1u mj-1[ f/e3f mi(e1v*e3v vn) ]
338	!! vorv = 1/e2v mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]
339	!! +1/e2v mi-1[ f/e3f mj(e2u*e3u un) ]
340	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
341	!! voru = 1/e1u mj-1(rotn) mj-1[ mi(e1v vn) ]
342	!! +1/e1u mj-1[ f mi(e1v vn) ]
343	!! vorv = 1/e2v mi-1(rotn) mi-1[ mj(e2u un) ]
344	!! +1/e2v mi-1[ f mj(e2u un) ]
345	!! Add this now trend to the general momentum trend (ua,va):
346	!! (ua,va) = (ua,va) + ( voru , vorv )
347	!!
348	!! ** Action : - Update (ua,va) arrays with the now vorticity term trend
349	!!
350	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
351	!!----------------------------------------------------------------------
352	!
353	INTEGER, INTENT(in) :: kt ! ocean timestep index
354	!
355	INTEGER :: ji, jj, jk ! dummy loop indices
356	REAL(wp) :: zfact1, zua, zcua, zx1, zy1 ! local scalars
357	REAL(wp) :: zfact2, zva, zcva, zx2, zy2 ! - -
358	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
359	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
360	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
361	REAL(KIND=jprb) :: zhook_handle
362
363	CHARACTER(LEN=*), PARAMETER :: RoutineName='VOR_MIX'
364
365	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
366
367	!!----------------------------------------------------------------------
368	!
369	IF( nn_timing == 1 ) CALL timing_start('vor_mix')
370	!
371	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww )
372	!
373	IF( kt == nit000 ) THEN
374	IF(lwp) WRITE(numout,*)
375	IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
376	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
377	ENDIF
378
379	zfact1 = 0.5 * 0.25 ! Local constant initialization
380	zfact2 = 0.5 * 0.5
381
382	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
383	! ! ===============
384	DO jk = 1, jpkm1 ! Horizontal slab
385	! ! ===============
386	!
387	! Relative and planetary potential vorticity and horizontal fluxes
388	! ----------------------------------------------------------------
389	IF( ln_sco ) THEN
390	IF( ln_dynadv_vec ) THEN
391	zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
392	ELSE
393	DO jj = 1, jpjm1
394	DO ji = 1, fs_jpim1 ! vector opt.
395	zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
396	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
397	& * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
398	END DO
399	END DO
400	ENDIF
401	zwz(:,:) = ff (:,:) / fse3f(:,:,jk)
402	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
403	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
404	ELSE
405	IF( ln_dynadv_vec ) THEN
406	zww(:,:) = rotn(:,:,jk)
407	ELSE
408	DO jj = 1, jpjm1
409	DO ji = 1, fs_jpim1 ! vector opt.
410	zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
411	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
412	& * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
413	END DO
414	END DO
415	ENDIF
416	zwz(:,:) = ff (:,:)
417	zwx(:,:) = e2u(:,:) * un(:,:,jk)
418	zwy(:,:) = e1v(:,:) * vn(:,:,jk)
419	ENDIF
420
421	! Compute and add the vorticity term trend
422	! ----------------------------------------
423	DO jj = 2, jpjm1
424	DO ji = fs_2, fs_jpim1 ! vector opt.
425	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
426	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
427	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
428	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
429	! enstrophy conserving formulation for relative vorticity term
430	zua = zfact1 * ( zww(ji ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
431	zva =-zfact1 * ( zww(ji-1,jj ) + zww(ji,jj) ) * ( zx1 + zx2 )
432	! energy conserving formulation for planetary vorticity term
433	zcua = zfact2 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
434	zcva =-zfact2 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
435	! mixed vorticity trend added to the momentum trends
436	ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
437	va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
438	END DO
439	END DO
440	! ! ===============
441	END DO ! End of slab
442	! ! ===============
443	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww )
444	!
445	IF( nn_timing == 1 ) CALL timing_stop('vor_mix')
446	!
447	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
448	END SUBROUTINE vor_mix
449
450
451	SUBROUTINE vor_ens( kt, kvor, pua, pva )
452	!!----------------------------------------------------------------------
453	!! * ROUTINE vor_ens *
454	!!
455	!! ** Purpose : Compute the now total vorticity trend and add it to
456	!! the general trend of the momentum equation.
457	!!
458	!! ** Method : Trend evaluated using now fields (centered in time)
459	!! and the Sadourny (1975) flux FORM formulation : conserves the
460	!! potential enstrophy of a horizontally non-divergent flow. the
461	!! trend of the vorticity term is given by:
462	!! * s-coordinate (ln_sco=T), the e3. are inside the derivative:
463	!! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ]
464	!! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ]
465	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
466	!! voru = 1/e1u mj-1[ rotn+f ] mj-1[ mi(e1v vn) ]
467	!! vorv = 1/e2v mi-1[ rotn+f ] mi-1[ mj(e2u un) ]
468	!! Add this trend to the general momentum trend (ua,va):
469	!! (ua,va) = (ua,va) + ( voru , vorv )
470	!!
471	!! ** Action : - Update (ua,va) arrays with the now vorticity term trend
472	!!
473	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
474	!!----------------------------------------------------------------------
475	!
476	INTEGER , INTENT(in ) :: kt ! ocean time-step index
477	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
478	! ! =nrvm (relative vorticity or metric)
479	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
480	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
481	!
482	INTEGER :: ji, jj, jk ! dummy loop indices
483	REAL(wp) :: zfact1, zuav, zvau ! temporary scalars
484	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
485	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
486	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
487	REAL(KIND=jprb) :: zhook_handle
488
489	CHARACTER(LEN=*), PARAMETER :: RoutineName='VOR_ENS'
490
491	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
492
493	!!----------------------------------------------------------------------
494	!
495	IF( nn_timing == 1 ) CALL timing_start('vor_ens')
496	!
497	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )
498	!
499	IF( kt == nit000 ) THEN
500	IF(lwp) WRITE(numout,*)
501	IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
502	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
503	ENDIF
504
505	zfact1 = 0.5 * 0.25 ! Local constant initialization
506
507	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
508	! ! ===============
509	DO jk = 1, jpkm1 ! Horizontal slab
510	! ! ===============
511	!
512	! Potential vorticity and horizontal fluxes
513	! -----------------------------------------
514	SELECT CASE( kvor ) ! vorticity considered
515	CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis)
516	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity
517	CASE ( 3 ) ! metric term
518	DO jj = 1, jpjm1
519	DO ji = 1, fs_jpim1 ! vector opt.
520	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
521	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
522	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
523	END DO
524	END DO
525	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity)
526	CASE ( 5 ) ! total (coriolis + metric)
527	DO jj = 1, jpjm1
528	DO ji = 1, fs_jpim1 ! vector opt.
529	zwz(ji,jj) = ( ff (ji,jj) &
530	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
531	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
532	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
533	& )
534	END DO
535	END DO
536	END SELECT
537	!
538	IF( ln_sco ) THEN
539	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
540	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
541	zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
542	zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
543	zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
544	END DO
545	END DO
546	ELSE
547	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
548	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
549	zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
550	zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
551	END DO
552	END DO
553	ENDIF
554	!
555	! Compute and add the vorticity term trend
556	! ----------------------------------------
557	DO jj = 2, jpjm1
558	DO ji = fs_2, fs_jpim1 ! vector opt.
559	zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
560	& + zwy(ji ,jj ) + zwy(ji+1,jj ) )
561	zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
562	& + zwx(ji ,jj ) + zwx(ji ,jj+1) )
563	pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) )
564	pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) )
565	END DO
566	END DO
567	! ! ===============
568	END DO ! End of slab
569	! ! ===============
570	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )
571	!
572	IF( nn_timing == 1 ) CALL timing_stop('vor_ens')
573	!
574	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
575	END SUBROUTINE vor_ens
576
577
578	SUBROUTINE vor_een( kt, kvor, pua, pva )
579	!!----------------------------------------------------------------------
580	!! * ROUTINE vor_een *
581	!!
582	!! ** Purpose : Compute the now total vorticity trend and add it to
583	!! the general trend of the momentum equation.
584	!!
585	!! ** Method : Trend evaluated using now fields (centered in time)
586	!! and the Arakawa and Lamb (1980) flux form formulation : conserves
587	!! both the horizontal kinetic energy and the potential enstrophy
588	!! when horizontal divergence is zero (see the NEMO documentation)
589	!! Add this trend to the general momentum trend (ua,va).
590	!!
591	!! ** Action : - Update (ua,va) with the now vorticity term trend
592	!!
593	!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
594	!!----------------------------------------------------------------------
595	!
596	INTEGER , INTENT(in ) :: kt ! ocean time-step index
597	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
598	! ! =nrvm (relative vorticity or metric)
599	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
600	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
601	!!
602	INTEGER :: ji, jj, jk ! dummy loop indices
603	INTEGER :: ierr ! local integer
604	REAL(wp) :: zfac12, zua, zva ! local scalars
605	REAL(wp) :: zmsk, ze3 ! local scalars
606	! ! 3D workspace
607	REAL(wp), POINTER , DIMENSION(:,: ) :: zwx, zwy, zwz
608	REAL(wp), POINTER , DIMENSION(:,: ) :: ztnw, ztne, ztsw, ztse
609	#if defined key_vvl
610	REAL(wp), POINTER , DIMENSION(:,:,:) :: ze3f ! 3D workspace (lk_vvl=T)
611	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
612	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
613	REAL(KIND=jprb) :: zhook_handle
614
615	CHARACTER(LEN=*), PARAMETER :: RoutineName='VOR_EEN'
616
617	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
618
619	#else
620	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ze3f ! lk_vvl=F, ze3f=1/e3f saved one for all
621	#endif
622	!!----------------------------------------------------------------------
623	!
624	IF( nn_timing == 1 ) CALL timing_start('vor_een')
625	!
626	CALL wrk_alloc( jpi, jpj, zwx , zwy , zwz )
627	CALL wrk_alloc( jpi, jpj, ztnw, ztne, ztsw, ztse )
628	#if defined key_vvl
629	CALL wrk_alloc( jpi, jpj, jpk, ze3f )
630	#endif
631	!
632	IF( kt == nit000 ) THEN
633	IF(lwp) WRITE(numout,*)
634	IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
635	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
636	#if ! defined key_vvl
637	IF( .NOT.ALLOCATED(ze3f) ) THEN
638	ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr )
639	IF( lk_mpp ) CALL mpp_sum ( ierr )
640	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
641	ENDIF
642	ze3f(:,:,:) = 0._wp
643	#endif
644	ENDIF
645
646	IF( kt == nit000 .OR. lk_vvl ) THEN ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points)
647
648	IF( ln_dynvor_een_old ) THEN ! original formulation
649	DO jk = 1, jpk
650	DO jj = 1, jpjm1
651	DO ji = 1, jpim1
652	ze3 = ( fse3t(ji,jj+1,jk)tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)tmask(ji+1,jj+1,jk) &
653	& + fse3t(ji,jj ,jk)tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)tmask(ji+1,jj ,jk) )
654	IF( ze3 /= 0._wp ) ze3f(ji,jj,jk) = 4.0_wp / ze3
655	END DO
656	END DO
657	END DO
658	ELSE ! new formulation from NEMO 3.6
659	DO jk = 1, jpk
660	DO jj = 1, jpjm1
661	DO ji = 1, jpim1
662	ze3 = ( fse3t(ji,jj+1,jk)tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)tmask(ji+1,jj+1,jk) &
663	& + fse3t(ji,jj ,jk)tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)tmask(ji+1,jj ,jk) )
664	zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
665	& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )
666	IF( ze3 /= 0._wp ) ze3f(ji,jj,jk) = zmsk / ze3
667	END DO
668	END DO
669	END DO
670	ENDIF
671
672	CALL lbc_lnk( ze3f, 'F', 1. )
673	ENDIF
674
675	zfac12 = 1._wp / 12._wp ! Local constant initialization
676
677
678	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
679	! ! ===============
680	DO jk = 1, jpkm1 ! Horizontal slab
681	! ! ===============
682
683	! Potential vorticity and horizontal fluxes
684	! -----------------------------------------
685	SELECT CASE( kvor ) ! vorticity considered
686	CASE ( 1 ) ! planetary vorticity (Coriolis)
687	zwz(:,:) = ff(:,:) * ze3f(:,:,jk)
688	CASE ( 2 ) ! relative vorticity
689	zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
690	CASE ( 3 ) ! metric term
691	DO jj = 1, jpjm1
692	DO ji = 1, fs_jpim1 ! vector opt.
693	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
694	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
695	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk)
696	END DO
697	END DO
698	CALL lbc_lnk( zwz, 'F', 1. )
699	CASE ( 4 ) ! total (relative + planetary vorticity)
700	zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
701	CASE ( 5 ) ! total (coriolis + metric)
702	DO jj = 1, jpjm1
703	DO ji = 1, fs_jpim1 ! vector opt.
704	zwz(ji,jj) = ( ff (ji,jj) &
705	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
706	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
707	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
708	& ) * ze3f(ji,jj,jk)
709	END DO
710	END DO
711	CALL lbc_lnk( zwz, 'F', 1. )
712	END SELECT
713
714	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
715	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
716
717	! Compute and add the vorticity term trend
718	! ----------------------------------------
719	jj = 2
720	ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0
721	DO ji = 2, jpi
722	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
723	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
724	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
725	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
726	END DO
727	DO jj = 3, jpj
728	DO ji = fs_2, jpi ! vector opt. ok because we start at jj = 3
729	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
730	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
731	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
732	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
733	END DO
734	END DO
735	DO jj = 2, jpjm1
736	DO ji = fs_2, fs_jpim1 ! vector opt.
737	zua = + zfac12 / e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
738	& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
739	zva = - zfac12 / e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
740	& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
741	pua(ji,jj,jk) = pua(ji,jj,jk) + zua
742	pva(ji,jj,jk) = pva(ji,jj,jk) + zva
743	END DO
744	END DO
745	! ! ===============
746	END DO ! End of slab
747	! ! ===============
748	CALL wrk_dealloc( jpi, jpj, zwx , zwy , zwz )
749	CALL wrk_dealloc( jpi, jpj, ztnw, ztne, ztsw, ztse )
750	#if defined key_vvl
751	CALL wrk_dealloc( jpi, jpj, jpk, ze3f )
752	#endif
753	!
754	IF( nn_timing == 1 ) CALL timing_stop('vor_een')
755	!
756	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
757	END SUBROUTINE vor_een
758
759
760	SUBROUTINE dyn_vor_init
761	!!---------------------------------------------------------------------
762	!! * ROUTINE dyn_vor_init *
763	!!
764	!! ** Purpose : Control the consistency between cpp options for
765	!! tracer advection schemes
766	!!----------------------------------------------------------------------
767	INTEGER :: ioptio ! local integer
768	INTEGER :: ji, jj, jk ! dummy loop indices
769	INTEGER :: ios ! Local integer output status for namelist read
770	!!
771	NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old
772	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
773	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
774	REAL(KIND=jprb) :: zhook_handle
775
776	CHARACTER(LEN=*), PARAMETER :: RoutineName='DYN_VOR_INIT'
777
778	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
779
780	!!----------------------------------------------------------------------
781
782	REWIND( numnam_ref ) ! Namelist namdyn_vor in reference namelist : Vorticity scheme options
783	READ ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
784	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist', lwp )
785
786	REWIND( numnam_cfg ) ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options
787	READ ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
788	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp )
789	IF(lwm) WRITE ( numond, namdyn_vor )
790
791	IF(lwp) THEN ! Namelist print
792	WRITE(numout,*)
793	WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
794	WRITE(numout,*) '~~~~~~~~~~~~'
795	WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme'
796	WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene
797	WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens
798	WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix
799	WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een
800	WRITE(numout,*) ' enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old
801	ENDIF
802
803	! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
804	! at angles with three ocean points and one land point
805	IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
806	DO jk = 1, jpk
807	DO jj = 2, jpjm1
808	DO ji = 2, jpim1
809	IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) &
810	fmask(ji,jj,jk) = 1._wp
811	END DO
812	END DO
813	END DO
814	!
815	CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
816	!
817	ENDIF
818
819	ioptio = 0 ! Control of vorticity scheme options
820	IF( ln_dynvor_ene ) ioptio = ioptio + 1
821	IF( ln_dynvor_ens ) ioptio = ioptio + 1
822	IF( ln_dynvor_mix ) ioptio = ioptio + 1
823	IF( ln_dynvor_een ) ioptio = ioptio + 1
824	IF( ln_dynvor_een_old ) ioptio = ioptio + 1
825	IF( lk_esopa ) ioptio = 1
826
827	IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
828
829	! ! Set nvor (type of scheme for vorticity)
830	IF( ln_dynvor_ene ) nvor = 0
831	IF( ln_dynvor_ens ) nvor = 1
832	IF( ln_dynvor_mix ) nvor = 2
833	IF( ln_dynvor_een .or. ln_dynvor_een_old ) nvor = 3
834	IF( lk_esopa ) nvor = -1
835
836	! ! Set ncor, nrvm, ntot (type of vorticity)
837	IF(lwp) WRITE(numout,*)
838	ncor = 1
839	IF( ln_dynadv_vec ) THEN
840	IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity'
841	nrvm = 2
842	ntot = 4
843	ELSE
844	IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term'
845	nrvm = 3
846	ntot = 5
847	ENDIF
848
849	IF(lwp) THEN ! Print the choice
850	WRITE(numout,*)
851	IF( nvor == 0 ) WRITE(numout,*) ' vorticity scheme : energy conserving scheme'
852	IF( nvor == 1 ) WRITE(numout,*) ' vorticity scheme : enstrophy conserving scheme'
853	IF( nvor == 2 ) WRITE(numout,*) ' vorticity scheme : mixed enstrophy/energy conserving scheme'
854	IF( nvor == 3 ) WRITE(numout,*) ' vorticity scheme : energy and enstrophy conserving scheme'
855	IF( nvor == -1 ) WRITE(numout,*) ' esopa test: use all lateral physics options'
856	ENDIF
857	!
858	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
859	END SUBROUTINE dyn_vor_init
860
861	!!==============================================================================
862	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: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 11738

Download in other formats: