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_tvd.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90 @ 4444

Last change on this file since 4444 was 4444, checked in by trackstand2, 10 years ago
Use of mbkmax in tra_adv_tvd
Property svn:keywords set to `Id`
File size: 22.7 KB

Line
1	MODULE traadv_tvd
2	!!==============================================================================
3	!! * MODULE traadv_tvd *
4	!! Ocean tracers: horizontal & vertical advective trend
5	!!==============================================================================
6	!! History : OPA ! 1995-12 (L. Mortier) Original code
7	!! ! 2000-01 (H. Loukos) adapted to ORCA
8	!! ! 2000-10 (MA Foujols E.Kestenare) include file not routine
9	!! ! 2000-12 (E. Kestenare M. Levy) fix bug in trtrd indexes
10	!! ! 2001-07 (E. Durand G. Madec) adaptation to ORCA config
11	!! 8.5 ! 2002-06 (G. Madec) F90: Free form and module
12	!! NEMO 1.0 ! 2004-01 (A. de Miranda, G. Madec, J.M. Molines ): advective bbl
13	!! 2.0 ! 2008-04 (S. Cravatte) add the i-, j- & k- trends computation
14	!! - ! 2009-11 (V. Garnier) Surface pressure gradient organization
15	!! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
16	!!----------------------------------------------------------------------
17
18	!!----------------------------------------------------------------------
19	!! tra_adv_tvd : update the tracer trend with the horizontal
20	!! and vertical advection trends using a TVD scheme
21	!! nonosc : compute monotonic tracer fluxes by a nonoscillatory
22	!! algorithm
23	!!----------------------------------------------------------------------
24	USE oce ! ocean dynamics and active tracers
25	USE dom_oce ! ocean space and time domain
26	USE trdmod_oce ! tracers trends
27	USE trdtra ! tracers trends
28	USE in_out_manager ! I/O manager
29	USE dynspg_oce ! choice/control of key cpp for surface pressure gradient
30	USE lib_mpp ! MPP library
31	USE lbclnk ! ocean lateral boundary condition (or mpp link)
32	USE diaptr ! poleward transport diagnostics
33	USE trc_oce ! share passive tracers/Ocean variables
34
35
36	IMPLICIT NONE
37	PRIVATE
38
39	PUBLIC tra_adv_tvd ! routine called by step.F90
40
41	LOGICAL :: l_trd ! flag to compute trends
42
43	!! * Control permutation of array indices
44	# include "oce_ftrans.h90"
45	# include "dom_oce_ftrans.h90"
46	# include "trc_oce_ftrans.h90"
47
48	!! * Substitutions
49	# include "domzgr_substitute.h90"
50	# include "vectopt_loop_substitute.h90"
51	!!----------------------------------------------------------------------
52	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
53	!! $Id$
54	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
55	!!----------------------------------------------------------------------
56	CONTAINS
57
58	SUBROUTINE tra_adv_tvd ( kt, cdtype, p2dt, pun, pvn, pwn, &
59	& ptb, ptn, pta, kjpt )
60	!!----------------------------------------------------------------------
61	!! * ROUTINE tra_adv_tvd *
62	!!
63	!! ** Purpose : Compute the now trend due to total advection of
64	!! tracers and add it to the general trend of tracer equations
65	!!
66	!! ** Method : TVD scheme, i.e. 2nd order centered scheme with
67	!! corrected flux (monotonic correction)
68	!! note: - this advection scheme needs a leap-frog time scheme
69	!!
70	!! ** Action : - update (pta) with the now advective tracer trends
71	!! - save the trends
72	!!----------------------------------------------------------------------
73	USE timing, ONLY: timing_start, timing_stop
74	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
75	USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace
76	USE wrk_nemo, ONLY: zwi => wrk_3d_12 , zwz => wrk_3d_13 ! 3D workspace
77
78	!! DCSE_NEMO: need additional directives for renamed module variables
79	!FTRANS zwx zwy zwi zwz :I :I :z
80
81	!
82	INTEGER , INTENT(in ) :: kt ! ocean time-step index
83	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
84	INTEGER , INTENT(in ) :: kjpt ! number of tracers
85	REAL(wp), DIMENSION( jpkorig ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
86
87	!! DCSE_NEMO: This style defeats ftrans
88	! REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
89	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
90	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
91
92	!FTRANS pun pvn pwn :I :I :z
93	!FTRANS ptb ptn pta :I :I :z :
94	REAL(wp), INTENT(in ) :: pun(jpi,jpj,jpkorig) ! ocean velocity component (u)
95	REAL(wp), INTENT(in ) :: pvn(jpi,jpj,jpkorig) ! ocean velocity component (v)
96	REAL(wp), INTENT(in ) :: pwn(jpi,jpj,jpkorig) ! ocean velocity component (w)
97	REAL(wp), INTENT(in ) :: ptb(jpi,jpj,jpkorig,kjpt) ! tracer fields (before)
98	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpkorig,kjpt) ! tracer fields (now)
99	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpkorig,kjpt) ! tracer trend
100
101	!
102	INTEGER :: ji, jj, jk, jn ! dummy loop indices
103	REAL(wp) :: z2dtt, zbtr, ztra ! local scalar
104	REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - -
105	REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - -
106	REAL(wp), DIMENSION (:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz
107	!FTRANS ztrdx ztrdy ztrdz :I :I :z
108
109	!!----------------------------------------------------------------------
110
111	CALL timing_start('tra_adv_tvd')
112
113	IF( wrk_in_use(3, 12,13) ) THEN
114	CALL ctl_stop('tra_adv_tvd: requested workspace arrays unavailable') ; RETURN
115	ENDIF
116
117	IF( kt == nit000 ) THEN
118	IF(lwp) WRITE(numout,*)
119	IF(lwp) WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype
120	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
121	!
122	l_trd = .FALSE.
123	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
124	ENDIF
125	!
126	IF( l_trd ) THEN
127	ALLOCATE( ztrdx(jpi,jpj,jpk) ) ; ztrdx(:,:,:) = 0.e0
128	ALLOCATE( ztrdy(jpi,jpj,jpk) ) ; ztrdy(:,:,:) = 0.e0
129	ALLOCATE( ztrdz(jpi,jpj,jpk) ) ; ztrdz(:,:,:) = 0.e0
130	END IF
131	!
132	#if defined key_z_first
133	DO jj = 1, jpj, 1
134	DO ji = 1, jpi, 1
135	zwi(ji,jj,1:mbkmax(ji,jj)) = 0.e0_wp
136	END DO
137	END DO
138	#else
139	zwi(:,:,1:jpk) = 0.e0_wp
140	#endif
141	!
142	! ! ===========
143	DO jn = 1, kjpt ! tracer loop
144	! ! ===========
145	! 1. Bottom value : flux set to zero
146	! ----------------------------------
147	! ARPDBG: FINISS was using jpk correct here given that it can be
148	! below the ocean floor??
149	#if defined key_z_first
150	DO jj = 1, jpj, 1
151	DO ji = 1, jpi, 1
152	zwx(ji,jj,mbkmax(ji,jj)) = 0.0_wp
153	zwy(ji,jj,mbkmax(ji,jj)) = 0.0_wp
154	zwz(ji,jj,mbkmax(ji,jj)) = 0.0_wp
155	zwi(ji,jj,mbkmax(ji,jj)) = 0.0_wp
156	END DO
157	END DO
158	#else
159	zwx(:,:,jpkf) = 0.e0 ; zwz(:,:,jpkf) = 0.e0
160	zwy(:,:,jpkf) = 0.e0 ; zwi(:,:,jpkf) = 0.e0
161	#endif
162	! 2. upstream advection with initial mass fluxes & intermediate update
163	! --------------------------------------------------------------------
164	! upstream tracer flux in the i and j direction
165	CALL timing_start('tvd_upstream')
166	#if defined key_z_first
167	DO jj = 1, jpjm1
168	DO ji = 1, jpim1
169	DO jk = 1, mbkmax(ji,jj)-1 ! jpkfm1
170	#else
171	DO jk = 1, jpkfm1
172	DO jj = 1, jpjm1
173	DO ji = 1, fs_jpim1 ! vector opt.
174	#endif
175	! upstream scheme
176	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
177	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
178	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
179	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
180	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) )
181	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) )
182	END DO
183	END DO
184	END DO
185	CALL timing_stop('tvd_upstream','section')
186
187	! upstream tracer flux in the k direction
188	! Surface value
189	CALL timing_start('tvd_upstreamk')
190	IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! volume variable
191	ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface
192	ENDIF
193	! Interior value
194	#if defined key_z_first
195	DO jj = 1, jpj
196	DO ji = 1, jpi
197	DO jk = 2, mbkmax(ji,jj)-1 ! jpkfm1
198	#else
199	DO jk = 2, jpkfm1
200	DO jj = 1, jpj
201	DO ji = 1, jpi
202	#endif
203	zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
204	zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
205	zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )
206	END DO
207	END DO
208	END DO
209	CALL timing_stop('tvd_upstreamk','section')
210
211	! total advective trend
212	CALL timing_start('tvd_tot')
213	#if defined key_z_first
214	DO jj = 2, jpjm1
215	DO ji = 2, jpim1
216	DO jk = 1, mbkmax(ji,jj)-1 ! jpkfm1
217	z2dtt = p2dt(jk)
218	#else
219	DO jk = 1, jpkfm1
220	z2dtt = p2dt(jk)
221	DO jj = 2, jpjm1
222	DO ji = fs_2, fs_jpim1 ! vector opt.
223	#endif
224	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
225	! total intermediate advective trends
226	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
227	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
228	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
229	! update and guess with monotonic sheme
230	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
231	zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)
232	END DO
233	END DO
234	END DO
235	CALL timing_stop('tvd_tot','section')
236
237	! ! Lateral boundary conditions on zwi (unchanged sign)
238	CALL timing_start('tvd_lbc')
239	CALL lbc_lnk( zwi, 'T', 1. )
240	CALL timing_stop('tvd_lbc','section')
241
242	! ! trend diagnostics (contribution of upstream fluxes)
243	IF( l_trd ) THEN
244	! store intermediate advective trends
245	ztrdx(:,:,1:jpkf) = zwx(:,:,1:jpkf) ; ztrdy(:,:,1:jpkf) = zwy(:,:,1:jpkf) ; ztrdz(:,:,1:jpkf) = zwz(:,:,1:jpkf)
246	END IF
247	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
248	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
249	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) )
250	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) )
251	ENDIF
252
253	! 3. antidiffusive flux : high order minus low order
254	! --------------------------------------------------
255	! antidiffusive flux on i and j
256	CALL timing_start('tvd_antidiff')
257	#if defined key_z_first
258	DO jj = 1, jpjm1
259	DO ji = 1, jpim1
260	DO jk = 1, mbkmax(ji,jj)-1 ! jpkfm1
261	#else
262	DO jk = 1, jpkfm1
263	DO jj = 1, jpjm1
264	DO ji = 1, fs_jpim1 ! vector opt.
265	#endif
266	zwx(ji,jj,jk) = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk)
267	zwy(ji,jj,jk) = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk)
268	END DO
269	END DO
270	END DO
271	CALL timing_stop('tvd_antidiff','section')
272
273	! antidiffusive flux on k
274	CALL timing_start('tvd_antidiffk')
275	#if defined key_z_first
276	DO jj = 1, jpj
277	DO ji = 1, jpi
278	zwz(ji,jj,1) = 0.e0 ! Surface value
279	DO jk = 2, mbkmax(ji,jj)-1 ! jpkfm1
280	#else
281	zwz(:,:,1) = 0.e0 ! Surface value
282	!
283	DO jk = 2, jpkfm1 ! Interior value
284	DO jj = 1, jpj
285	DO ji = 1, jpi
286	#endif
287	zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk)
288	END DO
289	END DO
290	END DO
291	CALL timing_stop('tvd_antidiffk','section')
292
293	CALL timing_start('tvd_lbc')
294	CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions
295	CALL lbc_lnk( zwz, 'W', 1. )
296	CALL timing_stop('tvd_lbc','section')
297
298	! 4. monotonicity algorithm
299	! -------------------------
300	CALL timing_start('tvd_nonosc')
301	CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
302	CALL timing_stop('tvd_nonosc','section')
303
304
305	! 5. final trend with corrected fluxes
306	! ------------------------------------
307	CALL timing_start('tvd_finaltr')
308	#if defined key_z_first
309	DO jj = 2, jpjm1
310	DO ji = 2, jpim1
311	DO jk = 1, mbkmax(ji,jj)-1 ! jpkfm1
312	#else
313	DO jk = 1, jpkfm1
314	DO jj = 2, jpjm1
315	DO ji = fs_2, fs_jpim1 ! vector opt.
316	#endif
317	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
318	! total advective trends
319	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
320	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
321	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
322	! add them to the general tracer trends
323	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
324	END DO
325	END DO
326	END DO
327
328	! ! trend diagnostics (contribution of upstream fluxes)
329	IF( l_trd ) THEN
330	ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed
331	ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed
332	ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed
333
334	CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, ztrdx, pun, ptn(:,:,:,jn) )
335	CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, ztrdy, pvn, ptn(:,:,:,jn) )
336	CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, ztrdz, pwn, ptn(:,:,:,jn) )
337	END IF
338	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
339	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
340	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:)
341	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:)
342	ENDIF
343	!
344	CALL timing_stop('tvd_finaltr','section')
345
346	END DO
347	!
348	IF( l_trd ) THEN
349	DEALLOCATE( ztrdx ) ; DEALLOCATE( ztrdy ) ; DEALLOCATE( ztrdz )
350	END IF
351	!
352	IF( wrk_not_released(3, 12,13) ) CALL ctl_stop('tra_adv_tvd: failed to release workspace arrays')
353	!
354	CALL timing_stop('tra_adv_tvd','section')
355
356	!! * Reset control of array index permutation
357	!FTRANS CLEAR
358	# include "oce_ftrans.h90"
359	# include "dom_oce_ftrans.h90"
360	# include "trc_oce_ftrans.h90"
361
362	END SUBROUTINE tra_adv_tvd
363
364
365	SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
366	!!---------------------------------------------------------------------
367	!! * ROUTINE nonosc *
368	!!
369	!! ** Purpose : compute monotonic tracer fluxes from the upstream
370	!! scheme and the before field by a nonoscillatory algorithm
371	!!
372	!! ** Method : ... ???
373	!! warning : pbef and paft must be masked, but the boundaries
374	!! conditions on the fluxes are not necessary zalezak (1979)
375	!! drange (1995) multi-dimensional forward-in-time and upstream-
376	!! in-space based differencing for fluid
377	!!----------------------------------------------------------------------
378	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
379	USE wrk_nemo, ONLY: zbetup => wrk_3d_8 , zbetdo => wrk_3d_9 ! 3D workspace
380	USE wrk_nemo, ONLY: zbup => wrk_3d_10 , zbdo => wrk_3d_11 ! - -
381	USE timing, ONLY: timing_start, timing_stop
382
383	!! DCSE_NEMO: need additional directives for renamed module variables
384	!FTRANS zbetup zbetdo zbup zbdo :I :I :z
385
386	!
387	REAL(wp), DIMENSION(jpkorig) , INTENT(in ) :: p2dt ! vertical profile of tracer time-step
388
389	!! DCSE_NEMO: This style defeats ftrans
390	! REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field
391	! REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions
392
393	!FTRANS pbef paft :I :I :z
394	!FTRANS paa pbb pcc :I :I :z
395	REAL(wp), INTENT(in ) :: pbef(jpi,jpj,jpkorig), paft(jpi,jpj,jpkorig) ! before & after field
396	REAL(wp), INTENT(inout) :: paa(jpi,jpj,jpkorig) ! monotonic fluxes in the 1st direction
397	REAL(wp), INTENT(inout) :: pbb(jpi,jpj,jpkorig) ! monotonic fluxes in the 2nd direction
398	REAL(wp), INTENT(inout) :: pcc(jpi,jpj,jpkorig) ! monotonic fluxes in the 3rd direction
399	!
400	INTEGER :: ji, jj, jk ! dummy loop indices
401	INTEGER :: ikm1 ! local integer
402	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt ! local scalars
403	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -
404	!!----------------------------------------------------------------------
405
406	IF( wrk_in_use(3, 8,9,10,11) ) THEN
407	CALL ctl_stop('nonosc: requested workspace array unavailable') ; RETURN
408	ENDIF
409
410	zbig = 1.e+40_wp
411	zrtrn = 1.e-15_wp
412	zbetup(:,:,jpkf) = 0._wp ; zbetdo(:,:,jpkf) = 0._wp
413
414
415	! Search local extrema
416	! --------------------
417	! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
418	zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ), &
419	& paft * tmask - zbig * ( 1.e0 - tmask ) )
420	zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ), &
421	& paft * tmask + zbig * ( 1.e0 - tmask ) )
422
423	#if defined key_z_first
424	DO jj = 2, jpjm1
425	DO ji = 2, jpim1
426	DO jk = 1, jpkfm1
427	ikm1 = MAX(jk-1,1)
428	z2dtt = p2dt(jk)
429	#else
430	DO jk = 1, jpkfm1
431	ikm1 = MAX(jk-1,1)
432	z2dtt = p2dt(jk)
433	DO jj = 2, jpjm1
434	DO ji = fs_2, fs_jpim1 ! vector opt.
435	#endif
436
437	! search maximum in neighbourhood
438	zup = MAX( zbup(ji ,jj ,jk ), &
439	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
440	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), &
441	& zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) )
442
443	! search minimum in neighbourhood
444	zdo = MIN( zbdo(ji ,jj ,jk ), &
445	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
446	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), &
447	& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
448
449	! positive part of the flux
450	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
451	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
452	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
453
454	! negative part of the flux
455	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
456	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
457	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
458
459	! up & down beta terms
460	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
461	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
462	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
463	END DO
464	END DO
465	END DO
466	CALL timing_start('tvd_lbc')
467	CALL lbc_lnk( zbetup, 'T', 1. ) ; CALL lbc_lnk( zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
468	CALL timing_stop('tvd_lbc','section')
469
470
471
472	! 3. monotonic flux in the i & j direction (paa & pbb)
473	! ----------------------------------------
474	#if defined key_z_first
475	DO jj = 2, jpjm1
476	DO ji = 2, jpim1
477	DO jk = 1, jpkfm1
478	#else
479	DO jk = 1, jpkfm1
480	DO jj = 2, jpjm1
481	DO ji = fs_2, fs_jpim1 ! vector opt.
482	#endif
483	zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
484	zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
485	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
486	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
487
488	zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
489	zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
490	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
491	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv )
492
493	! monotonic flux in the k direction, i.e. pcc
494	! -------------------------------------------
495	za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
496	zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
497	zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
498	pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1.e0 - zc) * zb )
499	END DO
500	END DO
501	END DO
502	CALL timing_start('tvd_lbc')
503	CALL lbc_lnk( paa, 'U', -1. ) ; CALL lbc_lnk( pbb, 'V', -1. ) ! lateral boundary condition (changed sign)
504	CALL timing_stop('tvd_lbc','section')
505
506	!
507	IF( wrk_not_released(3, 8,9,10,11) ) CALL ctl_stop('nonosc: failed to release workspace arrays')
508	!
509	END SUBROUTINE nonosc
510
511	!!======================================================================
512	END MODULE traadv_tvd

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

Download in other formats: