Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

traadv_tvd.F90 @ 6736

Last change on this file since 6736 was 6736, checked in by jamesharle, 8 years ago
FASTNEt code modifications
Property svn:keywords set to `Id`
File size: 19.6 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	!! 3.4.1! 2013-05 (H. Liu) add vertical PPM option (vppm)
17	!!----------------------------------------------------------------------
18
19	!!----------------------------------------------------------------------
20	!! tra_adv_tvd : update the tracer trend with the horizontal
21	!! and vertical advection trends using a TVD scheme
22	!! nonosc : compute monotonic tracer fluxes by a nonoscillatory
23	!! algorithm
24	!!----------------------------------------------------------------------
25	USE oce ! ocean dynamics and active tracers
26	USE dom_oce ! ocean space and time domain
27	USE trdmod_oce ! tracers trends
28	USE trdtra ! tracers trends
29	USE in_out_manager ! I/O manager
30	USE dynspg_oce ! choice/control of key cpp for surface pressure gradient
31	USE lib_mpp ! MPP library
32	USE lbclnk ! ocean lateral boundary condition (or mpp link)
33	USE diaptr ! poleward transport diagnostics
34	USE trc_oce ! share passive tracers/Ocean variables
35	USE wrk_nemo ! Memory Allocation
36	USE timing ! Timing
37	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
38
39	#if defined key_vppm
40	USE traadv_vppm ! vertical ppm scheme
41	!N.B. naac = 1 will not work here
42	#endif
43
44	IMPLICIT NONE
45	PRIVATE
46
47	PUBLIC tra_adv_tvd ! routine called by step.F90
48
49	LOGICAL :: l_trd ! flag to compute trends
50
51	!! * Substitutions
52	# include "domzgr_substitute.h90"
53	# include "vectopt_loop_substitute.h90"
54	!!----------------------------------------------------------------------
55	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
56	!! $Id$
57	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
58	!!----------------------------------------------------------------------
59	CONTAINS
60
61	SUBROUTINE tra_adv_tvd ( kt, kit000, cdtype, p2dt, pun, pvn, pwn, &
62	& ptb, ptn, pta, kjpt )
63	!!----------------------------------------------------------------------
64	!! * ROUTINE tra_adv_tvd *
65	!!
66	!! ** Purpose : Compute the now trend due to total advection of
67	!! tracers and add it to the general trend of tracer equations
68	!!
69	!! ** Method : TVD scheme, i.e. 2nd order centered scheme with
70	!! corrected flux (monotonic correction)
71	!! note: - this advection scheme needs a leap-frog time scheme
72	!!
73	!! ** Action : - update (pta) with the now advective tracer trends
74	!! - save the trends
75	!!----------------------------------------------------------------------
76	USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace
77	!
78	INTEGER , INTENT(in ) :: kt ! ocean time-step index
79	INTEGER , INTENT(in ) :: kit000 ! first time step index
80	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
81	INTEGER , INTENT(in ) :: kjpt ! number of tracers
82	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
83	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
84	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
85	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
86	!
87	INTEGER :: ji, jj, jk, jn ! dummy loop indices
88	REAL(wp) :: z2dtt, zbtr, ztra ! local scalar
89	REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - -
90	REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - -
91	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwz
92	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz
93	#if defined key_vppm
94	REAL(wp), POINTER, DIMENSION(:,:,:) :: hflux
95	#endif
96	!!----------------------------------------------------------------------
97	!
98	IF( nn_timing == 1 ) CALL timing_start('tra_adv_tvd')
99	!
100	CALL wrk_alloc( jpi, jpj, jpk, zwi, zwz )
101	#if defined key_vppm
102	CALL wrk_alloc( jpi, jpj, jpk, hflux )
103	#endif
104	!
105	IF( kt == kit000 ) THEN
106	IF(lwp) WRITE(numout,*)
107	#if defined key_vppm
108	IF(lwp) WRITE(numout,*) 'tra_adv_tvd : TVD_vPPM advection scheme on ', cdtype
109	#else
110	IF(lwp) WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype
111	#endif
112	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
113	!
114	l_trd = .FALSE.
115	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
116	ENDIF
117	!
118	IF( l_trd ) THEN
119	CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
120	ztrdx(:,:,:) = 0.e0 ; ztrdy(:,:,:) = 0.e0 ; ztrdz(:,:,:) = 0.e0
121	ENDIF
122	!
123	zwi(:,:,:) = 0.e0
124	!
125	! ! ===========
126	DO jn = 1, kjpt ! tracer loop
127	! ! ===========
128	! 1. Bottom value : flux set to zero
129	! ----------------------------------
130	zwx(:,:,jpk) = 0.e0 ; zwz(:,:,jpk) = 0.e0
131	zwy(:,:,jpk) = 0.e0 ; zwi(:,:,jpk) = 0.e0
132
133	! 2. upstream advection with initial mass fluxes & intermediate update
134	! --------------------------------------------------------------------
135	! upstream tracer flux in the i and j direction
136	DO jk = 1, jpkm1
137	DO jj = 1, jpjm1
138	DO ji = 1, fs_jpim1 ! vector opt.
139	! upstream scheme
140	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
141	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
142	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
143	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
144	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) )
145	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) )
146	END DO
147	END DO
148	END DO
149
150	! upstream tracer flux in the k direction
151	! Surface value
152	IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! volume variable
153	ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface
154	ENDIF
155	! Interior value
156	DO jk = 2, jpkm1
157	DO jj = 1, jpj
158	DO ji = 1, jpi
159	zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
160	zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
161	zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )
162	END DO
163	END DO
164	END DO
165
166	! total advective trend
167	DO jk = 1, jpkm1
168	z2dtt = p2dt(jk)
169	DO jj = 2, jpjm1
170	DO ji = fs_2, fs_jpim1 ! vector opt.
171	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
172	#if defined key_vppm
173	hflux(ji,jj,jk) = - (zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
174	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ))
175
176	ztra = zbtr * ( hflux(ji,jj,jk) - ( zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) )
177
178	hflux(ji,jj,jk) = hflux(ji,jj,jk) * tmask(ji,jj,jk)
179	#else
180	! total intermediate advective trends
181	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
182	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
183	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
184	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
185	#endif
186
187	! update and guess with monotonic sheme
188	zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)
189	END DO
190	END DO
191	END DO
192	! ! Lateral boundary conditions on zwi (unchanged sign)
193	CALL lbc_lnk( zwi, 'T', 1. )
194	!#if defined key_vppm
195	! CALL lbc_lnk( hflux, 'T', 1. ) ! this call seems unnecessary, H.Liu
196	!#endif
197
198	! ! trend diagnostics (contribution of upstream fluxes)
199	IF( l_trd ) THEN
200	! store intermediate advective trends
201	ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)
202	END IF
203	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
204	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
205	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) )
206	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) )
207	ENDIF
208
209	! 3. antidiffusive flux : high order minus low order
210	! --------------------------------------------------
211	! antidiffusive flux on i and j
212	DO jk = 1, jpkm1
213	DO jj = 1, jpjm1
214	DO ji = 1, fs_jpim1 ! vector opt.
215	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)
216	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)
217	END DO
218	END DO
219	END DO
220
221	! antidiffusive flux on k
222	zwz(:,:,1) = 0.e0 ! Surface value
223	!
224	DO jk = 2, jpkm1 ! Interior value
225	DO jj = 1, jpj
226	DO ji = 1, jpi
227	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)
228	END DO
229	END DO
230	END DO
231	CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions
232	CALL lbc_lnk( zwz, 'W', 1. )
233
234	! 4. monotonicity algorithm
235	! -------------------------
236	CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
237
238
239	! 5. final trend with corrected fluxes
240	! ------------------------------------
241	DO jk = 1, jpkm1
242	DO jj = 2, jpjm1
243	DO ji = fs_2, fs_jpim1 ! vector opt.
244	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
245	#if defined key_vppm
246	ztra = zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) + &
247	& zwy(ji,jj,jk) - zwy(ji ,jj-1,jk )
248	hflux(ji,jj,jk) = hflux(ji,jj,jk) - ztra * tmask(ji,jj,jk)
249	#else
250	! total advective trends
251	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
252	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
253	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
254	! add them to the general tracer trends
255	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
256	#endif
257	END DO
258	END DO
259	END DO
260	#if defined key_vppm
261	!CALL lbc_lnk( hflux, 'T', 1. ) ! This call seems unnecessary. H.Liu
262	CALL tra_adv_vppm(pta(:,:,:,jn), ptb(:,:,:,jn), hflux, z2dtt) ! pta has been updated during this call
263	#endif
264
265
266	! ! trend diagnostics (contribution of upstream fluxes)
267	IF( l_trd ) THEN
268	ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed
269	ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed
270	ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed
271
272	CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, ztrdx, pun, ptn(:,:,:,jn) )
273	CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, ztrdy, pvn, ptn(:,:,:,jn) )
274	CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, ztrdz, pwn, ptn(:,:,:,jn) )
275	END IF
276	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
277	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
278	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:)
279	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:)
280	ENDIF
281	!
282	END DO
283	!
284	CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz )
285	#if defined key_vppm
286	CALL wrk_dealloc( jpi, jpj, jpk, hflux )
287	#endif
288	IF( l_trd ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
289	!
290	IF( nn_timing == 1 ) CALL timing_stop('tra_adv_tvd')
291	!
292	END SUBROUTINE tra_adv_tvd
293
294
295	SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
296	!!---------------------------------------------------------------------
297	!! * ROUTINE nonosc *
298	!!
299	!! ** Purpose : compute monotonic tracer fluxes from the upstream
300	!! scheme and the before field by a nonoscillatory algorithm
301	!!
302	!! ** Method : ... ???
303	!! warning : pbef and paft must be masked, but the boundaries
304	!! conditions on the fluxes are not necessary zalezak (1979)
305	!! drange (1995) multi-dimensional forward-in-time and upstream-
306	!! in-space based differencing for fluid
307	!!----------------------------------------------------------------------
308	!
309	!!----------------------------------------------------------------------
310	REAL(wp), DIMENSION(jpk) , INTENT(in ) :: p2dt ! vertical profile of tracer time-step
311	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field
312	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions
313	!
314	INTEGER :: ji, jj, jk ! dummy loop indices
315	INTEGER :: ikm1 ! local integer
316	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt ! local scalars
317	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -
318	REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo
319	!!----------------------------------------------------------------------
320	!
321	IF( nn_timing == 1 ) CALL timing_start('nonosc')
322	!
323	CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
324	!
325
326	zbig = 1.e+40_wp
327	zrtrn = 1.e-15_wp
328	zbetup(:,:,jpk) = 0._wp ; zbetdo(:,:,jpk) = 0._wp
329
330
331	! Search local extrema
332	! --------------------
333	! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
334	zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ), &
335	& paft * tmask - zbig * ( 1.e0 - tmask ) )
336	zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ), &
337	& paft * tmask + zbig * ( 1.e0 - tmask ) )
338
339	DO jk = 1, jpkm1
340	ikm1 = MAX(jk-1,1)
341	z2dtt = p2dt(jk)
342	DO jj = 2, jpjm1
343	DO ji = fs_2, fs_jpim1 ! vector opt.
344
345	! search maximum in neighbourhood
346	zup = MAX( zbup(ji ,jj ,jk ), &
347	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
348	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), &
349	& zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) )
350
351	! search minimum in neighbourhood
352	zdo = MIN( zbdo(ji ,jj ,jk ), &
353	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
354	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), &
355	& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
356
357	! positive part of the flux
358	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
359	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
360	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
361
362	! negative part of the flux
363	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
364	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
365	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
366
367	! up & down beta terms
368	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
369	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
370	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
371	END DO
372	END DO
373	END DO
374	CALL lbc_lnk( zbetup, 'T', 1. ) ; CALL lbc_lnk( zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
375
376	! 3. monotonic flux in the i & j direction (paa & pbb)
377	! ----------------------------------------
378	DO jk = 1, jpkm1
379	DO jj = 2, jpjm1
380	DO ji = fs_2, fs_jpim1 ! vector opt.
381	zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
382	zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
383	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
384	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
385
386	zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
387	zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
388	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
389	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv )
390
391	! monotonic flux in the k direction, i.e. pcc
392	! -------------------------------------------
393	za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
394	zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
395	zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
396	pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1.e0 - zc) * zb )
397	END DO
398	END DO
399	END DO
400	CALL lbc_lnk( paa, 'U', -1. ) ; CALL lbc_lnk( pbb, 'V', -1. ) ! lateral boundary condition (changed sign)
401	!
402	CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
403	!
404	IF( nn_timing == 1 ) CALL timing_stop('nonosc')
405	!
406	END SUBROUTINE nonosc
407
408	!!======================================================================
409	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/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90 @ 6736

Download in other formats: