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

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

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

Download in other formats: