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.

traadv_qck.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 29.5 KB

Line
1	MODULE traadv_qck
2	!!==============================================================================
3	!! * MODULE traadv_qck *
4	!! Ocean tracers: horizontal & vertical advective trend
5	!!==============================================================================
6	!! History : 3.0 ! 2008-07 (G. Reffray) Original code
7	!! 3.3 ! 2010-05 (C.Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
8	!!----------------------------------------------------------------------
9
10	!!----------------------------------------------------------------------
11	!! tra_adv_qck : update the tracer trend with the horizontal advection
12	!! trends using a 3rd order finite difference scheme
13	!! tra_adv_qck_i : apply QUICK scheme in i-direction
14	!! tra_adv_qck_j : apply QUICK scheme in j-direction
15	!! tra_adv_cen2_k : 2nd centered scheme for the vertical advection
16	!!----------------------------------------------------------------------
17	USE oce ! ocean dynamics and active tracers
18	USE dom_oce ! ocean space and time domain
19	USE trdmod_oce ! ocean space and time domain
20	USE trdtra ! ocean tracers trends
21	USE trabbl ! advective term in the BBL
22	USE lib_mpp ! distribued memory computing
23	USE lbclnk ! ocean lateral boundary condition (or mpp link)
24	USE dynspg_oce ! surface pressure gradient variables
25	USE in_out_manager ! I/O manager
26	USE diaptr ! poleward transport diagnostics
27	USE trc_oce ! share passive tracers/Ocean variables
28
29	IMPLICIT NONE
30	PRIVATE
31
32	PUBLIC tra_adv_qck ! routine called by step.F90
33
34	LOGICAL :: l_trd ! flag to compute trends
35	REAL(wp) :: r1_6 = 1./ 6. ! 1/6 ratio
36
37	!! * Control permutation of array indices
38	# include "oce_ftrans.h90"
39	# include "dom_oce_ftrans.h90"
40	# include "trc_oce_ftrans.h90"
41
42	!! * Substitutions
43	# include "domzgr_substitute.h90"
44	# include "vectopt_loop_substitute.h90"
45	!!----------------------------------------------------------------------
46	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
47	!! $Id$
48	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
49	!!----------------------------------------------------------------------
50	CONTAINS
51
52	SUBROUTINE tra_adv_qck ( kt, cdtype, p2dt, pun, pvn, pwn, &
53	& ptb, ptn, pta, kjpt )
54	!!----------------------------------------------------------------------
55	!! * ROUTINE tra_adv_qck *
56	!!
57	!! ** Purpose : Compute the now trend due to the advection of tracers
58	!! and add it to the general trend of passive tracer equations.
59	!!
60	!! ** Method : The advection is evaluated by a third order scheme
61	!! For a positive velocity u : u(i)>0
62	!! \|--FU--\|--FC--\|--FD--\|------\|
63	!! i-1 i i+1 i+2
64	!!
65	!! For a negative velocity u : u(i)<0
66	!! \|------\|--FD--\|--FC--\|--FU--\|
67	!! i-1 i i+1 i+2
68	!! where FU is the second upwind point
69	!! FD is the first douwning point
70	!! FC is the central point (or the first upwind point)
71	!!
72	!! Flux(i) = u(i) * { 0.5(FC+FD) -0.5C(i)(FD-FC) -((1-C(i))/6)(FU+FD-2FC) }
73	!! with C(i)=\|u(i)\|dx(i)/dt (=Courant number)
74	!!
75	!! dt = 2*rdtra and the scalar values are tb and sb
76	!!
77	!! On the vertical, the simple centered scheme used ptn
78	!!
79	!! The fluxes are bounded by the ULTIMATE limiter to
80	!! guarantee the monotonicity of the solution and to
81	!! prevent the appearance of spurious numerical oscillations
82	!!
83	!! ** Action : - update (pta) with the now advective tracer trends
84	!! - save the trends
85	!!
86	!! ** Reference : Leonard (1979, 1991)
87	!!----------------------------------------------------------------------
88	INTEGER , INTENT(in ) :: kt ! ocean time-step index
89	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
90	INTEGER , INTENT(in ) :: kjpt ! number of tracers
91	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
92
93	!! DCSE_NEMO: This style defeats ftrans
94	! REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
95	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
96	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
97
98	!FTRANS pun pvn pwn :I :I :z
99	!FTRANS ptb ptn :I :I :z :
100	!FTRANS pta :I :I :z :
101	REAL(wp), INTENT(in ) :: pun(jpi,jpj,jpk) ! ocean velocity component (u)
102	REAL(wp), INTENT(in ) :: pvn(jpi,jpj,jpk) ! ocean velocity component (v)
103	REAL(wp), INTENT(in ) :: pwn(jpi,jpj,jpk) ! ocean velocity component (w)
104	REAL(wp), INTENT(in ) :: ptb(jpi,jpj,jpk,kjpt) ! tracer fields (before)
105	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpk,kjpt) ! tracer fields (now)
106	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpk,kjpt) ! tracer trend
107
108	!!----------------------------------------------------------------------
109
110	IF( kt == nit000 ) THEN
111	IF(lwp) WRITE(numout,*)
112	IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
113	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
114	IF(lwp) WRITE(numout,*)
115	!
116	l_trd = .FALSE.
117	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
118	ENDIF
119
120	! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
121	CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt )
122	CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt )
123
124	! II. The vertical fluxes are computed with the 2nd order centered scheme
125	CALL tra_adv_cen2_k( kt, cdtype, pwn, ptn, pta, kjpt )
126	!
127
128	!! * Reset control of array index permutation
129	!FTRANS CLEAR
130	# include "oce_ftrans.h90"
131	# include "dom_oce_ftrans.h90"
132	# include "trc_oce_ftrans.h90"
133
134	END SUBROUTINE tra_adv_qck
135
136
137	SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun, &
138	& ptb, ptn, pta, kjpt )
139	!!----------------------------------------------------------------------
140	!!
141	!!----------------------------------------------------------------------
142	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
143	USE oce , ONLY: zwx => ua ! ua used as workspace
144	USE wrk_nemo, ONLY: zfu => wrk_3d_1 , zfc => wrk_3d_2, zfd => wrk_3d_3 ! 3D workspace
145
146	!! DCSE_NEMO: need additional directives for renamed module variables
147	!FTRANS zwx zfu zfc zfd :I :I :z
148
149	!
150	INTEGER , INTENT(in ) :: kt ! ocean time-step index
151	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
152	INTEGER , INTENT(in ) :: kjpt ! number of tracers
153	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
154
155	!! DCSE_NEMO: This style defeats ftrans
156	! REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun ! i-velocity components
157	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
158	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
159
160	!FTRANS pun :I :I :z
161	!FTRANS ptb ptn pta :I :I :z :
162	REAL(wp), INTENT(in ) :: pun(jpi,jpj,jpk) ! i-velocity component
163	REAL(wp), INTENT(in ) :: ptb(jpi,jpj,jpk,kjpt) ! tracer field (before)
164	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpk,kjpt) ! tracer field (now)
165	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpk,kjpt) ! tracer trend
166
167	!!
168	INTEGER :: ji, jj, jk, jn ! dummy loop indices
169	REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk ! local scalars
170	!----------------------------------------------------------------------
171	!
172	IF( wrk_in_use(3, 1,2,3) ) THEN
173	CALL ctl_stop('tra_adv_qck_i: requested workspace arrays unavailable') ; RETURN
174	ENDIF
175	! ! ===========
176	DO jn = 1, kjpt ! tracer loop
177	! ! ===========
178	zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0
179	zfd(:,:,:) = 0.0 ; zwx(:,:,:) = 0.0
180	!
181	#if defined key_z_first
182	!--- Computation of the upstream and downstream value of the tracer and the mask
183	DO jj = 2, jpjm1
184	DO ji = 2, jpim1
185	DO jk = 1, jpkm1
186	#else
187	DO jk = 1, jpkm1
188	!
189	!--- Computation of the upstream and downstream value of the tracer and the mask
190	DO jj = 2, jpjm1
191	DO ji = fs_2, fs_jpim1 ! vector opt.
192	#endif
193	! Upstream in the x-direction for the tracer
194	zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)
195	! Downstream in the x-direction for the tracer
196	zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)
197	END DO
198	END DO
199	END DO
200	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) ! Lateral boundary conditions
201
202	!
203	! Horizontal advective fluxes
204	! ---------------------------
205	!
206	#if defined key_z_first
207	DO jj = 2, jpjm1
208	DO ji = 2, jpim1
209	DO jk = 1, jpkm1
210	#else
211	DO jk = 1, jpkm1
212	DO jj = 2, jpjm1
213	DO ji = fs_2, fs_jpim1 ! vector opt.
214	#endif
215	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
216	zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T
217	END DO
218	END DO
219	END DO
220	!
221	#if defined key_z_first
222	DO jj = 2, jpjm1
223	DO ji = 2, jpim1
224	DO jk = 1, jpkm1
225	zdt = p2dt(jk)
226	#else
227	DO jk = 1, jpkm1
228	zdt = p2dt(jk)
229	DO jj = 2, jpjm1
230	DO ji = fs_2, fs_jpim1 ! vector opt.
231	#endif
232	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
233	zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)
234	zwx(ji,jj,jk) = ABS( pun(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)
235	zfc(ji,jj,jk) = zdir * ptb(ji ,jj,jk,jn) + ( 1. - zdir ) * ptb(ji+1,jj,jk,jn) ! FC in the x-direction for T
236	zfd(ji,jj,jk) = zdir * ptb(ji+1,jj,jk,jn) + ( 1. - zdir ) * ptb(ji ,jj,jk,jn) ! FD in the x-direction for T
237	END DO
238	END DO
239	END DO
240	!--- Lateral boundary conditions
241	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
242	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
243
244	!--- QUICKEST scheme
245	CALL quickest( zfu, zfd, zfc, zwx )
246	!
247	! Mask at the T-points in the x-direction (mask=0 or mask=1)
248	#if defined key_z_first
249	DO jj = 2, jpjm1
250	DO ji = 2, jpim1
251	DO jk = 1, jpkm1
252	#else
253	DO jk = 1, jpkm1
254	DO jj = 2, jpjm1
255	DO ji = fs_2, fs_jpim1 ! vector opt.
256	#endif
257	zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
258	END DO
259	END DO
260	END DO
261	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ! Lateral boundary conditions
262
263	!
264	! Tracer flux on the x-direction
265	#if defined key_z_first
266	DO jj = 2, jpjm1
267	DO ji = 2, jpim1
268	DO jk = 1, jpkm1
269	#else
270	DO jk = 1, jpkm1
271	DO jj = 2, jpjm1
272	DO ji = fs_2, fs_jpim1 ! vector opt.
273	#endif
274	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
275	!--- If the second ustream point is a land point
276	!--- the flux is computed by the 1st order UPWIND scheme
277	zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
278	zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
279	zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
280	END DO
281	END DO
282	#if defined key_z_first
283	END DO
284	! Computation of the trend
285	DO jj = 2, jpjm1
286	DO ji = 2, jpim1
287	DO jk = 1, jpkm1
288	#else
289	!
290	! Computation of the trend
291	DO jj = 2, jpjm1
292	DO ji = fs_2, fs_jpim1 ! vector opt.
293	#endif
294	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
295	! horizontal advective trends
296	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
297	!--- add it to the general tracer trends
298	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
299	END DO
300	END DO
301	!
302	END DO
303	! ! trend diagnostics (contribution of upstream fluxes)
304	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
305	!
306	END DO
307	!
308	IF( wrk_not_released(3, 1,2,3) ) CALL ctl_stop('tra_adv_qck_i: failed to release workspace arrays')
309	!
310
311	!! * Reset control of array index permutation
312	!FTRANS CLEAR
313	# include "oce_ftrans.h90"
314	# include "dom_oce_ftrans.h90"
315	# include "trc_oce_ftrans.h90"
316
317	END SUBROUTINE tra_adv_qck_i
318
319
320	SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn, &
321	& ptb, ptn, pta, kjpt )
322	!!----------------------------------------------------------------------
323	!!
324	!!----------------------------------------------------------------------
325	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
326	USE oce , ONLY: zwy => ua ! ua used as workspace
327	USE wrk_nemo, ONLY: zfu => wrk_3d_1 , zfc => wrk_3d_2, zfd => wrk_3d_3 ! 3D workspace
328
329	!! DCSE_NEMO: need additional directives for renamed module variables
330	!FTRANS zwy zfu zfc zfd :I :I :z
331
332	!
333	INTEGER , INTENT(in ) :: kt ! ocean time-step index
334	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
335	INTEGER , INTENT(in ) :: kjpt ! number of tracers
336	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
337
338	!! DCSE_NEMO: This style defeats ftrans
339	! REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pvn ! j-velocity components
340	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
341	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
342
343	!FTRANS pvn :I :I :z
344	!FTRANS ptb ptn pta :I :I :z :
345	REAL(wp), INTENT(in ) :: pvn(jpi,jpj,jpk) ! j-velocity component
346	REAL(wp), INTENT(in ) :: ptb(jpi,jpj,jpk,kjpt) ! tracer field (before)
347	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpk,kjpt) ! tracer field (now)
348	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpk,kjpt) ! tracer trend
349
350	!!
351	INTEGER :: ji, jj, jk, jn ! dummy loop indices
352	REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk ! local scalars
353	!----------------------------------------------------------------------
354	!
355	IF(wrk_in_use(3, 1,2,3))THEN
356	CALL ctl_stop('tra_adv_qck_j: ERROR: requested workspace arrays unavailable')
357	RETURN
358	END IF
359	! ! ===========
360	DO jn = 1, kjpt ! tracer loop
361	! ! ===========
362	zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0
363	zfd(:,:,:) = 0.0 ; zwy(:,:,:) = 0.0
364	!
365	#if defined key_z_first
366	!--- Computation of the ustream and downstream value of the tracer and the mask
367	DO jj = 2, jpjm1
368	DO ji = 2, jpim1
369	DO jk = 1, jpkm1
370	#else
371	DO jk = 1, jpkm1
372	!
373	!--- Computation of the ustream and downstream value of the tracer and the mask
374	DO jj = 2, jpjm1
375	DO ji = fs_2, fs_jpim1 ! vector opt.
376	#endif
377	! Upstream in the x-direction for the tracer
378	zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
379	! Downstream in the x-direction for the tracer
380	zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
381	END DO
382	END DO
383	END DO
384	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) ! Lateral boundary conditions
385
386
387	!
388	! Horizontal advective fluxes
389	! ---------------------------
390	!
391	#if defined key_z_first
392	DO jj = 2, jpjm1
393	DO ji = 2, jpim1
394	DO jk = 1, jpkm1
395	#else
396	DO jk = 1, jpkm1
397	DO jj = 2, jpjm1
398	DO ji = fs_2, fs_jpim1 ! vector opt.
399	#endif
400	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
401	zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T
402	END DO
403	END DO
404	END DO
405	!
406	#if defined key_z_first
407	DO jj = 2, jpjm1
408	DO ji = 2, jpim1
409	DO jk = 1, jpkm1
410	zdt = p2dt(jk)
411	#else
412	DO jk = 1, jpkm1
413	zdt = p2dt(jk)
414	DO jj = 2, jpjm1
415	DO ji = fs_2, fs_jpim1 ! vector opt.
416	#endif
417	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
418	zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
419	zwy(ji,jj,jk) = ABS( pvn(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)
420	zfc(ji,jj,jk) = zdir * ptb(ji,jj ,jk,jn) + ( 1. - zdir ) * ptb(ji,jj+1,jk,jn) ! FC in the x-direction for T
421	zfd(ji,jj,jk) = zdir * ptb(ji,jj+1,jk,jn) + ( 1. - zdir ) * ptb(ji,jj ,jk,jn) ! FD in the x-direction for T
422	END DO
423	END DO
424	END DO
425
426	!--- Lateral boundary conditions
427	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
428	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
429
430	!--- QUICKEST scheme
431	CALL quickest( zfu, zfd, zfc, zwy )
432	!
433	! Mask at the T-points in the x-direction (mask=0 or mask=1)
434	#if defined key_z_first
435	DO jj = 2, jpjm1
436	DO ji = 2, jpim1
437	DO jk = 1, jpkm1
438	#else
439	DO jk = 1, jpkm1
440	DO jj = 2, jpjm1
441	DO ji = fs_2, fs_jpim1 ! vector opt.
442	#endif
443	zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
444	END DO
445	END DO
446	END DO
447	!--- Lateral boundary conditions
448	CALL lbc_lnk( zfu(:,:,:), 'T', 1. )
449	!
450	! Tracer flux on the x-direction
451	#if defined key_z_first
452	DO jj = 2, jpjm1
453	DO ji = 2, jpim1
454	DO jk = 1, jpkm1
455	#else
456	DO jk = 1, jpkm1
457	!
458	DO jj = 2, jpjm1
459	DO ji = fs_2, fs_jpim1 ! vector opt.
460	#endif
461	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
462	!--- If the second ustream point is a land point
463	!--- the flux is computed by the 1st order UPWIND scheme
464	zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
465	zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
466	zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
467	END DO
468	END DO
469	#if defined key_z_first
470	END DO
471	! Computation of the trend
472	DO jj = 2, jpjm1
473	DO ji = 2, jpim1
474	DO jk = 1, jpkm1
475	#else
476	!
477	! Computation of the trend
478	DO jj = 2, jpjm1
479	DO ji = fs_2, fs_jpim1 ! vector opt.
480	#endif
481	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
482	! horizontal advective trends
483	ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
484	!--- add it to the general tracer trends
485	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
486	END DO
487	END DO
488	!
489	END DO
490	! ! trend diagnostics (contribution of upstream fluxes)
491	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
492	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
493	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
494	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) )
495	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) )
496	ENDIF
497	!
498	END DO
499	!
500	IF( wrk_not_released(3, 1,2,3) ) CALL ctl_stop('tra_adv_qck_j: failed to release workspace arrays')
501	!
502
503	!! * Reset control of array index permutation
504	!FTRANS CLEAR
505	# include "oce_ftrans.h90"
506	# include "dom_oce_ftrans.h90"
507	# include "trc_oce_ftrans.h90"
508
509	END SUBROUTINE tra_adv_qck_j
510
511
512	SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn, &
513	& ptn, pta, kjpt )
514	!!----------------------------------------------------------------------
515	!!
516	!!----------------------------------------------------------------------
517	USE oce, ONLY: zwz => ua ! ua used as workspace
518
519	!! DCSE_NEMO: need additional directives for renamed module variables
520	!FTRANS zwz :I :I :z
521
522	!
523	INTEGER , INTENT(in ) :: kt ! ocean time-step index
524	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
525	INTEGER , INTENT(in ) :: kjpt ! number of tracers
526
527	!! DCSE_NEMO: This style defeats ftrans
528	! REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pwn ! vertical velocity
529	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptn ! tracer fields (now)
530	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
531
532	!FTRANS pwn :I :I :z
533	!FTRANS ptn pta :I :I :z :
534	REAL(wp), INTENT(in ) :: pwn(jpi,jpj,jpk) ! vertical velocity
535	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpk,kjpt) ! tracer fields (now)
536	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpk,kjpt) ! tracer trend
537
538	!
539	INTEGER :: ji, jj, jk, jn ! dummy loop indices
540	REAL(wp) :: zbtr , ztra ! local scalars
541	!!----------------------------------------------------------------------
542
543	! ! ===========
544	DO jn = 1, kjpt ! tracer loop
545	! ! ===========
546	! 1. Bottom value : flux set to zero
547	zwz(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero
548	!
549	! ! Surface value
550	IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! Variable volume : flux set to zero
551	ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn) ! Constant volume : advective flux through the surface
552	ENDIF
553	!
554	#if defined key_z_first
555	DO jj = 2, jpjm1
556	DO ji = 2, jpim1
557	DO jk = 2, jpkm1 ! Interior point: second order centered tracer flux at w-point
558	#else
559	DO jk = 2, jpkm1 ! Interior point: second order centered tracer flux at w-point
560	DO jj = 2, jpjm1
561	DO ji = fs_2, fs_jpim1 ! vector opt.
562	#endif
563	zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
564	END DO
565	END DO
566	END DO
567	!
568	#if defined key_z_first
569	DO jj = 2, jpjm1
570	DO ji = 2, jpim1
571	DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==!
572	#else
573	DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==!
574	DO jj = 2, jpjm1
575	DO ji = fs_2, fs_jpim1 ! vector opt.
576	#endif
577	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
578	! k- vertical advective trends
579	ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )
580	! added to the general tracer trends
581	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
582	END DO
583	END DO
584	END DO
585	! ! Save the vertical advective trends for diagnostic
586	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) )
587	!
588	END DO
589	!
590
591	!! * Reset control of array index permutation
592	!FTRANS CLEAR
593	# include "oce_ftrans.h90"
594	# include "dom_oce_ftrans.h90"
595	# include "trc_oce_ftrans.h90"
596
597	END SUBROUTINE tra_adv_cen2_k
598
599
600	SUBROUTINE quickest( pfu, pfd, pfc, puc )
601	!!----------------------------------------------------------------------
602	!!
603	!! ** Purpose : Computation of advective flux with Quickest scheme
604	!!
605	!! ** Method :
606	!!----------------------------------------------------------------------
607
608	!! DCSE_NEMO: This style defeats ftrans
609
610	! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfu ! second upwind point
611	! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfd ! first douwning point
612	! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfc ! the central point (or the first upwind point)
613	! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: puc ! input as Courant number ; output as flux
614
615	!FTRANS pfu pfd pfc puc :I :I :z
616	REAL(wp), INTENT(in ) :: pfu(jpi,jpj,jpk) ! second upwind point
617	REAL(wp), INTENT(in ) :: pfd(jpi,jpj,jpk) ! first douwning point
618	REAL(wp), INTENT(in ) :: pfc(jpi,jpj,jpk) ! the central point (or the first upwind point)
619	REAL(wp), INTENT(inout) :: puc(jpi,jpj,jpk) ! input as Courant number ; output as flux
620
621	!!
622	INTEGER :: ji, jj, jk ! dummy loop indices
623	REAL(wp) :: zcoef1, zcoef2, zcoef3 ! local scalars
624	REAL(wp) :: zc, zcurv, zfho ! - -
625	!----------------------------------------------------------------------
626
627	#if defined key_z_first
628	DO jj = 1, jpj
629	DO ji = 1, jpi
630	DO jk = 1, jpkm1
631	#else
632	DO jk = 1, jpkm1
633	DO jj = 1, jpj
634	DO ji = 1, jpi
635	#endif
636	zc = puc(ji,jj,jk) ! Courant number
637	zcurv = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
638	zcoef1 = 0.5 * ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
639	zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
640	zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
641	zfho = zcoef1 - zcoef2 - zcoef3 ! phi_f QUICKEST
642	!
643	zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
644	zcoef2 = ABS( zcoef1 )
645	zcoef3 = ABS( zcurv )
646	IF( zcoef3 >= zcoef2 ) THEN
647	zfho = pfc(ji,jj,jk)
648	ELSE
649	zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) ) ! phi_REF
650	IF( zcoef1 >= 0. ) THEN
651	zfho = MAX( pfc(ji,jj,jk), zfho )
652	zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) )
653	ELSE
654	zfho = MIN( pfc(ji,jj,jk), zfho )
655	zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) )
656	ENDIF
657	ENDIF
658	puc(ji,jj,jk) = zfho
659	END DO
660	END DO
661	END DO
662	!
663	END SUBROUTINE quickest
664
665	!!======================================================================
666	END MODULE traadv_qck

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

Download in other formats: