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

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 20.8 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 wrk_nemo, ONLY: wrk_in_use, wrk_not_released
74	USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace
75	USE wrk_nemo, ONLY: zwi => wrk_3d_12 , zwz => wrk_3d_13 ! 3D workspace
76
77	!! DCSE_NEMO: need additional directives for renamed module variables
78	!FTRANS zwx zwy zwi zwz :I :I :z
79
80	!
81	INTEGER , INTENT(in ) :: kt ! ocean time-step index
82	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
83	INTEGER , INTENT(in ) :: kjpt ! number of tracers
84	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
85
86	!! DCSE_NEMO: This style defeats ftrans
87	! REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
88	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
89	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
90
91	!FTRANS pun pvn pwn :I :I :z
92	!FTRANS ptb ptn pta :I :I :z :
93	REAL(wp), INTENT(in ) :: pun(jpi,jpj,jpk) ! ocean velocity component (u)
94	REAL(wp), INTENT(in ) :: pvn(jpi,jpj,jpk) ! ocean velocity component (v)
95	REAL(wp), INTENT(in ) :: pwn(jpi,jpj,jpk) ! ocean velocity component (w)
96	REAL(wp), INTENT(in ) :: ptb(jpi,jpj,jpk,kjpt) ! tracer fields (before)
97	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpk,kjpt) ! tracer fields (now)
98	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpk,kjpt) ! tracer trend
99
100	!
101	INTEGER :: ji, jj, jk, jn ! dummy loop indices
102	REAL(wp) :: z2dtt, zbtr, ztra ! local scalar
103	REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - -
104	REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - -
105	REAL(wp), DIMENSION (:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz
106	!FTRANS ztrdx ztrdy ztrdz :I :I :z
107
108	!!----------------------------------------------------------------------
109
110	IF( wrk_in_use(3, 12,13) ) THEN
111	CALL ctl_stop('tra_adv_tvd: requested workspace arrays unavailable') ; RETURN
112	ENDIF
113
114	IF( kt == nit000 ) THEN
115	IF(lwp) WRITE(numout,*)
116	IF(lwp) WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype
117	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
118	!
119	l_trd = .FALSE.
120	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
121	ENDIF
122	!
123	IF( l_trd ) THEN
124	ALLOCATE( ztrdx(jpi,jpj,jpk) ) ; ztrdx(:,:,:) = 0.e0
125	ALLOCATE( ztrdy(jpi,jpj,jpk) ) ; ztrdy(:,:,:) = 0.e0
126	ALLOCATE( ztrdz(jpi,jpj,jpk) ) ; ztrdz(:,:,:) = 0.e0
127	END IF
128	!
129	zwi(:,:,:) = 0.e0
130	!
131	! ! ===========
132	DO jn = 1, kjpt ! tracer loop
133	! ! ===========
134	! 1. Bottom value : flux set to zero
135	! ----------------------------------
136	zwx(:,:,jpk) = 0.e0 ; zwz(:,:,jpk) = 0.e0
137	zwy(:,:,jpk) = 0.e0 ; zwi(:,:,jpk) = 0.e0
138
139	! 2. upstream advection with initial mass fluxes & intermediate update
140	! --------------------------------------------------------------------
141	! upstream tracer flux in the i and j direction
142	#if defined key_z_first
143	DO jj = 1, jpjm1
144	DO ji = 1, jpim1
145	DO jk = 1, jpkm1
146	#else
147	DO jk = 1, jpkm1
148	DO jj = 1, jpjm1
149	DO ji = 1, fs_jpim1 ! vector opt.
150	#endif
151	! upstream scheme
152	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
153	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
154	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
155	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
156	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) )
157	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) )
158	END DO
159	END DO
160	END DO
161
162	! upstream tracer flux in the k direction
163	! Surface value
164	IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! volume variable
165	ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface
166	ENDIF
167	! Interior value
168	#if defined key_z_first
169	DO jj = 1, jpj
170	DO ji = 1, jpi
171	DO jk = 2, jpkm1
172	#else
173	DO jk = 2, jpkm1
174	DO jj = 1, jpj
175	DO ji = 1, jpi
176	#endif
177	zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
178	zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
179	zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )
180	END DO
181	END DO
182	END DO
183
184	! total advective trend
185	#if defined key_z_first
186	DO jj = 2, jpjm1
187	DO ji = 2, jpim1
188	DO jk = 1, jpkm1
189	z2dtt = p2dt(jk)
190	#else
191	DO jk = 1, jpkm1
192	z2dtt = p2dt(jk)
193	DO jj = 2, jpjm1
194	DO ji = fs_2, fs_jpim1 ! vector opt.
195	#endif
196	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
197	! total intermediate advective trends
198	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
199	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
200	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
201	! update and guess with monotonic sheme
202	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
203	zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)
204	END DO
205	END DO
206	END DO
207	! ! Lateral boundary conditions on zwi (unchanged sign)
208	CALL lbc_lnk( zwi, 'T', 1. )
209
210	! ! trend diagnostics (contribution of upstream fluxes)
211	IF( l_trd ) THEN
212	! store intermediate advective trends
213	ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)
214	END IF
215	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
216	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
217	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) )
218	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) )
219	ENDIF
220
221	! 3. antidiffusive flux : high order minus low order
222	! --------------------------------------------------
223	! antidiffusive flux on i and j
224	#if defined key_z_first
225	DO jj = 1, jpjm1
226	DO ji = 1, jpim1
227	DO jk = 1, jpkm1
228	#else
229	DO jk = 1, jpkm1
230	DO jj = 1, jpjm1
231	DO ji = 1, fs_jpim1 ! vector opt.
232	#endif
233	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)
234	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)
235	END DO
236	END DO
237	END DO
238
239	! antidiffusive flux on k
240	#if defined key_z_first
241	DO jj = 1, jpj
242	DO ji = 1, jpi
243	zwz(ji,jj,1) = 0.e0 ! Surface value
244	DO jk = 2, jpkm1
245	#else
246	zwz(:,:,1) = 0.e0 ! Surface value
247	!
248	DO jk = 2, jpkm1 ! Interior value
249	DO jj = 1, jpj
250	DO ji = 1, jpi
251	#endif
252	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)
253	END DO
254	END DO
255	END DO
256	CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions
257	CALL lbc_lnk( zwz, 'W', 1. )
258
259	! 4. monotonicity algorithm
260	! -------------------------
261	CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
262
263
264	! 5. final trend with corrected fluxes
265	! ------------------------------------
266	#if defined key_z_first
267	DO jj = 2, jpjm1
268	DO ji = 2, jpim1
269	DO jk = 1, jpkm1
270	#else
271	DO jk = 1, jpkm1
272	DO jj = 2, jpjm1
273	DO ji = fs_2, fs_jpim1 ! vector opt.
274	#endif
275	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
276	! total advective trends
277	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
278	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
279	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
280	! add them to the general tracer trends
281	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
282	END DO
283	END DO
284	END DO
285
286	! ! trend diagnostics (contribution of upstream fluxes)
287	IF( l_trd ) THEN
288	ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed
289	ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed
290	ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed
291
292	CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, ztrdx, pun, ptn(:,:,:,jn) )
293	CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, ztrdy, pvn, ptn(:,:,:,jn) )
294	CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, ztrdz, pwn, ptn(:,:,:,jn) )
295	END IF
296	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
297	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
298	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:)
299	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:)
300	ENDIF
301	!
302	END DO
303	!
304	IF( l_trd ) THEN
305	DEALLOCATE( ztrdx ) ; DEALLOCATE( ztrdy ) ; DEALLOCATE( ztrdz )
306	END IF
307	!
308	IF( wrk_not_released(3, 12,13) ) CALL ctl_stop('tra_adv_tvd: failed to release workspace arrays')
309	!
310
311	!! * Reset control of array index permutation
312	!FTRANS CLEAR
313	# include "oce_ftrans.h90"
314	# include "dom_oce_ftrans.h90"
315	# include "trc_oce_ftrans.h90"
316
317	END SUBROUTINE tra_adv_tvd
318
319
320	SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
321	!!---------------------------------------------------------------------
322	!! * ROUTINE nonosc *
323	!!
324	!! ** Purpose : compute monotonic tracer fluxes from the upstream
325	!! scheme and the before field by a nonoscillatory algorithm
326	!!
327	!! ** Method : ... ???
328	!! warning : pbef and paft must be masked, but the boundaries
329	!! conditions on the fluxes are not necessary zalezak (1979)
330	!! drange (1995) multi-dimensional forward-in-time and upstream-
331	!! in-space based differencing for fluid
332	!!----------------------------------------------------------------------
333	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
334	USE wrk_nemo, ONLY: zbetup => wrk_3d_8 , zbetdo => wrk_3d_9 ! 3D workspace
335	USE wrk_nemo, ONLY: zbup => wrk_3d_10 , zbdo => wrk_3d_11 ! - -
336
337	!! DCSE_NEMO: need additional directives for renamed module variables
338	!FTRANS zbetup zbetdo zbup zbdo :I :I :z
339
340	!
341	REAL(wp), DIMENSION(jpk) , INTENT(in ) :: p2dt ! vertical profile of tracer time-step
342
343	!! DCSE_NEMO: This style defeats ftrans
344	! REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field
345	! REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions
346
347	!FTRANS pbef paft :I :I :z
348	!FTRANS paa pbb pcc :I :I :z
349	REAL(wp), INTENT(in ) :: pbef(jpi,jpj,jpk), paft(jpi,jpj,jpk) ! before & after field
350	REAL(wp), INTENT(inout) :: paa(jpi,jpj,jpk) ! monotonic fluxes in the 1st direction
351	REAL(wp), INTENT(inout) :: pbb(jpi,jpj,jpk) ! monotonic fluxes in the 2nd direction
352	REAL(wp), INTENT(inout) :: pcc(jpi,jpj,jpk) ! monotonic fluxes in the 3rd direction
353	!
354	INTEGER :: ji, jj, jk ! dummy loop indices
355	INTEGER :: ikm1 ! local integer
356	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt ! local scalars
357	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -
358	!!----------------------------------------------------------------------
359
360	IF( wrk_in_use(3, 8,9,10,11) ) THEN
361	CALL ctl_stop('nonosc: requested workspace array unavailable') ; RETURN
362	ENDIF
363
364	zbig = 1.e+40_wp
365	zrtrn = 1.e-15_wp
366	zbetup(:,:,jpk) = 0._wp ; zbetdo(:,:,jpk) = 0._wp
367
368
369	! Search local extrema
370	! --------------------
371	! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
372	zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ), &
373	& paft * tmask - zbig * ( 1.e0 - tmask ) )
374	zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ), &
375	& paft * tmask + zbig * ( 1.e0 - tmask ) )
376
377	#if defined key_z_first
378	DO jj = 2, jpjm1
379	DO ji = 2, jpim1
380	DO jk = 1, jpkm1
381	ikm1 = MAX(jk-1,1)
382	z2dtt = p2dt(jk)
383	#else
384	DO jk = 1, jpkm1
385	ikm1 = MAX(jk-1,1)
386	z2dtt = p2dt(jk)
387	DO jj = 2, jpjm1
388	DO ji = fs_2, fs_jpim1 ! vector opt.
389	#endif
390
391	! search maximum in neighbourhood
392	zup = MAX( zbup(ji ,jj ,jk ), &
393	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
394	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), &
395	& zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) )
396
397	! search minimum in neighbourhood
398	zdo = MIN( zbdo(ji ,jj ,jk ), &
399	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
400	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), &
401	& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
402
403	! positive part of the flux
404	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
405	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
406	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
407
408	! negative part of the flux
409	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
410	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
411	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
412
413	! up & down beta terms
414	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
415	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
416	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
417	END DO
418	END DO
419	END DO
420	CALL lbc_lnk( zbetup, 'T', 1. ) ; CALL lbc_lnk( zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
421
422
423
424	! 3. monotonic flux in the i & j direction (paa & pbb)
425	! ----------------------------------------
426	#if defined key_z_first
427	DO jj = 2, jpjm1
428	DO ji = 2, jpim1
429	DO jk = 1, jpkm1
430	#else
431	DO jk = 1, jpkm1
432	DO jj = 2, jpjm1
433	DO ji = fs_2, fs_jpim1 ! vector opt.
434	#endif
435	zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
436	zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
437	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
438	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
439
440	zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
441	zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
442	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
443	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv )
444
445	! monotonic flux in the k direction, i.e. pcc
446	! -------------------------------------------
447	za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
448	zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
449	zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
450	pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1.e0 - zc) * zb )
451	END DO
452	END DO
453	END DO
454	CALL lbc_lnk( paa, 'U', -1. ) ; CALL lbc_lnk( pbb, 'V', -1. ) ! lateral boundary condition (changed sign)
455	!
456	IF( wrk_not_released(3, 8,9,10,11) ) CALL ctl_stop('nonosc: failed to release workspace arrays')
457	!
458	END SUBROUTINE nonosc
459
460	!!======================================================================
461	END MODULE traadv_tvd

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

Download in other formats: