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 NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_qck.F90 @ 10880

Last change on this file since 10880 was 10880, checked in by davestorkey, 5 years ago

branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps:

Move time indices from dom_oce.F90 to step.F90.
Implement time dimension for passive tracers with independent set of time indices in trcstp.F90.
Update all the traadv and trcadv modules.

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

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

Download in other formats: