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/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 9094

Last change on this file since 9094 was 9094, checked in by cetlod, 6 years ago
Use of lbclnk_multi in subdir LDF & TRA
Property svn:keywords set to `Id`
File size: 23.4 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 trc_oce ! share passive tracers/Ocean variables
20	USE trd_oce ! trends: ocean variables
21	USE trdtra ! trends manager: tracers
22	USE diaptr ! poleward transport diagnostics
23	!
24	USE lib_mpp ! distribued memory computing
25	USE lbclnk ! ocean lateral boundary condition (or mpp link)
26	USE in_out_manager ! I/O manager
27	USE timing ! Timing
28	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
29
30	IMPLICIT NONE
31	PRIVATE
32
33	PUBLIC tra_adv_qck ! routine called by step.F90
34
35	REAL(wp) :: r1_6 = 1./ 6. ! 1/6 ratio
36
37	LOGICAL :: l_trd ! flag to compute trends
38	LOGICAL :: l_ptr ! flag to compute poleward transport
39
40
41	!! * Substitutions
42	# include "vectopt_loop_substitute.h90"
43	!!----------------------------------------------------------------------
44	!! NEMO/OPA 4.0 , NEMO Consortium (2017)
45	!! $Id$
46	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
47	!!----------------------------------------------------------------------
48	CONTAINS
49
50	SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn, &
51	& ptb, ptn, pta, kjpt )
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 : u(i)>0
60	!! \|--FU--\|--FC--\|--FD--\|------\|
61	!! i-1 i i+1 i+2
62	!!
63	!! For a negative velocity u : u(i)<0
64	!! \|------\|--FD--\|--FC--\|--FU--\|
65	!! i-1 i i+1 i+2
66	!! where FU is the second upwind point
67	!! FD is the first douwning point
68	!! FC is the central point (or the first upwind point)
69	!!
70	!! Flux(i) = u(i) * { 0.5(FC+FD) -0.5C(i)(FD-FC) -((1-C(i))/6)(FU+FD-2FC) }
71	!! with C(i)=\|u(i)\|dx(i)/dt (=Courant number)
72	!!
73	!! dt = 2*rdtra and the scalar values are tb and sb
74	!!
75	!! On the vertical, the simple centered scheme used ptn
76	!!
77	!! The fluxes are bounded by the ULTIMATE limiter to
78	!! guarantee the monotonicity of the solution and to
79	!! prevent the appearance of spurious numerical oscillations
80	!!
81	!! ** Action : - update pta with the now advective tracer trends
82	!! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
83	!! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
84	!!
85	!! ** Reference : Leonard (1979, 1991)
86	!!----------------------------------------------------------------------
87	INTEGER , INTENT(in ) :: kt ! ocean time-step index
88	INTEGER , INTENT(in ) :: kit000 ! first 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) , INTENT(in ) :: p2dt ! tracer time-step
92	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
93	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
94	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
95	!!----------------------------------------------------------------------
96	!
97	IF( ln_timing ) CALL timing_start('tra_adv_qck')
98	!
99	IF( kt == kit000 ) THEN
100	IF(lwp) WRITE(numout,*)
101	IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
102	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
103	IF(lwp) WRITE(numout,*)
104	ENDIF
105	!
106	l_trd = .FALSE.
107	l_ptr = .FALSE.
108	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
109	IF( cdtype == 'TRA' .AND. ln_diaptr ) l_ptr = .TRUE.
110	!
111	!
112	! ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
113	CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt )
114	CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt )
115
116	! ! vertical fluxes are computed with the 2nd order centered scheme
117	CALL tra_adv_cen2_k( kt, cdtype, pwn, ptn, pta, kjpt )
118	!
119	IF( ln_timing ) CALL timing_stop('tra_adv_qck')
120	!
121	END SUBROUTINE tra_adv_qck
122
123
124	SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun, &
125	& ptb, ptn, pta, kjpt )
126	!!----------------------------------------------------------------------
127	!!
128	!!----------------------------------------------------------------------
129	INTEGER , INTENT(in ) :: kt ! ocean time-step index
130	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
131	INTEGER , INTENT(in ) :: kjpt ! number of tracers
132	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
133	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun ! i-velocity components
134	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
135	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
136	!!
137	INTEGER :: ji, jj, jk, jn ! dummy loop indices
138	REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk ! local scalars
139	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zfu, zfc, zfd
140	!----------------------------------------------------------------------
141	!
142	! ! ===========
143	DO jn = 1, kjpt ! tracer loop
144	! ! ===========
145	zfu(:,:,:) = 0._wp ; zfc(:,:,:) = 0._wp
146	zfd(:,:,:) = 0._wp ; zwx(:,:,:) = 0._wp
147	!
148	!!gm why not using a SHIFT instruction...
149	DO jk = 1, jpkm1 !--- Computation of the ustream and downstream value of the tracer and the mask
150	DO jj = 2, jpjm1
151	DO ji = fs_2, fs_jpim1 ! vector opt.
152	zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn) ! Upstream in the x-direction for the tracer
153	zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn) ! Downstream in the x-direction for the tracer
154	END DO
155	END DO
156	END DO
157	CALL lbc_lnk_multi( zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. ) ! Lateral boundary conditions
158
159	!
160	! Horizontal advective fluxes
161	! ---------------------------
162	DO jk = 1, jpkm1
163	DO jj = 2, jpjm1
164	DO ji = fs_2, fs_jpim1 ! vector opt.
165	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
166	zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T
167	END DO
168	END DO
169	END DO
170	!
171	DO jk = 1, jpkm1
172	DO jj = 2, jpjm1
173	DO ji = fs_2, fs_jpim1 ! vector opt.
174	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
175	zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk)
176	zwx(ji,jj,jk) = ABS( pun(ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)
177	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
178	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
179	END DO
180	END DO
181	END DO
182	!--- Lateral boundary conditions
183	CALL lbc_lnk_multi( zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1., zwx(:,:,:), 'T', 1. )
184
185	!--- QUICKEST scheme
186	CALL quickest( zfu, zfd, zfc, zwx )
187	!
188	! Mask at the T-points in the x-direction (mask=0 or mask=1)
189	DO jk = 1, jpkm1
190	DO jj = 2, jpjm1
191	DO ji = fs_2, fs_jpim1 ! vector opt.
192	zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
193	END DO
194	END DO
195	END DO
196	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ! Lateral boundary conditions
197
198	!
199	! Tracer flux on the x-direction
200	DO jk = 1, jpkm1
201	!
202	DO jj = 2, jpjm1
203	DO ji = fs_2, fs_jpim1 ! vector opt.
204	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
205	!--- If the second ustream point is a land point
206	!--- the flux is computed by the 1st order UPWIND scheme
207	zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
208	zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
209	zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
210	END DO
211	END DO
212	END DO
213	!
214	CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
215	!
216	! Computation of the trend
217	DO jk = 1, jpkm1
218	DO jj = 2, jpjm1
219	DO ji = fs_2, fs_jpim1 ! vector opt.
220	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
221	! horizontal advective trends
222	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
223	!--- add it to the general tracer trends
224	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
225	END DO
226	END DO
227	END DO
228	! ! trend diagnostics
229	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
230	!
231	END DO
232	!
233	END SUBROUTINE tra_adv_qck_i
234
235
236	SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn, &
237	& ptb, ptn, pta, kjpt )
238	!!----------------------------------------------------------------------
239	!!
240	!!----------------------------------------------------------------------
241	INTEGER , INTENT(in ) :: kt ! ocean time-step index
242	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
243	INTEGER , INTENT(in ) :: kjpt ! number of tracers
244	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
245	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pvn ! j-velocity components
246	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
247	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
248	!!
249	INTEGER :: ji, jj, jk, jn ! dummy loop indices
250	REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk ! local scalars
251	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwy, zfu, zfc, zfd ! 3D workspace
252	!----------------------------------------------------------------------
253	!
254	! ! ===========
255	DO jn = 1, kjpt ! tracer loop
256	! ! ===========
257	zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0
258	zfd(:,:,:) = 0.0 ; zwy(:,:,:) = 0.0
259	!
260	DO jk = 1, jpkm1
261	!
262	!--- Computation of the ustream and downstream value of the tracer and the mask
263	DO jj = 2, jpjm1
264	DO ji = fs_2, fs_jpim1 ! vector opt.
265	! Upstream in the x-direction for the tracer
266	zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
267	! Downstream in the x-direction for the tracer
268	zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
269	END DO
270	END DO
271	END DO
272	CALL lbc_lnk_multi( zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. ) ! Lateral boundary conditions
273
274
275	!
276	! Horizontal advective fluxes
277	! ---------------------------
278	!
279	DO jk = 1, jpkm1
280	DO jj = 2, jpjm1
281	DO ji = fs_2, fs_jpim1 ! vector opt.
282	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
283	zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T
284	END DO
285	END DO
286	END DO
287	!
288	DO jk = 1, jpkm1
289	DO jj = 2, jpjm1
290	DO ji = fs_2, fs_jpim1 ! vector opt.
291	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
292	zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk)
293	zwy(ji,jj,jk) = ABS( pvn(ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)
294	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
295	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
296	END DO
297	END DO
298	END DO
299
300	!--- Lateral boundary conditions
301	CALL lbc_lnk_multi( zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1., zwy(:,:,:), 'T', 1. )
302
303	!--- QUICKEST scheme
304	CALL quickest( zfu, zfd, zfc, zwy )
305	!
306	! Mask at the T-points in the x-direction (mask=0 or mask=1)
307	DO jk = 1, jpkm1
308	DO jj = 2, jpjm1
309	DO ji = fs_2, fs_jpim1 ! vector opt.
310	zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
311	END DO
312	END DO
313	END DO
314	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) !--- Lateral boundary conditions
315	!
316	! Tracer flux on the x-direction
317	DO jk = 1, jpkm1
318	!
319	DO jj = 2, jpjm1
320	DO ji = fs_2, fs_jpim1 ! vector opt.
321	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
322	!--- If the second ustream point is a land point
323	!--- the flux is computed by the 1st order UPWIND scheme
324	zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
325	zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
326	zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
327	END DO
328	END DO
329	END DO
330	!
331	CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
332	!
333	! Computation of the trend
334	DO jk = 1, jpkm1
335	DO jj = 2, jpjm1
336	DO ji = fs_2, fs_jpim1 ! vector opt.
337	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
338	! horizontal advective trends
339	ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
340	!--- add it to the general tracer trends
341	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
342	END DO
343	END DO
344	END DO
345	! ! trend diagnostics
346	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
347	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
348	IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
349	!
350	END DO
351	!
352	END SUBROUTINE tra_adv_qck_j
353
354
355	SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn, &
356	& ptn, pta, kjpt )
357	!!----------------------------------------------------------------------
358	!!
359	!!----------------------------------------------------------------------
360	INTEGER , INTENT(in ) :: kt ! ocean time-step index
361	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
362	INTEGER , INTENT(in ) :: kjpt ! number of tracers
363	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pwn ! vertical velocity
364	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptn ! before and now tracer fields
365	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
366	!
367	INTEGER :: ji, jj, jk, jn ! dummy loop indices
368	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwz ! 3D workspace
369	!!----------------------------------------------------------------------
370	!
371	zwz(:,:, 1 ) = 0._wp ! surface & bottom values set to zero for all tracers
372	zwz(:,:,jpk) = 0._wp
373	!
374	! ! ===========
375	DO jn = 1, kjpt ! tracer loop
376	! ! ===========
377	!
378	DO jk = 2, jpkm1 !* Interior point (w-masked 2nd order centered flux)
379	DO jj = 2, jpjm1
380	DO ji = fs_2, fs_jpim1 ! vector opt.
381	zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
382	END DO
383	END DO
384	END DO
385	IF( ln_linssh ) THEN !* top value (only in linear free surf. as zwz is multiplied by wmask)
386	IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean)
387	DO jj = 1, jpj
388	DO ji = 1, jpi
389	zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn) ! linear free surface
390	END DO
391	END DO
392	ELSE ! no ocean cavities (only ocean surface)
393	zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
394	ENDIF
395	ENDIF
396	!
397	DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==!
398	DO jj = 2, jpjm1
399	DO ji = fs_2, fs_jpim1 ! vector opt.
400	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) &
401	& * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
402	END DO
403	END DO
404	END DO
405	! ! Send trends for diagnostic
406	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
407	!
408	END DO
409	!
410	END SUBROUTINE tra_adv_cen2_k
411
412
413	SUBROUTINE quickest( pfu, pfd, pfc, puc )
414	!!----------------------------------------------------------------------
415	!!
416	!! ** Purpose : Computation of advective flux with Quickest scheme
417	!!
418	!! ** Method :
419	!!----------------------------------------------------------------------
420	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfu ! second upwind point
421	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfd ! first douwning point
422	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfc ! the central point (or the first upwind point)
423	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: puc ! input as Courant number ; output as flux
424	!!
425	INTEGER :: ji, jj, jk ! dummy loop indices
426	REAL(wp) :: zcoef1, zcoef2, zcoef3 ! local scalars
427	REAL(wp) :: zc, zcurv, zfho ! - -
428	!----------------------------------------------------------------------
429	!
430	IF( ln_timing ) CALL timing_start('quickest')
431	!
432	DO jk = 1, jpkm1
433	DO jj = 1, jpj
434	DO ji = 1, jpi
435	zc = puc(ji,jj,jk) ! Courant number
436	zcurv = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
437	zcoef1 = 0.5 * ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
438	zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
439	zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
440	zfho = zcoef1 - zcoef2 - zcoef3 ! phi_f QUICKEST
441	!
442	zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
443	zcoef2 = ABS( zcoef1 )
444	zcoef3 = ABS( zcurv )
445	IF( zcoef3 >= zcoef2 ) THEN
446	zfho = pfc(ji,jj,jk)
447	ELSE
448	zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) ) ! phi_REF
449	IF( zcoef1 >= 0. ) THEN
450	zfho = MAX( pfc(ji,jj,jk), zfho )
451	zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) )
452	ELSE
453	zfho = MIN( pfc(ji,jj,jk), zfho )
454	zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) )
455	ENDIF
456	ENDIF
457	puc(ji,jj,jk) = zfho
458	END DO
459	END DO
460	END DO
461	!
462	IF( ln_timing ) CALL timing_stop('quickest')
463	!
464	END SUBROUTINE quickest
465
466	!!======================================================================
467	END MODULE traadv_qck

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

Download in other formats: