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

Context Navigation

source: branches/dev_004_VVL/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 1388

Last change on this file since 1388 was 1231, checked in by ctlod, 15 years ago
replace the traadv_qck.F90 module with a new one entirely rewritten, see ticket: #251
File size: 23.5 KB

Line
1	MODULE traadv_qck
2	!!==============================================================================
3	!! * MODULE traadv_qck *
4	!! Ocean active tracers: horizontal & vertical advective trend
5	!!==============================================================================
6	!! History : 9.0 ! 2008-07 (G. Reffray) Original code
7	!!----------------------------------------------------------------------
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_zon :
14	!! tra_adv_qck_mer :
15	!! tra_adv_cen2_ver : 2nd centered scheme for the vertical advection
16	!!----------------------------------------------------------------------
17	!! * Modules used
18	USE oce ! ocean dynamics and active tracers
19	USE dom_oce ! ocean space and time domain
20	USE trdmod ! ocean active tracers trends
21	USE trdmod_oce ! ocean variables trends
22	USE trabbl ! advective term in the BBL
23	USE lib_mpp ! distribued memory computing
24	USE lbclnk ! ocean lateral boundary condition (or mpp link)
25	USE dynspg_oce ! surface pressure gradient variables
26	USE in_out_manager ! I/O manager
27	USE diaptr ! poleward transport diagnostics
28	USE prtctl ! Print control
29
30	IMPLICIT NONE
31	PRIVATE
32
33	!! * Accessibility
34	PUBLIC tra_adv_qck ! routine called by step.F90
35
36	!! * Module variables
37	REAL(wp), DIMENSION(jpi,jpj), SAVE :: btr2
38	REAL(wp) , SAVE :: cst
39	!!----------------------------------------------------------------------
40	!! * Substitutions
41	# include "domzgr_substitute.h90"
42	# include "vectopt_loop_substitute.h90"
43	!!----------------------------------------------------------------------
44	!! OPA 9.0 , LOCEAN-IPSL (2008)
45	!! $Id$
46	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
47	!!----------------------------------------------------------------------
48
49	CONTAINS
50
51	SUBROUTINE tra_adv_qck( kt, pun, pvn, pwn )
52	!!----------------------------------------------------------------------
53	!! * ROUTINE tra_adv_qck *
54	!!
55	!! ** Purpose : Compute the now trend due to the advection of tracers
56	!! and add it to the general trend of passive tracer equations.
57	!!
58	!! ** Method : The advection is evaluated by a third order scheme
59	!! For a positive velocity u :
60	!!
61	!!
62	!! i-1 i i+1 i+2
63	!!
64	!! \|--FU--\|--FC--\|--FD--\|------\|
65	!! u(i)>0
66	!!
67	!! For a negative velocity u :
68	!!
69	!! \|------\|--FD--\|--FC--\|--FU--\|
70	!! u(i)<0
71	!!
72	!! FU is the second upwind point
73	!! FD is the first douwning point
74	!! FC is the central point (or the first upwind point)
75	!!
76	!! Flux(i) = u(i) * {0.5(FC+FD) -0.5C(i)(FD-FC) -((1-C(i))/6)(FU+FD-2FC)}
77	!! with C(i)=\|u(i)\|dx(i)/dt (Courant number)
78	!!
79	!! dt = 2*rdtra and the scalar values are tb and sb
80	!!
81	!! On the vertical, the simple centered scheme used tn and sn
82	!!
83	!! The fluxes are bounded by the ULTIMATE limiter
84	!! to guarantee the monotonicity of the solution and to
85	!! prevent the appearance of spurious numerical oscillations
86	!!
87	!! ** Action : - update (ta,sa) with the now advective tracer trends
88	!! - save the trends ('key_trdtra')
89	!!
90	!! ** Reference : Leonard (1979, 1991)
91	!!----------------------------------------------------------------------
92	USE oce, ONLY : ztrdt => ua ! use ua as workspace
93	USE oce, ONLY : ztrds => va ! use va as workspace
94	!! * Arguments
95	INTEGER , INTENT(in) :: kt ! ocean time-step index
96	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! effective ocean velocity, u_component
97	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! effective ocean velocity, v_component
98	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! effective ocean velocity, w_component
99	!!
100	INTEGER :: ji, jj, jk ! dummy loop indices
101	REAL(wp) :: z_hdivn_x, z_hdivn_y, z_hdivn ! temporary scalars
102	REAL(wp) :: zbtr, z2 ! " "
103	!!----------------------------------------------------------------------
104
105	IF( kt == nit000 ) THEN
106	IF(lwp) WRITE(numout,*)
107	IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme'
108	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
109	IF(lwp) WRITE(numout,*)
110	btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
111	cst = 1. / 6.
112	ENDIF
113
114	IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1.
115	ELSE ; z2 = 2.
116	ENDIF
117
118	! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
119	!---------------------------------------------------------------------------
120
121	CALL tra_adv_qck_zon( kt, pun, tb, tn, ta, ztrdt, z2)
122	CALL tra_adv_qck_zon( kt, pun, sb, sn, sa, ztrds, z2)
123
124	IF( l_trdtra ) CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)
125
126	CALL tra_adv_qck_mer( kt, pvn, tb, tn, ta, ztrdt, pht_adv, z2)
127	CALL tra_adv_qck_mer( kt, pvn, sb, sn, sa, ztrds, pst_adv, z2)
128
129	IF( l_trdtra ) THEN
130	CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)
131	!
132	! Save the horizontal up-to-date ta/sa trends
133	ztrdt(:,:,:) = ta(:,:,:)
134	ztrds(:,:,:) = sa(:,:,:)
135	END IF
136
137	IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' qck had - Ta: ', mask1=tmask, &
138	& tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
139
140	! II. The vertical fluxes are computed with the 2nd order centered scheme
141	!-------------------------------------------------------------------------
142	!
143	CALL tra_adv_cen2_ver( pwn, tn, ta )
144	CALL tra_adv_cen2_ver( pwn, sn, sa )
145	!
146	!Save the vertical advective trends for diagnostic
147	! ----------------------------------------------------
148	IF( l_trdtra ) THEN
149	! Recompute the vertical advection zta & zsa trends computed
150	! at the step 2. above in making the difference between the new
151	! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract
152	! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
153
154	DO jk = 1, jpkm1
155	DO jj = 2, jpjm1
156	DO ji = fs_2, fs_jpim1 ! vector opt.
157	#if defined key_zco
158	zbtr = btr2(ji,jj)
159	z_hdivn_x = e2u(ji,jj)pun(ji,jj,jk) - e2u(ji-1,jj)pun(ji-1,jj,jk)
160	z_hdivn_y = e1v(ji,jj)pvn(ji,jj,jk) - e1v(ji,jj-1)pvn(ji,jj-1,jk)
161	#else
162	zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
163	z_hdivn_x = e2u(ji,jj)fse3u(ji,jj,jk)pun(ji,jj,jk) - e2u(ji-1,jj)fse3u(ji-1,jj,jk)pun(ji-1,jj,jk)
164	z_hdivn_y = e1v(ji,jj)fse3v(ji,jj,jk)pvn(ji,jj,jk) - e1v(ji,jj-1)fse3v(ji,jj-1,jk)pvn(ji,jj-1,jk)
165	#endif
166	z_hdivn = (z_hdivn_x + z_hdivn_y) * zbtr
167	ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn
168	ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn
169	END DO
170	END DO
171	END DO
172	CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)
173	ENDIF
174
175	IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' qck zad - Ta: ', mask1=tmask, &
176	& tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
177	!
178	END SUBROUTINE tra_adv_qck
179
180
181	SUBROUTINE tra_adv_qck_zon ( kt, pun, tra, tran, traa, ztrdtra, z2 )
182	!!----------------------------------------------------------------------
183	!!
184	!!----------------------------------------------------------------------
185	!! * Arguments
186	INTEGER, INTENT(in) :: kt ! ocean time-step index
187	REAL, INTENT(in) :: z2
188	REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pun, tra, tran ! horizontal effective velocity
189	REAL(wp), INTENT(out) , DIMENSION(jpi,jpj,jpk) :: ztrdtra
190	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: traa
191	!
192	INTEGER :: ji, jj, jk
193	REAL(wp) :: za, zbtr, dir, dx, dt ! temporary scalars
194	REAL(wp) :: z_hdivn_x
195	REAL(wp), DIMENSION(jpi,jpj) :: zmask, zupst, zdwst, zc_cfl
196	REAL(wp), DIMENSION(jpi,jpj) :: zfu, zfc, zfd, zfho, zmskl, zsc_e
197	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux
198	!----------------------------------------------------------------------
199	! I. Part 1 : x-direction
200	!----------------------------------------------------------------------
201
202	zfu (:,jpj) = 0.e0 ; zfc (:,jpj) = 0.e0
203	zfd (:,jpj) = 0.e0 ; zc_cfl(:,jpj) = 0.e0
204	zsc_e (:,jpj) = 0.e0 ; zmskl (:,jpj) = 0.e0
205	zfho (:,jpj) = 0.e0
206	! ===============
207	DO jk = 1, jpkm1 ! Horizontal slab
208	! ! ===============
209	!--- Computation of the ustream and downstream value of the tracer and the mask
210	DO jj = 2, jpjm1
211	DO ji = 2, fs_jpim1 ! vector opt.
212	! Upstream in the x-direction for the tracer
213	zupst(ji,jj)=tra(ji-1,jj,jk)
214	! Downstream in the x-direction for the tracer
215	zdwst(ji,jj)=tra(ji+1,jj,jk)
216	! Mask at the T-points in the x-direction (mask=0 or mask=1)
217	zmask(ji,jj)=tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2
218	ENDDO
219	ENDDO
220	!
221	!--- Lateral boundary conditions
222	CALL lbc_lnk( zupst(:,:), 'T', 1. )
223	CALL lbc_lnk( zdwst(:,:), 'T', 1. )
224	CALL lbc_lnk( zmask(:,:), 'T', 1. )
225	!
226	! Horizontal advective fluxes
227	! ---------------------------
228	!
229	dt = z2 * rdttra(jk)
230	!--- tracer flux at u-points
231	DO jj = 1, jpjm1
232	DO ji = 1, jpi
233	#if defined key_zco
234	zsc_e(ji,jj) = e2u(ji,jj)
235	#else
236	zsc_e(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk)
237	#endif
238	dir = 0.5 + sign(0.5,pun(ji,jj,jk)) ! if pun>0 : dir = 1 otherwise dir = 0
239	dx = dir * e1t(ji,jj) + (1-dir)* e1t(ji+1,jj)
240	zc_cfl (ji,jj) = ABS(pun(ji,jj,jk))*dt/dx ! (0<zc_cfl<1 : Courant number on x-direction)
241
242	zfu(ji,jj) = dirzupst(ji ,jj )+(1-dir)zdwst(ji+1,jj ) ! FU in the x-direction for T
243	zfc(ji,jj) = dirtra (ji ,jj,jk)+(1-dir)tra (ji+1,jj,jk) ! FC in the x-direction for T
244	zfd(ji,jj) = dirtra (ji+1,jj,jk)+(1-dir)tra (ji ,jj,jk) ! FD in the x-direction for T
245	zmskl(ji,jj) = dirzmask(ji ,jj) +(1-dir)zmask(ji+1,jj)
246	ENDDO
247	ENDDO
248	!
249	!--- QUICKEST scheme
250	! Tracer flux on the x-direction
251	CALL quickest(zfu,zfd,zfc,zfho,zc_cfl)
252	!--- If the second ustream point is a land point
253	!--- the flux is computed by the 1st order UPWIND scheme
254	zfho(:,:) = zmskl(:,:)zfho(:,:) + (1.-zmskl(:,:))zfc(:,:)
255	!--- Computation of fluxes
256	zflux(:,:,jk) = zsc_e(:,:)pun(:,:,jk)zfho(:,:)
257	!
258	!--- Tracer flux divergence at t-point added to the general trend
259	DO jj = 2, jpjm1
260	DO ji = fs_2, fs_jpim1 ! vector opt.
261	!--- horizontal advective trends
262	#if defined key_zco
263	zbtr = btr2(ji,jj)
264	#else
265	zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
266	#endif
267	za = - zbtr * ( zflux(ji,jj,jk) - zflux(ji-1,jj,jk) )
268	!--- add it to the general tracer trends
269	traa(ji,jj,jk) = traa(ji,jj,jk) + za
270	END DO
271	END DO
272	! ! ===============
273	END DO ! End of slab
274	! ! ===============
275	!
276	! Save the horizontal advective trends for diagnostic
277	! -----------------------------------------------------
278	IF( l_trdtra ) THEN
279	! T/S ZONAL advection trends
280	ztrdtra(:,:,:) = 0.e0
281	!
282	DO jk = 1, jpkm1
283	DO jj = 2, jpjm1
284	DO ji = fs_2, fs_jpim1 ! vector opt.
285	!-- Compute zonal divergence by splitting hdivn (see divcur.F90)
286	! N.B. This computation is not valid along OBCs (if any)
287	#if defined key_zco
288	zbtr = btr2(ji,jj)
289	z_hdivn_x = ( e2u(ji ,jj) * pun(ji ,jj,jk) &
290	& - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr
291	#else
292	zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
293	z_hdivn_x = ( e2u(ji ,jj) * fse3u(ji ,jj,jk) * pun(ji ,jj,jk) &
294	& - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr
295	#endif
296	ztrdtra(ji,jj,jk) = - zbtr * ( zflux(ji,jj,jk) - zflux(ji-1,jj,jk) ) + tran(ji,jj,jk) * z_hdivn_x
297	END DO
298	END DO
299	END DO
300	END IF
301
302	END SUBROUTINE tra_adv_qck_zon
303
304
305	SUBROUTINE tra_adv_qck_mer ( kt, pvn, tra, tran, traa, ztrdtra, trd_adv, z2 )
306	!!----------------------------------------------------------------------
307	!!
308	!!----------------------------------------------------------------------
309	!! * Arguments
310	INTEGER, INTENT(in) :: kt ! ocean time-step index
311	REAL, INTENT(in) :: z2
312	REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pvn, tra, tran ! horizontal effective velocity
313	REAL(wp), INTENT(out) , DIMENSION(jpj) :: trd_adv
314	REAL(wp), INTENT(out) , DIMENSION(jpi,jpj,jpk) :: ztrdtra
315	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: traa
316	!
317	INTEGER :: ji, jj, jk
318	REAL(wp) :: za, zbtr, dir, dx, dt ! temporary scalars
319	REAL(wp) :: z_hdivn_y
320	REAL(wp), DIMENSION(jpi,jpj) :: zmask, zupst, zdwst, zc_cfl
321	REAL(wp), DIMENSION(jpi,jpj) :: zfu, zfc, zfd, zfho, zmskl, zsc_e
322	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux
323	!----------------------------------------------------------------------
324	! II. Part 2 : y-direction
325	!----------------------------------------------------------------------
326
327	zfu (:,jpj) = 0.e0 ; zfc (:,jpj) = 0.e0
328	zfd (:,jpj) = 0.e0 ; zc_cfl(:,jpj) = 0.e0
329	zsc_e (:,jpj) = 0.e0 ; zmskl (:,jpj) = 0.e0
330	zfho (:,jpj) = 0.e0
331
332	! ===============
333	DO jk = 1, jpkm1 ! Horizontal slab
334	! ! ===============
335	!--- Computation of the ustream and downstream value of the tracer and the mask
336	DO jj = 2, jpjm1
337	DO ji = 2, fs_jpim1 ! vector opt.
338	! Upstream in the x-direction for the tracer
339	zupst(ji,jj)=tra(ji,jj-1,jk)
340	! Downstream in the x-direction for the tracer
341	zdwst(ji,jj)=tra(ji,jj+1,jk)
342	! Mask at the T-points in the x-direction (mask=0 or mask=1)
343	zmask(ji,jj)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2
344	ENDDO
345	ENDDO
346	!
347	!--- Lateral boundary conditions
348	CALL lbc_lnk( zupst(:,:), 'T', 1. )
349	CALL lbc_lnk( zdwst(:,:), 'T', 1. )
350	CALL lbc_lnk( zmask(:,:), 'T', 1. )
351	!
352	! Horizontal advective fluxes
353	! ---------------------------
354	!
355	dt = z2 * rdttra(jk)
356	!--- tracer flux at v-points
357	DO jj = 1, jpjm1
358	DO ji = 1, jpi
359	#if defined key_zco
360	zsc_e(ji,jj) = e1v(ji,jj)
361	#else
362	zsc_e(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk)
363	#endif
364	dir = 0.5 + sign(0.5,pvn(ji,jj,jk)) ! if pvn>0 : dir = 1 otherwise dir = 0
365	dx = dir * e2t(ji,jj) + (1-dir)* e2t(ji,jj+1)
366	zc_cfl(ji,jj) = ABS(pvn(ji,jj,jk))*dt/dx ! (0<zc_cfl<1 : Courant number on y-direction)
367
368	zfu(ji,jj) = dirzupst(ji,jj )+(1-dir)zdwst(ji,jj+1 ) ! FU in the y-direction for T
369	zfc(ji,jj) = dirtra (ji,jj ,jk)+(1-dir)tra (ji,jj+1,jk) ! FC in the y-direction for T
370	zfd(ji,jj) = dirtra (ji,jj+1,jk)+(1-dir)tra (ji,jj ,jk) ! FD in the y-direction for T
371	zmskl(ji,jj) = dirzmask(ji,jj )+(1-dir)zmask(ji,jj+1)
372	ENDDO
373	ENDDO
374	!
375	!--- QUICKEST scheme
376	! Tracer flux on the y-direction
377	CALL quickest(zfu,zfd,zfc,zfho,zc_cfl)
378	!--- If the second ustream point is a land point
379	!--- the flux is computed by the 1st order UPWIND scheme
380	zfho(:,:) = zmskl(:,:)zfho(:,:) + (1.-zmskl(:,:))zfc(:,:)
381	!--- Computation of fluxes
382	zflux(:,:,jk) = zsc_e(:,:)pvn(:,:,jk)zfho(:,:)
383	!
384	!--- Tracer flux divergence at t-point added to the general trend
385	DO jj = 2, jpjm1
386	DO ji = fs_2, fs_jpim1 ! vector opt.
387	!--- horizontal advective trends
388	#if defined key_zco
389	zbtr = btr2(ji,jj)
390	#else
391	zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
392	#endif
393	za = - zbtr * ( zflux(ji,jj,jk) - zflux(ji,jj-1,jk) )
394	!--- add it to the general tracer trends
395	traa(ji,jj,jk) = traa(ji,jj,jk) + za
396	END DO
397	END DO
398	! ! ===============
399	END DO ! End of slab
400	! ! ===============
401	!
402	! Save the horizontal advective trends for diagnostic
403	! -----------------------------------------------------
404	IF( l_trdtra ) THEN
405	! T/S MERIDIONAL advection trends
406	DO jk = 1, jpkm1
407	DO jj = 2, jpjm1
408	DO ji = fs_2, fs_jpim1 ! vector opt.
409	!-- Compute merid. divergence by splitting hdivn (see divcur.F90)
410	! N.B. This computation is not valid along OBCs (if any)
411	#if defined key_zco
412	zbtr = btr2(ji,jj)
413	z_hdivn_y = ( e1v(ji,jj ) * pvn(ji,jj ,jk) &
414	& - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr
415	#else
416	zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
417	z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) &
418	& - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr
419	#endif
420	ztrdtra(ji,jj,jk) = - zbtr * ( zflux(ji,jj,jk) - zflux(ji,jj-1,jk) ) + tran(ji,jj,jk) * z_hdivn_y
421	END DO
422	END DO
423	END DO
424	END IF
425
426	! "zonal" mean advective heat and salt transport
427	! ----------------------------------------------
428
429	IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
430	IF( lk_zco ) THEN
431	DO jk = 1, jpkm1
432	DO jj = 2, jpjm1
433	DO ji = fs_2, fs_jpim1 ! vector opt.
434	zflux(ji,jj,jk) = zflux(ji,jj,jk) * fse3v(ji,jj,jk)
435	END DO
436	END DO
437	END DO
438	ENDIF
439	trd_adv(:) = ptr_vj( zflux(:,:,:) )
440	ENDIF
441
442	END SUBROUTINE tra_adv_qck_mer
443
444
445	SUBROUTINE tra_adv_cen2_ver ( pwn, tra , traa )
446	!!----------------------------------------------------------------------
447	!!
448	!!----------------------------------------------------------------------
449	!! * Arguments
450	REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pwn ! horizontal effective velocity
451	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: tra, traa
452	!
453	INTEGER ji, jj, jk
454	REAL(wp) :: za, ze3tr ! temporary scalars
455	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux
456	!
457	! Vertical advection
458	! ------------------
459	!
460	! 1. Vertical advective fluxes
461	! ----------------------------
462	!
463	! Bottom value : flux set to zero
464	zflux(:,:,jpk) = 0.e0
465	!
466	! Surface value
467	IF( lk_dynspg_rl .OR. lk_vvl ) THEN
468	! rigid lid : flux set to zero
469	zflux(:,:, 1 ) = 0.e0
470	ELSE
471	! free surface-constant volume
472	zflux(:,:, 1 ) = pwn(:,:,1) * tra(:,:,1)
473	ENDIF
474	!
475	! Second order centered tracer flux at w-point
476	!
477	DO jk = 2, jpkm1
478	DO jj = 2, jpjm1
479	DO ji = fs_2, fs_jpim1 ! vector opt.
480	zflux(ji,jj,jk) = pwn(ji,jj,jk)0.5( tra(ji,jj,jk-1) + tra(ji,jj,jk) )
481	END DO
482	END DO
483	END DO
484	!
485	! 2. Tracer flux divergence at t-point added to the general trend
486	! ---------------------------------------------------------------
487	!
488	DO jk = 1, jpkm1
489	DO jj = 2, jpjm1
490	DO ji = fs_2, fs_jpim1 ! vector opt.
491	ze3tr = 1. / fse3t(ji,jj,jk)
492	! vertical advective trends
493	za = - ze3tr * ( zflux(ji,jj,jk) - zflux(ji,jj,jk+1) )
494	! add it to the general tracer trends
495	traa(ji,jj,jk) = traa(ji,jj,jk) + za
496	END DO
497	END DO
498	END DO
499	!
500	END SUBROUTINE tra_adv_cen2_ver
501
502
503	SUBROUTINE quickest(fu,fd,fc,fho,fc_cfl)
504	!!----------------------------------------------------------------------
505	!!
506	!!----------------------------------------------------------------------
507	!! * Arguments
508	REAL(wp), INTENT(in) , DIMENSION(jpi,jpj) :: fu, fd, fc, fc_cfl
509	REAL(wp), INTENT(out) , DIMENSION(jpi,jpj) :: fho
510	REAL(wp) , DIMENSION(jpi,jpj) :: zcurv, zcoef1, zcoef2, zcoef3 ! temporary scalars
511	!
512	zcurv (:,:) = fd(:,:) + fu(:,:) - 2.*fc(:,:)
513	zcoef1(:,:) = 0.5*( fc(:,:) + fd(:,:) )
514	zcoef2(:,:) = 0.5fc_cfl(:,:)( fd(:,:) - fc(:,:) )
515	zcoef3(:,:) = ( ( 1. - ( fc_cfl(:,:)fc_cfl(:,:) ) )cst )*zcurv(:,:)
516	fho (:,:) = zcoef1(:,:) - zcoef2(:,:) - zcoef3(:,:) ! phi_f QUICKEST
517	!
518	zcoef1(:,:) = fd(:,:) - fu(:,:) ! DEL
519	zcoef2(:,:) = ABS( zcoef1(:,:) ) ! ABS(DEL)
520	zcoef3(:,:) = ABS( zcurv(:,:) ) ! ABS(CURV)
521	!
522	WHERE ( zcoef3(:,:) >= zcoef2(:,:) )
523	fho(:,:) = fc(:,:)
524	ELSEWHERE
525	zcoef3(:,:) = fu(:,:) + ( ( fc(:,:) - fu(:,:) )/MAX(fc_cfl(:,:),1.e-9) ) ! phi_REF
526	WHERE ( zcoef1(:,:) >= 0.e0 )
527	fho(:,:) = MAX(fc(:,:),fho(:,:))
528	fho(:,:) = MIN(fho(:,:),MIN(zcoef3(:,:),fd(:,:)))
529	ELSEWHERE
530	fho(:,:) = MIN(fc(:,:),fho(:,:))
531	fho(:,:) = MAX(fho(:,:),MAX(zcoef3(:,:),fd(:,:)))
532	ENDWHERE
533	ENDWHERE
534	!
535	END SUBROUTINE quickest
536
537	!!======================================================================
538	END MODULE traadv_qck

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

Download in other formats: