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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynvor.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/DYN/dynvor.F90 @ 1104

Last change on this file since 1104 was 1104, checked in by ctlod, 16 years ago
trunk: change the second ztrdu(:,:,:) field into ztrdv(:,:,:) when calling trd_mod subroutine for the planetary vorticity trend in dynvor.F90, see ticket: #20
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 37.0 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 : 1.0 ! 89-12 (P. Andrich) vor_ens: Original code
8	!! 5.0 ! 91-11 (G. Madec) vor_ene, vor_mix: Original code
9	!! 6.0 ! 96-01 (G. Madec) s-coord, suppress work arrays
10	!! 8.5 ! 02-08 (G. Madec) F90: Free form and module
11	!! 8.5 ! 04-02 (G. Madec) vor_een: Original code
12	!! 9.0 ! 03-08 (G. Madec) vor_ctl: Original code
13	!! 9.0 ! 05-11 (G. Madec) dyn_vor: Original code (new step architecture)
14	!! 9.0 ! 06-11 (G. Madec) flux form advection: add metric term
15	!!----------------------------------------------------------------------
16
17	!!----------------------------------------------------------------------
18	!! dyn_vor : Update the momentum trend with the vorticity trend
19	!! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T)
20	!! vor_ene : energy conserving scheme (ln_dynvor_ene=T)
21	!! vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=T)
22	!! vor_een : energy and enstrophy conserving (ln_dynvor_een=T)
23	!! vor_ctl : set and control of the different vorticity option
24	!!----------------------------------------------------------------------
25	USE oce ! ocean dynamics and tracers
26	USE dom_oce ! ocean space and time domain
27	USE dynadv ! momentum advection (use ln_dynadv_vec value)
28	USE trdmod ! ocean dynamics trends
29	USE trdmod_oce ! ocean variables trends
30	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
31	USE prtctl ! Print control
32	USE in_out_manager ! I/O manager
33
34	IMPLICIT NONE
35	PRIVATE
36
37	PUBLIC dyn_vor ! routine called by step.F90
38
39	!!* Namelist nam_dynvor: vorticity term
40	LOGICAL, PUBLIC :: ln_dynvor_ene = .FALSE. !: energy conserving scheme
41	LOGICAL, PUBLIC :: ln_dynvor_ens = .TRUE. !: enstrophy conserving scheme
42	LOGICAL, PUBLIC :: ln_dynvor_mix = .FALSE. !: mixed scheme
43	LOGICAL, PUBLIC :: ln_dynvor_een = .FALSE. !: energy and enstrophy conserving scheme
44
45	INTEGER :: nvor = 0 ! type of vorticity trend used
46	INTEGER :: ncor = 1 ! coriolis
47	INTEGER :: nrvm = 2 ! =2 relative vorticity ; =3 metric term
48	INTEGER :: ntot = 4 ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term
49
50	!! * Substitutions
51	# include "domzgr_substitute.h90"
52	# include "vectopt_loop_substitute.h90"
53	!!----------------------------------------------------------------------
54	!! OPA 9.0 , LOCEAN-IPSL (2006)
55	!! $Header$
56	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
57	!!----------------------------------------------------------------------
58
59	CONTAINS
60
61	SUBROUTINE dyn_vor( kt )
62	!!----------------------------------------------------------------------
63	!!
64	!! ** Purpose : compute the lateral ocean tracer physics.
65	!!
66	!! ** Action : - Update (ua,va) with the now vorticity term trend
67	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
68	!! and planetary vorticity trends) ('key_trddyn')
69	!!----------------------------------------------------------------------
70	USE oce, ONLY : ztrdu => ta ! use ta as 3D workspace
71	USE oce, ONLY : ztrdv => sa ! use sa as 3D workspace
72	!!
73	INTEGER, INTENT( in ) :: kt ! ocean time-step index
74	!!----------------------------------------------------------------------
75
76	IF( kt == nit000 ) CALL vor_ctl ! initialisation & control of options
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	INTEGER , INTENT(in ) :: kt ! ocean time-step index
207	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
208	! ! =nrvm (relative vorticity or metric)
209	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
210	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
211	!!
212	INTEGER :: ji, jj, jk ! dummy loop indices
213	REAL(wp) :: zx1, zy1, zfact2 ! temporary scalars
214	REAL(wp) :: zx2, zy2 ! " "
215	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! temporary 2D workspace
216	!!----------------------------------------------------------------------
217
218	IF( kt == nit000 ) THEN
219	IF(lwp) WRITE(numout,*)
220	IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
221	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
222	ENDIF
223
224	! Local constant initialization
225	zfact2 = 0.5 * 0.5
226
227	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
228	! ! ===============
229	DO jk = 1, jpkm1 ! Horizontal slab
230	! ! ===============
231	! Potential vorticity and horizontal fluxes
232	! -----------------------------------------
233	SELECT CASE( kvor ) ! vorticity considered
234	CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis)
235	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity
236	CASE ( 3 ) ! metric term
237	DO jj = 1, jpjm1
238	DO ji = 1, fs_jpim1 ! vector opt.
239	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
240	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
241	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
242	END DO
243	END DO
244	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity)
245	CASE ( 5 ) ! total (coriolis + metric)
246	DO jj = 1, jpjm1
247	DO ji = 1, fs_jpim1 ! vector opt.
248	zwz(ji,jj) = ( ff (ji,jj) &
249	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
250	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
251	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
252	& )
253	END DO
254	END DO
255	END SELECT
256
257	IF( ln_sco ) THEN
258	zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
259	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
260	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
261	ELSE
262	zwx(:,:) = e2u(:,:) * un(:,:,jk)
263	zwy(:,:) = e1v(:,:) * vn(:,:,jk)
264	ENDIF
265
266	! Compute and add the vorticity term trend
267	! ----------------------------------------
268	DO jj = 2, jpjm1
269	DO ji = fs_2, fs_jpim1 ! vector opt.
270	zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
271	zy2 = zwy(ji,jj ) + zwy(ji+1,jj )
272	zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
273	zx2 = zwx(ji ,jj) + zwx(ji ,jj+1)
274	pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
275	pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
276	END DO
277	END DO
278	! ! ===============
279	END DO ! End of slab
280	! ! ===============
281	END SUBROUTINE vor_ene
282
283
284	SUBROUTINE vor_mix( kt )
285	!!----------------------------------------------------------------------
286	!! * ROUTINE vor_mix *
287	!!
288	!! ** Purpose : Compute the now total vorticity trend and add it to
289	!! the general trend of the momentum equation.
290	!!
291	!! ** Method : Trend evaluated using now fields (centered in time)
292	!! Mixte formulation : conserves the potential enstrophy of a hori-
293	!! zontally non-divergent flow for (rotzu x uh), the relative vor-
294	!! ticity term and the horizontal kinetic energy for (f x uh), the
295	!! coriolis term. the now trend of the vorticity term is given by:
296	!! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
297	!! voru = 1/e1u mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]
298	!! +1/e1u mj-1[ f/e3f mi(e1v*e3v vn) ]
299	!! vorv = 1/e2v mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]
300	!! +1/e2v mi-1[ f/e3f mj(e2u*e3u un) ]
301	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
302	!! voru = 1/e1u mj-1(rotn) mj-1[ mi(e1v vn) ]
303	!! +1/e1u mj-1[ f mi(e1v vn) ]
304	!! vorv = 1/e2v mi-1(rotn) mi-1[ mj(e2u un) ]
305	!! +1/e2v mi-1[ f mj(e2u un) ]
306	!! Add this now trend to the general momentum trend (ua,va):
307	!! (ua,va) = (ua,va) + ( voru , vorv )
308	!!
309	!! ** Action : - Update (ua,va) arrays with the now vorticity term trend
310	!! - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
311	!! and planetary vorticity trends) ('key_trddyn')
312	!!
313	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
314	!!----------------------------------------------------------------------
315	INTEGER, INTENT(in) :: kt ! ocean timestep index
316	!!
317	INTEGER :: ji, jj, jk ! dummy loop indices
318	REAL(wp) :: zfact1, zua, zcua, zx1, zy1 ! temporary scalars
319	REAL(wp) :: zfact2, zva, zcva, zx2, zy2 ! " "
320	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz, zww ! temporary 3D workspace
321	!!----------------------------------------------------------------------
322
323	IF( kt == nit000 ) THEN
324	IF(lwp) WRITE(numout,*)
325	IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
326	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
327	ENDIF
328
329	! Local constant initialization
330	zfact1 = 0.5 * 0.25
331	zfact2 = 0.5 * 0.5
332
333	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
334	! ! ===============
335	DO jk = 1, jpkm1 ! Horizontal slab
336	! ! ===============
337
338	! Relative and planetary potential vorticity and horizontal fluxes
339	! ----------------------------------------------------------------
340	IF( ln_sco ) THEN
341	IF( ln_dynadv_vec ) THEN
342	zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
343	ELSE
344	DO jj = 1, jpjm1
345	DO ji = 1, fs_jpim1 ! vector opt.
346	zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
347	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
348	& * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
349	END DO
350	END DO
351	ENDIF
352	zwz(:,:) = ff (:,:) / fse3f(:,:,jk)
353	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
354	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
355	ELSE
356	IF( ln_dynadv_vec ) THEN
357	zww(:,:) = rotn(:,:,jk)
358	ELSE
359	DO jj = 1, jpjm1
360	DO ji = 1, fs_jpim1 ! vector opt.
361	zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
362	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
363	& * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
364	END DO
365	END DO
366	ENDIF
367	zwz(:,:) = ff (:,:)
368	zwx(:,:) = e2u(:,:) * un(:,:,jk)
369	zwy(:,:) = e1v(:,:) * vn(:,:,jk)
370	ENDIF
371
372	! Compute and add the vorticity term trend
373	! ----------------------------------------
374	DO jj = 2, jpjm1
375	DO ji = fs_2, fs_jpim1 ! vector opt.
376	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
377	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
378	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
379	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
380	! enstrophy conserving formulation for relative vorticity term
381	zua = zfact1 * ( zww(ji ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
382	zva =-zfact1 * ( zww(ji-1,jj ) + zww(ji,jj) ) * ( zx1 + zx2 )
383	! energy conserving formulation for planetary vorticity term
384	zcua = zfact2 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
385	zcva =-zfact2 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
386	! mixed vorticity trend added to the momentum trends
387	ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
388	va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
389	END DO
390	END DO
391	! ! ===============
392	END DO ! End of slab
393	! ! ===============
394	END SUBROUTINE vor_mix
395
396
397	SUBROUTINE vor_ens( kt, kvor, pua, pva )
398	!!----------------------------------------------------------------------
399	!! * ROUTINE vor_ens *
400	!!
401	!! ** Purpose : Compute the now total vorticity trend and add it to
402	!! the general trend of the momentum equation.
403	!!
404	!! ** Method : Trend evaluated using now fields (centered in time)
405	!! and the Sadourny (1975) flux FORM formulation : conserves the
406	!! potential enstrophy of a horizontally non-divergent flow. the
407	!! trend of the vorticity term is given by:
408	!! * s-coordinate (ln_sco=T), the e3. are inside the derivative:
409	!! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ]
410	!! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ]
411	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
412	!! voru = 1/e1u mj-1[ rotn+f ] mj-1[ mi(e1v vn) ]
413	!! vorv = 1/e2v mi-1[ rotn+f ] mi-1[ mj(e2u un) ]
414	!! Add this trend to the general momentum trend (ua,va):
415	!! (ua,va) = (ua,va) + ( voru , vorv )
416	!!
417	!! ** Action : - Update (ua,va) arrays with the now vorticity term trend
418	!! - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
419	!! and planetary vorticity trends) ('key_trddyn')
420	!!
421	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
422	!!----------------------------------------------------------------------
423	INTEGER , INTENT(in ) :: kt ! ocean time-step index
424	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
425	! ! =nrvm (relative vorticity or metric)
426	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
427	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
428	!!
429	INTEGER :: ji, jj, jk ! dummy loop indices
430	REAL(wp) :: zfact1, zuav, zvau ! temporary scalars
431	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! temporary 3D workspace
432	!!----------------------------------------------------------------------
433
434	IF( kt == nit000 ) THEN
435	IF(lwp) WRITE(numout,*)
436	IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
437	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
438	ENDIF
439
440	! Local constant initialization
441	zfact1 = 0.5 * 0.25
442
443	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
444	! ! ===============
445	DO jk = 1, jpkm1 ! Horizontal slab
446	! ! ===============
447	! Potential vorticity and horizontal fluxes
448	! -----------------------------------------
449	SELECT CASE( kvor ) ! vorticity considered
450	CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis)
451	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity
452	CASE ( 3 ) ! metric term
453	DO jj = 1, jpjm1
454	DO ji = 1, fs_jpim1 ! vector opt.
455	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
456	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
457	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
458	END DO
459	END DO
460	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity)
461	CASE ( 5 ) ! total (coriolis + metric)
462	DO jj = 1, jpjm1
463	DO ji = 1, fs_jpim1 ! vector opt.
464	zwz(ji,jj) = ( ff (ji,jj) &
465	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
466	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
467	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
468	& )
469	END DO
470	END DO
471	END SELECT
472
473	IF( ln_sco ) THEN
474	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
475	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
476	zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
477	zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
478	zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
479	END DO
480	END DO
481	ELSE
482	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
483	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
484	zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
485	zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
486	END DO
487	END DO
488	ENDIF
489
490	! Compute and add the vorticity term trend
491	! ----------------------------------------
492	DO jj = 2, jpjm1
493	DO ji = fs_2, fs_jpim1 ! vector opt.
494	zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
495	& + zwy(ji ,jj ) + zwy(ji+1,jj ) )
496	zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
497	& + zwx(ji ,jj ) + zwx(ji ,jj+1) )
498	pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) )
499	pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) )
500	END DO
501	END DO
502	! ! ===============
503	END DO ! End of slab
504	! ! ===============
505	END SUBROUTINE vor_ens
506
507
508	SUBROUTINE vor_een( kt, kvor, pua, pva )
509	!!----------------------------------------------------------------------
510	!! * ROUTINE vor_een *
511	!!
512	!! ** Purpose : Compute the now total vorticity trend and add it to
513	!! the general trend of the momentum equation.
514	!!
515	!! ** Method : Trend evaluated using now fields (centered in time)
516	!! and the Arakawa and Lamb (19XX) flux form formulation : conserves
517	!! both the horizontal kinetic energy and the potential enstrophy
518	!! when horizontal divergence is zero.
519	!! The trend of the vorticity term is given by:
520	!! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
521	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
522	!! Add this trend to the general momentum trend (ua,va):
523	!! (ua,va) = (ua,va) + ( voru , vorv )
524	!!
525	!! ** Action : - Update (ua,va) with the now vorticity term trend
526	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
527	!! and planetary vorticity trends) ('key_trddyn')
528	!!
529	!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
530	!!----------------------------------------------------------------------
531	INTEGER , INTENT(in ) :: kt ! ocean time-step index
532	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
533	! ! =nrvm (relative vorticity or metric)
534	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
535	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
536	!!
537	INTEGER :: ji, jj, jk ! dummy loop indices
538	REAL(wp) :: zfac12, zua, zva ! temporary scalars
539	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! temporary 2D workspace
540	REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse ! temporary 3D workspace
541	REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE :: ze3f
542	!!----------------------------------------------------------------------
543
544	IF( kt == nit000 ) THEN
545	IF(lwp) WRITE(numout,*)
546	IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
547	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
548
549	DO jk = 1, jpk
550	DO jj = 1, jpjm1
551	DO ji = 1, jpim1
552	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) &
553	& + fse3t(ji,jj ,jk)tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)tmask(ji+1,jj ,jk) ) * 0.25
554	IF( ze3f(ji,jj,jk) /= 0.e0 ) ze3f(ji,jj,jk) = 1.e0 / ze3f(ji,jj,jk)
555	END DO
556	END DO
557	END DO
558	CALL lbc_lnk( ze3f, 'F', 1. )
559	ENDIF
560
561	! Local constant initialization
562	zfac12 = 1.e0 / 12.e0
563
564
565	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
566	! ! ===============
567	DO jk = 1, jpkm1 ! Horizontal slab
568	! ! ===============
569
570	! Potential vorticity and horizontal fluxes
571	! -----------------------------------------
572	SELECT CASE( kvor ) ! vorticity considered
573	CASE ( 1 ) ; zwz(:,:) = ff(:,:) * ze3f(:,:,jk) ! planetary vorticity (Coriolis)
574	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) ! relative vorticity
575	CASE ( 3 ) ! metric term
576	DO jj = 1, jpjm1
577	DO ji = 1, fs_jpim1 ! vector opt.
578	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
579	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
580	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk)
581	END DO
582	END DO
583	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) ! total (relative + planetary vorticity)
584	CASE ( 5 ) ! total (coriolis + metric)
585	DO jj = 1, jpjm1
586	DO ji = 1, fs_jpim1 ! vector opt.
587	zwz(ji,jj) = ( ff (ji,jj) &
588	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
589	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
590	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
591	& ) * ze3f(ji,jj,jk)
592	END DO
593	END DO
594	END SELECT
595
596	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
597	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
598
599	! Compute and add the vorticity term trend
600	! ----------------------------------------
601	jj=2
602	ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0
603	DO ji = 2, jpi
604	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
605	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
606	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
607	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
608	END DO
609	DO jj = 3, jpj
610	DO ji = fs_2, jpi ! vector opt.
611	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
612	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
613	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
614	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
615	END DO
616	END DO
617	DO jj = 2, jpjm1
618	DO ji = fs_2, fs_jpim1 ! vector opt.
619	zua = + zfac12 / e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
620	& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
621	zva = - zfac12 / e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
622	& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
623	pua(ji,jj,jk) = pua(ji,jj,jk) + zua
624	pva(ji,jj,jk) = pva(ji,jj,jk) + zva
625	END DO
626	END DO
627	! ! ===============
628	END DO ! End of slab
629	! ! ===============
630	END SUBROUTINE vor_een
631
632
633	SUBROUTINE vor_ctl
634	!!---------------------------------------------------------------------
635	!! * ROUTINE vor_ctl *
636	!!
637	!! ** Purpose : Control the consistency between cpp options for
638	!! tracer advection schemes
639	!!----------------------------------------------------------------------
640	INTEGER :: ioptio ! temporary integer
641	NAMELIST/nam_dynvor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een
642	!!----------------------------------------------------------------------
643
644	REWIND ( numnam ) ! Read Namelist nam_dynvor : Vorticity scheme options
645	READ ( numnam, nam_dynvor )
646
647	IF(lwp) THEN ! Namelist print
648	WRITE(numout,*)
649	WRITE(numout,*) 'dyn:vor_ctl : vorticity term : read namelist and control the consistency'
650	WRITE(numout,*) '~~~~~~~~~~~'
651	WRITE(numout,*) ' Namelist nam_dynvor : oice of the vorticity term scheme'
652	WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene
653	WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens
654	WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix
655	WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een
656	ENDIF
657
658	ioptio = 0 ! Control of vorticity scheme options
659	IF( ln_dynvor_ene ) ioptio = ioptio + 1
660	IF( ln_dynvor_ens ) ioptio = ioptio + 1
661	IF( ln_dynvor_mix ) ioptio = ioptio + 1
662	IF( ln_dynvor_een ) ioptio = ioptio + 1
663	IF( lk_esopa ) ioptio = 1
664
665	IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
666
667	! ! Set nvor (type of scheme for vorticity)
668	IF( ln_dynvor_ene ) nvor = 0
669	IF( ln_dynvor_ens ) nvor = 1
670	IF( ln_dynvor_mix ) nvor = 2
671	IF( ln_dynvor_een ) nvor = 3
672	IF( lk_esopa ) nvor = -1
673
674	! ! Set ncor, nrvm, ntot (type of vorticity)
675	IF(lwp) WRITE(numout,*)
676	ncor = 1
677	IF( ln_dynadv_vec ) THEN
678	IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity'
679	nrvm = 2
680	ntot = 4
681	ELSE
682	IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term'
683	nrvm = 3
684	ntot = 5
685	ENDIF
686
687	IF(lwp) THEN ! Print the choice
688	WRITE(numout,*)
689	IF( nvor == 0 ) WRITE(numout,*) ' vorticity scheme : energy conserving scheme'
690	IF( nvor == 1 ) WRITE(numout,*) ' vorticity scheme : enstrophy conserving scheme'
691	IF( nvor == 2 ) WRITE(numout,*) ' vorticity scheme : mixed enstrophy/energy conserving scheme'
692	IF( nvor == 3 ) WRITE(numout,*) ' vorticity scheme : energy and enstrophy conserving scheme'
693	IF( nvor == -1 ) WRITE(numout,*) ' esopa test: use all lateral physics options'
694	ENDIF
695	!
696	END SUBROUTINE vor_ctl
697
698	!!==============================================================================
699	END MODULE dynvor

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

Download in other formats: