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 trunk/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 703

Last change on this file since 703 was 699, checked in by smasson, 17 years ago
insert revision Id
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 26.1 KB

Line
1	MODULE traadv_qck
2	!!==============================================================================
3	!! * MODULE traadv_qck *
4	!! Ocean active tracers: horizontal & vertical advective trend
5	!!==============================================================================
6
7	!!----------------------------------------------------------------------
8	!! tra_adv_qck : update the tracer trend with the horizontal
9	!! advection trends using a 3st order
10	!! finite difference scheme
11	!! The vertical advection scheme is the 2nd centered scheme
12	!!----------------------------------------------------------------------
13	!! * Modules used
14	USE oce ! ocean dynamics and active tracers
15	USE dom_oce ! ocean space and time domain
16	USE trdmod ! ocean active tracers trends
17	USE trdmod_oce ! ocean variables trends
18	USE flxrnf !
19	USE trabbl ! advective term in the BBL
20	USE ocfzpt !
21	USE lib_mpp
22	USE lbclnk ! ocean lateral boundary condition (or mpp link)
23	USE in_out_manager ! I/O manager
24	USE diaptr ! poleward transport diagnostics
25	USE dynspg_oce !
26	USE prtctl ! Print control
27
28	IMPLICIT NONE
29	PRIVATE
30
31	!! * Accessibility
32	PUBLIC tra_adv_qck ! routine called by step.F90
33
34	!! * Module variables
35	REAL(wp), DIMENSION(jpi,jpj), SAVE :: &
36	zbtr2
37	REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE :: &
38	sl
39	REAL(wp) :: &
40	cst1, cst2, dt, coef1 ! temporary scalars
41	INTEGER :: &
42	ji, jj, jk ! dummy loop indices
43	!!----------------------------------------------------------------------
44	!! * Substitutions
45	# include "domzgr_substitute.h90"
46	# include "vectopt_loop_substitute.h90"
47	!!----------------------------------------------------------------------
48	!! OPA 9.0 , LOCEAN-IPSL (2005)
49	!! $Id$
50	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
51	!!----------------------------------------------------------------------
52
53	CONTAINS
54
55	#if ! defined key_mpp_omp
56	!!----------------------------------------------------------------------
57	!! Default option : quickest advection scheme (k-j-i loop)
58	!!----------------------------------------------------------------------
59
60	SUBROUTINE tra_adv_qck( kt, pun, pvn, pwn )
61	!!----------------------------------------------------------------------
62	!! * ROUTINE tra_adv_qck *
63	!!
64	!! ** Purpose : Compute the now trend due to the advection of tracers
65	!! and add it to the general trend of passive tracer equations.
66	!!
67	!! ** Method : The advection is evaluated by a third order scheme
68	!! For a positive velocity u :
69	!!
70	!!
71	!! i-1 i i+1 i+2
72	!!
73	!! \|--FU--\|--FC--\|--FD--\|------\|
74	!! u(i)>0
75	!!
76	!! For a negative velocity u :
77	!!
78	!! \|------\|--FD--\|--FC--\|--FU--\|
79	!! u(i)<0
80	!!
81	!! FU is the second upwind point
82	!! FD is the first douwning point
83	!! FC is the central point (or the first upwind point)
84	!!
85	!! Flux(i) = u(i) * {0.5(FC+FD) -0.5C(i)(FD-FC) -((1-C(i)å?)/6)(FU+FD-2FC)}
86	!! with C(i)=\|u(i)\|dx(i)/dt (Courant number)
87	!!
88	!! dt = 2*rdtra and the scalar values are tb and sb
89	!!
90	!! On the vertical, the simple centered scheme used tn and sn
91	!!
92	!! The fluxes are bounded by the ULTIMATE limiter
93	!! to guarantee the monotonicity of the solution and to
94	!! prevent the appearance of spurious numerical oscillations
95	!!
96	!! ** Action : - update (ta,sa) with the now advective tracer trends
97	!! - save the trends in (ttrdh,strdh) ('key_trdtra')
98	!!
99	!! ** Reference : Leonard (1979, 1991)
100	!! History :
101	!! 9.0 ! 06-09 (G. Reffray) Original code
102	!!----------------------------------------------------------------------
103	!! * Arguments
104	INTEGER, INTENT( in ) :: kt ! ocean time-step index
105	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! effective ocean velocity, u_component
106	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! effective ocean velocity, v_component
107	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! effective ocean velocity, w_component
108	!!
109	REAL(wp) :: z2 ! temporary scalar
110	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztrdt, ztrds ! temporary 3D workspace
111
112	IF( kt == nit000 ) THEN
113	IF(lwp) WRITE(numout,*)
114	IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3st order quickest advection scheme'
115	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ Vector optimization case'
116	IF(lwp) WRITE(numout,*)
117
118	zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
119	cst1 = 1./12.
120	cst2 = 2./3.
121	IF (l_trdtra ) THEN
122	CALL ctl_warn( ' Trends not yet implemented for PPM advection scheme ' )
123	ENDIF
124	ENDIF
125
126	IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1.
127	ELSE ; z2 = 2.
128	ENDIF
129
130	! Save ta and sa trends
131	IF( l_trdtra ) THEN ! to be done
132	ztrdt(:,:,:) = ta(:,:,:)
133	ztrds(:,:,:) = sa(:,:,:)
134	l_adv = 'qst'
135	ENDIF
136
137	! I. Slope estimation at the T-point for the limiter ULTIMATE
138	! SL = Sum(1/C_out) with C_out : Courant number for the outflows
139	!---------------------------------------------------------------
140
141	sl(:,:,:) = 100.
142
143	! ===============
144	DO jk = 1, jpkm1 ! Horizontal slab
145	! ! ===============
146	dt = z2 * rdttra(jk)
147	DO jj = 2, jpjm1
148	DO ji = 2, fs_jpim1 ! vector opt.
149	coef1 = 1.e-12
150	IF (pun(ji-1,jj ,jk ).LT.0.) coef1 = coef1 + ABS(pun(ji-1,jj ,jk ))*dt/e1t(ji,jj)
151	IF (pun(ji ,jj ,jk ).GT.0.) coef1 = coef1 + ABS(pun(ji ,jj ,jk ))*dt/e1t(ji,jj)
152	IF (pvn(ji ,jj-1,jk ).LT.0.) coef1 = coef1 + ABS(pvn(ji ,jj-1,jk ))*dt/e2t(ji,jj)
153	IF (pvn(ji ,jj ,jk ).GT.0.) coef1 = coef1 + ABS(pvn(ji ,jj ,jk ))*dt/e2t(ji,jj)
154	IF (pwn(ji ,jj ,jk+1).LT.0.) coef1 = coef1 + ABS(pwn(ji ,jj ,jk+1))*dt/(fse3t(ji,jj,jk))
155	IF (pwn(ji ,jj ,jk ).GT.0.) coef1 = coef1 + ABS(pwn(ji ,jj ,jk ))*dt/(fse3t(ji,jj,jk))
156	sl(ji,jj,jk) = 1./coef1
157	sl(ji,jj,jk) = MIN(100.,sl(ji,jj,jk))
158	sl(ji,jj,jk) = MAX(1. ,sl(ji,jj,jk))
159	ENDDO
160	ENDDO
161	ENDDO
162	!--- Lateral boundary conditions on the limiter slope
163	CALL lbc_lnk( sl(:,:,:), 'T', 1. )
164
165	! II. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
166	!---------------------------------------------------------------------------
167
168	CALL tra_adv_qck_hor( kt , pun, pvn, tb , ta , pht_adv , z2)
169	CALL tra_adv_qck_hor( kt , pun, pvn, sb , sa , pst_adv , z2)
170
171	! Save the horizontal advective trends for diagnostic
172	! ---------------------------------------------------
173	! IF( l_trdtra ) THEN ! to be done
174	! ! T/S ZONAL advection trends
175	! ENDIF
176
177	IF(ln_ctl) THEN
178	CALL prt_ctl(tab3d_1=ta, clinfo1=' centered2 had - Ta: ', mask1=tmask, &
179	& tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
180	ENDIF
181
182	! III. The vertical fluxes are computed with the 2nd order centered scheme
183	!-------------------------------------------------------------------------
184
185	CALL tra_adv_qck_ver( pwn, tn , ta, z2 )
186	CALL tra_adv_qck_ver( pwn, sn , sa, z2 )
187
188	! Save the vertical advective trends for diagnostic
189	! -------------------------------------------------
190	! IF( l_trdtra ) THEN ! to be done
191	! Recompute the vertical advection zta & zsa trends computed
192	! at the step 2. above in making the difference between the new
193	! trends and the previous one: ta()/sa - ztdta()/ztdsa() and substract
194	! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
195	! ENDIF
196
197	IF(ln_ctl) THEN
198	CALL prt_ctl(tab3d_1=ta, clinfo1=' centered2 zad - Ta: ', mask1=tmask, &
199	& tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
200	ENDIF
201
202	END SUBROUTINE tra_adv_qck
203
204	SUBROUTINE tra_adv_qck_hor ( kt , pun, pvn, tra , traa , phtra_adv ,z2 )
205	!!----------------------------------------------------------------------
206	!!
207	!!----------------------------------------------------------------------
208	!! * Arguments
209	INTEGER, INTENT( in ) :: kt ! ocean time-step index
210	REAL, INTENT( in ) :: z2
211	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun, pvn ! horizontal effective velocity
212
213	REAL(wp), INTENT ( out ), DIMENSION(jpj) :: &
214	phtra_adv
215
216	REAL(wp), INTENT ( inout ), DIMENSION(jpi,jpj,jpk) :: &
217	tra, traa
218
219	REAL(wp) :: &
220	za, zbtr, e1, e2, c, dir, fu, fc, fd, & ! temporary scalars
221	coef2, coef3, fho, mask, dx
222
223	REAL(wp), DIMENSION(jpi,jpj) :: &
224	zee
225
226	REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
227	zmask, zlap, dwst, lim
228
229
230
231	!----------------------------------------------------------------------
232	! 0. Initialization (should ot be needed on the whole array ???)
233	!----------------------------------------------------------------------
234
235	zmask = 0.0
236	zlap = 0.0
237	dwst = 0.0
238	lim = 0.0
239
240	!----------------------------------------------------------------------
241	! I. Part 1 : x-direction
242	!----------------------------------------------------------------------
243
244	! ===============
245	DO jk = 1, jpkm1 ! Horizontal slab
246	! ! ===============
247	! Initialization of metric arrays (for z- or s-coordinates)
248	! ---------------------------------------------------------
249	DO jj = 1, jpjm1
250	DO ji = 1, fs_jpim1 ! vector opt.
251	#if defined key_zco
252	! z-coordinates, no vertical scale factors
253	zee(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
254	#else
255	! vertical scale factor are used
256	zee(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
257	#endif
258	END DO
259	END DO
260
261	! Laplacian of tracers (at before time step)
262	! ------------------------------------------
263	!--- First derivative (gradient)
264	DO jj = 1, jpjm1
265	DO ji = 1, fs_jpim1 ! vector opt.
266	zmask(ji,jj,jk) = zee(ji,jj) * ( tra(ji+1,jj ,jk) - tra(ji,jj,jk) )
267	END DO
268	END DO
269	DO jj = 2, jpjm1
270	DO ji = fs_2, fs_jpim1 ! vector opt.
271	#if defined key_zco
272	zee(ji,jj) = e1t(ji,jj) / e2t(ji,jj)
273	#else
274	zee(ji,jj) = e1t(ji,jj) / (e2t(ji,jj) * fse3t(ji,jj,jk))
275	#endif
276	zlap(ji,jj,jk) = zee(ji,jj) * ( zmask(ji,jj,jk) - zmask(ji-1,jj,jk) )
277	END DO
278	END DO
279	!--- Function lim=FU+SL*(FC-FU) used by the limiter
280	!--- Computation of the ustream and downstream lim at the T-points
281	DO jj = 2, jpjm1
282	DO ji = 2, fs_jpim1 ! vector opt.
283	! Upstream in the x-direction for the tracer
284	zmask(ji,jj,jk)=tra(ji-1,jj,jk)+sl(ji,jj,jk)*(tra(ji,jj,jk)-tra(ji-1,jj,jk))
285	! Downstream in the x-direction for the tracer
286	dwst (ji,jj,jk)=tra(ji+1,jj,jk)+sl(ji,jj,jk)*(tra(ji,jj,jk)-tra(ji+1,jj,jk))
287	ENDDO
288	ENDDO
289	END DO
290	!--- Lateral boundary conditions on the laplacian (unchanged sgn)
291	CALL lbc_lnk( zlap(:,:,:), 'T', 1. )
292	!--- Lateral boundary conditions for the lim function
293	CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) ; CALL lbc_lnk( dwst(:,:,:), 'T', 1. )
294	! ===============
295	DO jk = 1, jpkm1 ! Horizontal slab
296	! ! ===============
297	! --- lim at the U-points in function of the direction of the flow
298	! ----------------------------------------------------------------
299	DO jj = 1, jpjm1
300	DO ji = 1, fs_jpim1 ! vector opt.
301	dir = 0.5 + sign(0.5,pun(ji,jj,jk)) ! if pun>0 : diru = 1 otherwise diru = 0
302	lim(ji,jj,jk)=dirzmask(ji,jj,jk)+(1-dir)dwst(ji+1,jj,jk)
303	! Mask at the T-points in the x-direction (mask=0 or mask=1)
304	zmask(ji,jj,jk)=tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2
305	END DO
306	END DO
307	END DO
308	!--- Lateral boundary conditions for the mask
309	CALL lbc_lnk( zmask(:,:,:), 'T', 1. )
310
311	! Horizontal advective fluxes
312	! ---------------------------
313	! ===============
314	DO jk = 1, jpkm1 ! Horizontal slab
315	! ===============
316	dt = z2 * rdttra(jk)
317	!--- tracer flux at u and v-points
318	DO jj = 1, jpjm1
319	DO ji = 1, fs_jpim1 ! vector opt.
320	#if defined key_zco
321	e2 = e2u(ji,jj)
322	#else
323	e2 = e2u(ji,jj) * fse3u(ji,jj,jk)
324	#endif
325	dir = 0.5 + sign(0.5,pun(ji,jj,jk)) ! if pun>0 : diru = 1 otherwise diru = 0
326
327	dx = dir * e1t(ji,jj) + (1-dir)* e1t(ji+1,jj)
328	c = ABS(pun(ji,jj,jk))*dt/dx ! (0<cx<1 : Courant number on x-direction)
329
330	fu = lim(ji,jj,jk) ! FU + sl(FC-FU) in the x-direction for T
331	fc = dirtra(ji ,jj,jk)+(1-dir)tra(ji+1,jj,jk) ! FC in the x-direction for T
332	fd = dirtra(ji+1,jj,jk)+(1-dir)tra(ji ,jj,jk) ! FD in the x-direction for T
333
334	!--- QUICKEST scheme
335	! Temperature on the x-direction
336	coef1 = 0.5*(fc+fd)
337	coef2 = 0.5c(fd-fc)
338	coef3 = ((1.-(cc))/6.)(dirzlap(ji,jj,jk) + (1-dir)zlap(ji+1,jj,jk) )
339	fho = coef1-coef2-coef3
340	fho = bound(fu,fd,fc,fho)
341	!--- If the second ustream point is a land point
342	!--- the flux is computed by the 1st order UPWIND scheme
343	mask=dirzmask(ji,jj,jk)+(1-dir)zmask(ji+1,jj,jk)
344	fho = maskfho + (1-mask)fc
345	dwst(ji,jj,jk)=e2pun(ji,jj,jk)fho
346	END DO
347	END DO
348
349	!--- Tracer flux divergence at t-point added to the general trend
350	DO jj = 2, jpjm1
351	DO ji = fs_2, fs_jpim1 ! vector opt.
352	!--- horizontal advective trends
353	#if defined key_zco
354	zbtr = zbtr2(ji,jj)
355	#else
356	zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
357	#endif
358	za = - zbtr * ( dwst(ji,jj,jk) - dwst(ji-1,jj ,jk) )
359	!--- add it to the general tracer trends
360	traa(ji,jj,jk) = traa(ji,jj,jk) + za
361	END DO
362	END DO
363	! ! ===============
364	END DO ! End of slab
365	! ! ===============
366	!----------------------------------------------------------------------
367	! I. Part 2 : y-direction
368	!----------------------------------------------------------------------
369	! ==============
370	DO jk = 1, jpkm1 ! Horizontal slab
371	! ! ===============
372	! Initialization of metric arrays (for z- or s-coordinates)
373	! ---------------------------------------------------------
374	DO jj = 1, jpjm1
375	DO ji = 1, fs_jpim1 ! vector opt.
376	#if defined key_zco
377	! z-coordinates, no vertical scale factors
378	zee(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
379	#else
380	! s-coordinates, vertical scale factor are used
381	zee(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
382	#endif
383	END DO
384	END DO
385
386	! Laplacian of tracers (at before time step)
387	! ------------------------------------------
388	!--- First derivative (gradient)
389	DO jj = 1, jpjm1
390	DO ji = 1, fs_jpim1 ! vector opt.
391	zmask(ji,jj,jk) = zee(ji,jj) * ( tra(ji ,jj+1,jk) - tra(ji,jj,jk) )
392	END DO
393	END DO
394	!--- Second derivative (divergence)
395	DO jj = 2, jpjm1
396	DO ji = fs_2, fs_jpim1 ! vector opt.
397	#if defined key_zco
398	zee(ji,jj) = e2t(ji,jj) / e1t(ji,jj)
399	#else
400	zee(ji,jj) = e2t(ji,jj) / (e1t(ji,jj) * fse3t(ji,jj,jk))
401	#endif
402	zlap(ji,jj,jk) = zee(ji,jj) * ( zmask(ji,jj,jk) - zmask(ji,jj-1,jk) )
403	END DO
404	END DO
405	!--- Function lim=FU+SL*(FC-FU) used by the limiter
406	!--- Computation of the ustream and downstream lim at the T-points
407	DO jj = 2, jpjm1
408	DO ji = 2, fs_jpim1 ! vector opt.
409	! Upstream in the y-direction for the tracer
410	zmask(ji,jj,jk)=tra(ji,jj-1,jk)+sl(ji,jj,jk)*(tra(ji,jj,jk)-tra(ji,jj-1,jk))
411	! Downstream in the y-direction for the tracer
412	dwst (ji,jj,jk)=tra(ji,jj+1,jk)+sl(ji,jj,jk)*(tra(ji,jj,jk)-tra(ji,jj+1,jk))
413	ENDDO
414	ENDDO
415	END DO
416	!--- Lateral boundary conditions on the laplacian (unchanged sgn)
417	CALL lbc_lnk( zlap(:,:,:), 'T', 1. )
418	!--- Lateral boundary conditions for the lim function
419	CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) ; CALL lbc_lnk( dwst(:,:,:), 'T', 1. )
420
421	DO jk = 1, jpkm1 ! Horizontal slab
422	! ! ===============
423	! --- lim at the V-points in function of the direction of the flow
424	! ----------------------------------------------------------------
425	DO jj = 1, jpjm1
426	DO ji = 1, fs_jpim1 ! vector opt.
427	dir = 0.5 + sign(0.5,pvn(ji,jj,jk)) ! if pvn>0 : dirv = 1 otherwise dirv = 0
428	lim(ji,jj,jk)=dirzmask(ji,jj,jk)+(1-dir)dwst(ji,jj+1,jk)
429	! Mask at the T-points in the y-direction (mask=0 or mask=1)
430	zmask(ji,jj,jk)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2
431	END DO
432	END DO
433	END DO
434	!--- Lateral boundary conditions for the mask
435	CALL lbc_lnk( zmask(:,:,:), 'T', 1. )
436
437	! Horizontal advective fluxes
438	! -------------------------------
439	! ===============
440	DO jk = 1, jpkm1 ! Horizontal slab
441	! ===============
442	dt = z2 * rdttra(jk)
443	!--- tracer flux at u and v-points
444	DO jj = 1, jpjm1
445	DO ji = 1, fs_jpim1 ! vector opt.
446	#if defined key_zco
447	e1 = e1v(ji,jj)
448	#else
449	e1 = e1v(ji,jj) * fse3v(ji,jj,jk)
450	#endif
451	dir = 0.5 + sign(0.5,pvn(ji,jj,jk)) ! if pvn>0 : dirv = 1 otherwise dirv = 0
452
453	dx = dir * e2t(ji,jj) + (1-dir)* e2t(ji,jj+1)
454	c = ABS(pvn(ji,jj,jk))*dt/dx ! (0<cy<1 : Courant number on y-direction)
455
456	fu = lim(ji,jj,jk) ! FU + sl(FC-FU) in the y-direction for T
457	fc = dirtra(ji,jj ,jk)+(1-dir)tra(ji,jj+1,jk) ! FC in the y-direction for T
458	fd = dirtra(ji,jj+1,jk)+(1-dir)tra(ji,jj ,jk) ! FD in the y-direction for T
459
460	!--- QUICKEST scheme
461	! Temperature on the y-direction
462	coef1 = 0.5*(fc+fd)
463	coef2 = 0.5c(fd-fc)
464	coef3 = ((1.-(cc))/6.)(dirzlap(ji,jj,jk) + (1-dir)zlap(ji,jj+1,jk) )
465	fho = coef1-coef2-coef3
466	fho = bound(fu,fd,fc,fho)
467	!--- If the second ustream point is a land point
468	!--- the flux is computed by the 1st order UPWIND scheme
469	mask=dirzmask(ji,jj,jk)+(1-dir)zmask(ji,jj+1,jk)
470	fho = maskfho + (1-mask)fc
471	dwst(ji,jj,jk)=e1pvn(ji,jj,jk)fho
472	END DO
473	END DO
474
475	!--- Tracer flux divergence at t-point added to the general trend
476	DO jj = 2, jpjm1
477	DO ji = fs_2, fs_jpim1 ! vector opt.
478	!--- horizontal advective trends
479	#if defined key_zco
480	zbtr = zbtr2(ji,jj)
481	#else
482	zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
483	#endif
484	za = - zbtr * ( dwst(ji,jj,jk) - dwst(ji ,jj-1,jk) )
485	!--- add it to the general tracer trends
486	traa(ji,jj,jk) = traa(ji,jj,jk) + za
487	END DO
488	END DO
489	! ! ===============
490	END DO ! End of slab
491	! ! ===============
492
493	! "zonal" mean advective heat and salt transport
494	IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
495	#if defined key_zco
496	DO jk = 1, jpkm1
497	DO jj = 2, jpjm1
498	DO ji = fs_2, fs_jpim1 ! vector opt.
499	dwst(ji,jj,jk) = dwst(ji,jj,jk) * fse3v(ji,jj,jk)
500	END DO
501	END DO
502	END DO
503	phtra_adv(:) = ptr_vj( dwst(:,:,:) )
504	#else
505	phtra_adv(:) = ptr_vj( dwst(:,:,:) )
506	# endif
507	ENDIF
508
509	END SUBROUTINE tra_adv_qck_hor
510
511	SUBROUTINE tra_adv_qck_ver ( pwn, tra , traa, z2 )
512	!!----------------------------------------------------------------------
513	!!
514	!!----------------------------------------------------------------------
515	!! * Arguments
516
517	REAL(wp), INTENT ( in ) :: z2
518	REAL(wp), INTENT ( in ), DIMENSION(jpi,jpj,jpk) :: &
519	pwn
520	REAL(wp), INTENT ( inout ), DIMENSION(jpi,jpj,jpk) :: &
521	tra, traa
522
523	REAL(wp) :: &
524	za, ze3tr, dt, dir, fc, fd ! temporary scalars
525
526	! Vertical advection
527	! ------------------
528
529	! 1. Vertical advective fluxes
530	! ----------------------------
531
532	!Bottom value : flux set to zero
533	sl(:,:,jpk) = 0.e0
534
535	! Surface value
536	IF( lk_dynspg_rl .OR. lk_vvl ) THEN
537	! rigid lid : flux set to zero
538	sl(:,:, 1 ) = 0.e0
539	ELSE
540	! free surface-constant volume
541	sl(:,:, 1 ) = pwn(:,:,1) * tra(:,:,1)
542	ENDIF
543
544	! Second order centered tracer flux at w-point
545
546	DO jk = 2, jpkm1
547	dt = z2 * rdttra(jk)
548	DO jj = 2, jpjm1
549	DO ji = fs_2, fs_jpim1 ! vector opt.
550	dir = 0.5 + sign(0.5,pwn(ji,jj,jk)) ! if pwn>0 : dirw = 1 otherwise dirw = 0
551	fc = dirtra(ji,jj,jk )fse3t(ji,jj,jk-1)+(1-dir)tra(ji,jj,jk-1)fse3t(ji,jj,jk ) ! FC in the z-direction for T
552	fd = dirtra(ji,jj,jk-1)fse3t(ji,jj,jk )+(1-dir)tra(ji,jj,jk )fse3t(ji,jj,jk-1) ! FD in the z-direction for T
553	!--- Second order centered scheme
554	sl(ji,jj,jk)=pwn(ji,jj,jk)*(fc+fd)/(fse3t(ji,jj,jk-1)+fse3t(ji,jj,jk))
555	END DO
556	END DO
557	END DO
558
559	! 2. Tracer flux divergence at t-point added to the general trend
560	! ---------------------------------------------------------------
561
562	DO jk = 1, jpkm1
563	DO jj = 2, jpjm1
564	DO ji = fs_2, fs_jpim1 ! vector opt.
565	ze3tr = 1. / fse3t(ji,jj,jk)
566	! vertical advective trends
567	za = - ze3tr * ( sl(ji,jj,jk) - sl(ji,jj,jk+1) )
568	! add it to the general tracer trends
569	traa(ji,jj,jk) = traa(ji,jj,jk) + za
570	END DO
571	END DO
572	END DO
573
574	END SUBROUTINE tra_adv_qck_ver
575
576	REAL FUNCTION bound(fu,fd,fc,fho)
577	real :: fu,fd,fc,fho,fref1,fref2
578	fref1 = fu
579	fref2 = MAX(MIN(fc,fd),MIN(MAX(fc,fd),fref1))
580	bound = MAX(MIN(fho,fc),MIN(MAX(fho,fc),fref2))
581	END FUNCTION
582
583	#else
584	!!----------------------------------------------------------------------
585	!! 'key_mpp_omp' : quickest advection (k- and j-slabs)
586	!!----------------------------------------------------------------------
587	SUBROUTINE tra_adv_qck( kt, pun, pvn, pwn )
588	!!----------------------------------------------------------------------
589	!! * Arguments
590	INTEGER, INTENT( in ) :: kt ! ocean time-step index
591	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! effective ocean velocity, u_component
592	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! effective ocean velocity, v_component
593	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! effective ocean velocity, w_component
594	!!----------------------------------------------------------------------
595	IF(lwp) WRITE(numout,*)
596	IF(lwp) WRITE(numout,*) 'tra_adv_ qck:3st order quickest advection scheme'
597	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ Vector optimization case'
598	IF(lwp) WRITE(numout,*) 'WITH AUTOTASKING =>this routine doesn t exist for the moment'
599	IF(lwp) WRITE(numout,*) ' EMPTY ROUTINE!!!!!!'
600
601	END SUBROUTINE tra_adv_qck
602
603	#endif
604
605	!!======================================================================
606	END MODULE traadv_qck

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

Download in other formats: