Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

traadv_tvd.F90 @ 1970

Last change on this file since 1970 was 1970, checked in by acc, 14 years ago
ticket #684 step 5: Add in changes from the trunk between revisions 1821 and 1879.
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 22.5 KB

Line
1	MODULE traadv_tvd
2	!!==============================================================================
3	!! * MODULE traadv_tvd *
4	!! Ocean active tracers: horizontal & vertical advective trend
5	!!==============================================================================
6	!! History : ! 95-12 (L. Mortier) Original code
7	!! ! 00-01 (H. Loukos) adapted to ORCA
8	!! ! 00-10 (MA Foujols E.Kestenare) include file not routine
9	!! ! 00-12 (E. Kestenare M. Levy) fix bug in trtrd indexes
10	!! ! 01-07 (E. Durand G. Madec) adaptation to ORCA config
11	!! 8.5 ! 02-06 (G. Madec) F90: Free form and module
12	!! 9.0 ! 04-01 (A. de Miranda, G. Madec, J.M. Molines ): advective bbl
13	!! 9.0 ! 08-04 (S. Cravatte) add the i-, j- & k- trends computation
14	!! " " ! 05-11 (V. Garnier) Surface pressure gradient organization
15	!!----------------------------------------------------------------------
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 ! ocean active tracers trends
27	USE trdmod_oce ! ocean variables trends
28	USE in_out_manager ! I/O manager
29	USE dynspg_oce ! choice/control of key cpp for surface pressure gradient
30	USE trabbl ! Advective term of BBL
31	USE lib_mpp
32	USE lbclnk ! ocean lateral boundary condition (or mpp link)
33	USE diaptr ! poleward transport diagnostics
34	USE prtctl ! Print control
35
36
37	IMPLICIT NONE
38	PRIVATE
39
40	PUBLIC tra_adv_tvd ! routine called by step.F90
41
42	!! * Substitutions
43	# include "domzgr_substitute.h90"
44	# include "vectopt_loop_substitute.h90"
45	!!----------------------------------------------------------------------
46	!! OPA 9.0 , LOCEAN-IPSL (2006)
47	!! $Id$
48	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
49	!!----------------------------------------------------------------------
50
51	CONTAINS
52
53	SUBROUTINE tra_adv_tvd( kt, pun, pvn, pwn )
54	!!----------------------------------------------------------------------
55	!! * ROUTINE tra_adv_tvd *
56	!!
57	!! ** Purpose : Compute the now trend due to total advection of
58	!! tracers and add it to the general trend of tracer equations
59	!!
60	!! ** Method : TVD scheme, i.e. 2nd order centered scheme with
61	!! corrected flux (monotonic correction)
62	!! note: - this advection scheme needs a leap-frog time scheme
63	!!
64	!! ** Action : - update (ta,sa) with the now advective tracer trends
65	!! - save the trends in (ztrdt,ztrds) ('key_trdtra')
66	!!----------------------------------------------------------------------
67	USE oce , ztrdt => ua ! use ua as workspace
68	USE oce , ztrds => va ! use va as workspace
69	!!
70	INTEGER , INTENT(in) :: kt ! ocean time-step index
71	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! ocean velocity u-component
72	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! ocean velocity v-component
73	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! ocean velocity w-component
74	!!
75	INTEGER :: ji, jj, jk ! dummy loop indices
76	REAL(wp) :: & ! temporary scalar
77	ztat, zsat, & ! " "
78	z_hdivn_x, z_hdivn_y, z_hdivn
79	REAL(wp) :: &
80	z2dtt, zbtr, zeu, zev, & ! temporary scalar
81	zew, z2, zbtr1, & ! temporary scalar
82	zfp_ui, zfp_vj, zfp_wk, & ! " "
83	zfm_ui, zfm_vj, zfm_wk ! " "
84	REAL(wp), DIMENSION (jpi,jpj,jpk) :: zti, ztu, ztv, ztw ! temporary workspace
85	REAL(wp), DIMENSION (jpi,jpj,jpk) :: zsi, zsu, zsv, zsw ! " "
86	!!----------------------------------------------------------------------
87
88	zti(:,:,:) = 0.e0 ; zsi(:,:,:) = 0.e0
89
90	IF( kt == nit000 .AND. lwp ) THEN
91	WRITE(numout,*)
92	WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme'
93	WRITE(numout,*) '~~~~~~~~~~~'
94	ENDIF
95
96	IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1.
97	ELSE ; z2 = 2.
98	ENDIF
99
100	! 1. Bottom value : flux set to zero
101	! ---------------
102	ztu(:,:,jpk) = 0.e0 ; zsu(:,:,jpk) = 0.e0
103	ztv(:,:,jpk) = 0.e0 ; zsv(:,:,jpk) = 0.e0
104	ztw(:,:,jpk) = 0.e0 ; zsw(:,:,jpk) = 0.e0
105	zti(:,:,jpk) = 0.e0 ; zsi(:,:,jpk) = 0.e0
106
107
108	! 2. upstream advection with initial mass fluxes & intermediate update
109	! --------------------------------------------------------------------
110	! upstream tracer flux in the i and j direction
111	DO jk = 1, jpkm1
112	DO jj = 1, jpjm1
113	DO ji = 1, fs_jpim1 ! vector opt.
114	zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
115	zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
116	! upstream scheme
117	zfp_ui = zeu + ABS( zeu )
118	zfm_ui = zeu - ABS( zeu )
119	zfp_vj = zev + ABS( zev )
120	zfm_vj = zev - ABS( zev )
121	ztu(ji,jj,jk) = zfp_ui * tb(ji,jj,jk) + zfm_ui * tb(ji+1,jj ,jk)
122	ztv(ji,jj,jk) = zfp_vj * tb(ji,jj,jk) + zfm_vj * tb(ji ,jj+1,jk)
123	zsu(ji,jj,jk) = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj ,jk)
124	zsv(ji,jj,jk) = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji ,jj+1,jk)
125	END DO
126	END DO
127	END DO
128
129	! upstream tracer flux in the k direction
130	! Surface value
131	IF( lk_vvl ) THEN
132	! variable volume : flux set to zero
133	ztw(:,:,1) = 0.e0
134	zsw(:,:,1) = 0.e0
135	ELSE
136	! free surface-constant volume
137	DO jj = 1, jpj
138	DO ji = 1, jpi
139	zew = e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,1)
140	ztw(ji,jj,1) = zew * tb(ji,jj,1)
141	zsw(ji,jj,1) = zew * sb(ji,jj,1)
142	END DO
143	END DO
144	ENDIF
145
146	! Interior value
147	DO jk = 2, jpkm1
148	DO jj = 1, jpj
149	DO ji = 1, jpi
150	zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
151	zfp_wk = zew + ABS( zew )
152	zfm_wk = zew - ABS( zew )
153	ztw(ji,jj,jk) = zfp_wk * tb(ji,jj,jk) + zfm_wk * tb(ji,jj,jk-1)
154	zsw(ji,jj,jk) = zfp_wk * sb(ji,jj,jk) + zfm_wk * sb(ji,jj,jk-1)
155	END DO
156	END DO
157	END DO
158
159	! total advective trend
160	DO jk = 1, jpkm1
161	z2dtt = z2 * rdttra(jk)
162	DO jj = 2, jpjm1
163	DO ji = fs_2, fs_jpim1 ! vector opt.
164	zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
165	! total intermediate advective trends
166	ztat = - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk ) &
167	& + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk ) &
168	& + ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) * zbtr
169	zsat = - ( zsu(ji,jj,jk) - zsu(ji-1,jj ,jk ) &
170	& + zsv(ji,jj,jk) - zsv(ji ,jj-1,jk ) &
171	& + zsw(ji,jj,jk) - zsw(ji ,jj ,jk+1) ) * zbtr
172	! update and guess with monotonic sheme
173	ta(ji,jj,jk) = ta(ji,jj,jk) + ztat
174	sa(ji,jj,jk) = sa(ji,jj,jk) + zsat
175	zti (ji,jj,jk) = ( tb(ji,jj,jk) + z2dtt * ztat ) * tmask(ji,jj,jk)
176	zsi (ji,jj,jk) = ( sb(ji,jj,jk) + z2dtt * zsat ) * tmask(ji,jj,jk)
177	END DO
178	END DO
179	END DO
180
181	! "zonal" mean advective heat and salt transport
182	IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
183	pht_adv(:) = ptr_vj( ztv(:,:,:) )
184	pst_adv(:) = ptr_vj( zsv(:,:,:) )
185	ENDIF
186
187	! Save the intermediate i / j / k advective trends for diagnostics
188	! -------------------------------------------------------------------
189	! Warning : We should use zun instead of un in the computations below, but we
190	! also use hdivn which is computed with un, vn (check ???). So we use un, vn
191	! for consistency. Results are therefore approximate with key_trabbl_adv.
192
193	IF( l_trdtra ) THEN
194	ztrdt(:,:,:) = 0.e0 ; ztrds(:,:,:) = 0.e0
195	!
196	! T/S ZONAL advection trends
197	DO jk = 1, jpkm1
198	DO jj = 2, jpjm1
199	DO ji = fs_2, fs_jpim1 ! vector opt.
200	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
201	ztrdt(ji,jj,jk) = - ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr
202	ztrds(ji,jj,jk) = - ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) ) * zbtr
203	END DO
204	END DO
205	END DO
206	CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt) ! save the trends
207	!
208	! T/S MERIDIONAL advection trends
209	DO jk = 1, jpkm1
210	DO jj = 2, jpjm1
211	DO ji = fs_2, fs_jpim1 ! vector opt.
212	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
213	ztrdt(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zbtr
214	ztrds(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji,jj-1,jk) ) * zbtr
215	END DO
216	END DO
217	END DO
218	CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt) ! save the trends
219	!
220	! T/S VERTICAL advection trends
221	DO jk = 1, jpkm1
222	DO jj = 2, jpjm1
223	DO ji = fs_2, fs_jpim1 ! vector opt.
224	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
225	ztrdt(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
226	ztrds(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr
227	END DO
228	END DO
229	END DO
230	CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt) ! save the trends
231	!
232	ENDIF
233
234	! Lateral boundary conditions on zti, zsi (unchanged sign)
235	CALL lbc_lnk( zti, 'T', 1. )
236	CALL lbc_lnk( zsi, 'T', 1. )
237
238
239	! 3. antidiffusive flux : high order minus low order
240	! --------------------------------------------------
241	! antidiffusive flux on i and j
242	DO jk = 1, jpkm1
243	DO jj = 1, jpjm1
244	DO ji = 1, fs_jpim1 ! vector opt.
245	zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
246	zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
247	ztu(ji,jj,jk) = zeu * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) ) - ztu(ji,jj,jk)
248	zsu(ji,jj,jk) = zeu * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) ) - zsu(ji,jj,jk)
249	ztv(ji,jj,jk) = zev * ( tn(ji,jj,jk) + tn(ji,jj+1,jk) ) - ztv(ji,jj,jk)
250	zsv(ji,jj,jk) = zev * ( sn(ji,jj,jk) + sn(ji,jj+1,jk) ) - zsv(ji,jj,jk)
251	END DO
252	END DO
253	END DO
254
255	! antidiffusive flux on k
256	! Surface value
257	ztw(:,:,1) = 0.e0
258	zsw(:,:,1) = 0.e0
259
260	! Interior value
261	DO jk = 2, jpkm1
262	DO jj = 1, jpj
263	DO ji = 1, jpi
264	zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
265	ztw(ji,jj,jk) = zew * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) - ztw(ji,jj,jk)
266	zsw(ji,jj,jk) = zew * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) - zsw(ji,jj,jk)
267	END DO
268	END DO
269	END DO
270
271	! Lateral bondary conditions
272	CALL lbc_lnk( ztu, 'U', -1. ) ; CALL lbc_lnk( zsu, 'U', -1. )
273	CALL lbc_lnk( ztv, 'V', -1. ) ; CALL lbc_lnk( zsv, 'V', -1. )
274	CALL lbc_lnk( ztw, 'W', 1. ) ; CALL lbc_lnk( zsw, 'W', 1. )
275
276	! 4. monotonicity algorithm
277	! -------------------------
278	CALL nonosc( tb, ztu, ztv, ztw, zti, z2 )
279	CALL nonosc( sb, zsu, zsv, zsw, zsi, z2 )
280
281
282	! 5. final trend with corrected fluxes
283	! ------------------------------------
284	DO jk = 1, jpkm1
285	DO jj = 2, jpjm1
286	DO ji = fs_2, fs_jpim1 ! vector opt.
287	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
288	! total advective trends
289	ztat = - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk ) &
290	& + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk ) &
291	& + ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) * zbtr
292	zsat = - ( zsu(ji,jj,jk) - zsu(ji-1,jj ,jk ) &
293	& + zsv(ji,jj,jk) - zsv(ji ,jj-1,jk ) &
294	& + zsw(ji,jj,jk) - zsw(ji ,jj ,jk+1) ) * zbtr
295	! add them to the general tracer trends
296	ta(ji,jj,jk) = ta(ji,jj,jk) + ztat
297	sa(ji,jj,jk) = sa(ji,jj,jk) + zsat
298	END DO
299	END DO
300	END DO
301
302
303	! Save the advective trends for diagnostics
304	! --------------------------------------------
305
306	IF( l_trdtra ) THEN
307	ztrdt(:,:,:) = 0.e0 ; ztrds(:,:,:) = 0.e0
308	!
309	! T/S ZONAL advection trends
310	DO jk = 1, jpkm1
311	DO jj = 2, jpjm1
312	DO ji = fs_2, fs_jpim1 ! vector opt.
313	!-- Compute zonal divergence by splitting hdivn (see divcur.F90)
314	! N.B. This computation is not valid along OBCs (if any)
315	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
316	z_hdivn_x = ( e2u(ji ,jj) * fse3u(ji ,jj,jk) * pun(ji ,jj,jk) &
317	& - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr
318	!-- Compute T/S zonal advection trends
319	ztrdt(ji,jj,jk) = - ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_x
320	ztrds(ji,jj,jk) = - ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_x
321	END DO
322	END DO
323	END DO
324	CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt, cnbpas='bis') ! <<< ADD TO PREVIOUSLY COMPUTED
325	!
326	! T/S MERIDIONAL advection trends
327	DO jk = 1, jpkm1
328	DO jj = 2, jpjm1
329	DO ji = fs_2, fs_jpim1 ! vector opt.
330	!-- Compute merid. divergence by splitting hdivn (see divcur.F90)
331	! N.B. This computation is not valid along OBCs (if any)
332	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
333	z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) &
334	& - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr
335	!-- Compute T/S meridional advection trends
336	ztrdt(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_y
337	ztrds(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji,jj-1,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_y
338	END DO
339	END DO
340	END DO
341	CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt, cnbpas='bis') ! <<< ADD TO PREVIOUSLY COMPUTED
342	!
343	! T/S VERTICAL advection trends
344	DO jk = 1, jpkm1
345	DO jj = 2, jpjm1
346	DO ji = fs_2, fs_jpim1 ! vector opt.
347	zbtr1 = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
348	#if defined key_zco
349	zbtr = zbtr1
350	z_hdivn_x = e2u(ji,jj)pun(ji,jj,jk) - e2u(ji-1,jj)pun(ji-1,jj,jk)
351	z_hdivn_y = e1v(ji,jj)pvn(ji,jj,jk) - e1v(ji,jj-1)pvn(ji,jj-1,jk)
352	#else
353	zbtr = zbtr1 / fse3t(ji,jj,jk)
354	z_hdivn_x = e2u(ji,jj)fse3u(ji,jj,jk)pun(ji,jj,jk) - e2u(ji-1,jj)fse3u(ji-1,jj,jk)pun(ji-1,jj,jk)
355	z_hdivn_y = e1v(ji,jj)fse3v(ji,jj,jk)pvn(ji,jj,jk) - e1v(ji,jj-1)fse3v(ji,jj-1,jk)pvn(ji,jj-1,jk)
356	#endif
357	z_hdivn = (z_hdivn_x + z_hdivn_y) * zbtr
358	zbtr = zbtr1 / fse3t(ji,jj,jk)
359	ztrdt(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr - tn(ji,jj,jk) * z_hdivn
360	ztrds(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr - sn(ji,jj,jk) * z_hdivn
361	END DO
362	END DO
363	END DO
364	CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt, cnbpas='bis') ! <<< ADD TO PREVIOUSLY COMPUTED
365	!
366	ENDIF
367
368	IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' tvd adv - Ta: ', mask1=tmask, &
369	& tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
370
371	! "zonal" mean advective heat and salt transport
372	IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
373	pht_adv(:) = ptr_vj( ztv(:,:,:) ) + pht_adv(:)
374	pst_adv(:) = ptr_vj( zsv(:,:,:) ) + pst_adv(:)
375	ENDIF
376	!
377	END SUBROUTINE tra_adv_tvd
378
379
380	SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, prdt )
381	!!---------------------------------------------------------------------
382	!! * ROUTINE nonosc *
383	!!
384	!! ** Purpose : compute monotonic tracer fluxes from the upstream
385	!! scheme and the before field by a nonoscillatory algorithm
386	!!
387	!! ** Method : ... ???
388	!! warning : pbef and paft must be masked, but the boundaries
389	!! conditions on the fluxes are not necessary zalezak (1979)
390	!! drange (1995) multi-dimensional forward-in-time and upstream-
391	!! in-space based differencing for fluid
392	!!----------------------------------------------------------------------
393	REAL(wp), INTENT( in ) :: prdt ! ???
394	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT( inout ) :: &
395	pbef, & ! before field
396	paft, & ! after field
397	paa, & ! monotonic flux in the i direction
398	pbb, & ! monotonic flux in the j direction
399	pcc ! monotonic flux in the k direction
400	!!
401	INTEGER :: ji, jj, jk ! dummy loop indices
402	INTEGER :: ikm1
403	REAL(wp), DIMENSION (jpi,jpj,jpk) :: zbetup, zbetdo
404	REAL(wp), DIMENSION (jpi,jpj,jpk) :: zbup, zbdo
405	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt
406	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv
407	REAL(wp) :: zup, zdo
408	!!----------------------------------------------------------------------
409
410	zbig = 1.e+40
411	zrtrn = 1.e-15
412	zbetup(:,:,jpk) = 0.e0 ; zbetdo(:,:,jpk) = 0.e0
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	DO jk = 1, jpkm1
424	ikm1 = MAX(jk-1,1)
425	z2dtt = prdt * rdttra(jk)
426	DO jj = 2, jpjm1
427	DO ji = fs_2, fs_jpim1 ! vector opt.
428
429	! search maximum in neighbourhood
430	zup = MAX( zbup(ji ,jj ,jk ), &
431	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
432	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), &
433	& zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) )
434
435	! search minimum in neighbourhood
436	zdo = MIN( zbdo(ji ,jj ,jk ), &
437	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
438	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), &
439	& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
440
441	! positive part of the flux
442	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
443	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
444	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
445
446	! negative part of the flux
447	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
448	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
449	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
450
451	! up & down beta terms
452	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
453	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
454	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
455	END DO
456	END DO
457	END DO
458
459	! lateral boundary condition on zbetup & zbetdo (unchanged sign)
460	CALL lbc_lnk( zbetup, 'T', 1. )
461	CALL lbc_lnk( zbetdo, 'T', 1. )
462
463
464	! 3. monotonic flux in the i & j direction (paa & pbb)
465	! ----------------------------------------
466	DO jk = 1, jpkm1
467	DO jj = 2, jpjm1
468	DO ji = fs_2, fs_jpim1 ! vector opt.
469	zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
470	zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
471	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
472	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
473
474	zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
475	zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
476	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
477	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv )
478
479	! monotonic flux in the k direction, i.e. pcc
480	! -------------------------------------------
481	za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
482	zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
483	zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
484	pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1.e0 - zc) * zb )
485	END DO
486	END DO
487	END DO
488
489	! lateral boundary condition on paa, pbb, pcc
490	CALL lbc_lnk( paa, 'U', -1. ) ! changed sign
491	CALL lbc_lnk( pbb, 'V', -1. ) ! changed sign
492	!
493	END SUBROUTINE nonosc
494
495	!!======================================================================
496	END MODULE traadv_tvd

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/DEV_r1784_mid_year_merge_2010/NEMO/OPA_SRC/TRA/traadv_tvd.F90 @ 1970

Download in other formats: