Context Navigation

dynvor.F90 @ 2690

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

Download in other formats: