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 @ 789

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

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

Download in other formats: