Context Navigation

advtrac.F90 @ 4100

Last change on this file since 4100 was 1999, checked in by Laurent Fairhead, 11 years ago
Merged trunk changes r1920:1997 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 14.3 KB

Line
1	! $Id: advtrac.F90 1999 2014-03-20 09:57:19Z oboucher $
2
3	SUBROUTINE advtrac(pbaru,pbarv , p, masse,q,iapptrac,teta, flxw, pk)
4	! Auteur : F. Hourdin
5	!
6	! Modif. P. Le Van (20/12/97)
7	! F. Codron (10/99)
8	! D. Le Croller (07/2001)
9	! M.A Filiberti (04/2002)
10	!
11	USE infotrac, ONLY: nqtot, iadv
12	USE control_mod, ONLY: iapp_tracvl, day_step
13
14
15	IMPLICIT NONE
16	!
17	include "dimensions.h"
18	include "paramet.h"
19	include "comconst.h"
20	include "comvert.h"
21	include "comdissip.h"
22	include "comgeom2.h"
23	include "logic.h"
24	include "temps.h"
25	include "ener.h"
26	include "description.h"
27	include "iniprint.h"
28
29	!-------------------------------------------------------------------
30	! Arguments
31	!-------------------------------------------------------------------
32	INTEGER,INTENT(OUT) :: iapptrac
33	REAL,INTENT(IN) :: pbaru(ip1jmp1,llm)
34	REAL,INTENT(IN) :: pbarv(ip1jm,llm)
35	REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
36	REAL,INTENT(IN) :: masse(ip1jmp1,llm)
37	REAL,INTENT(IN) :: p( ip1jmp1,llmp1 )
38	REAL,INTENT(IN) :: teta(ip1jmp1,llm)
39	REAL,INTENT(IN) :: pk(ip1jmp1,llm)
40	REAL,INTENT(OUT) :: flxw(ip1jmp1,llm)
41	!-------------------------------------------------------------------
42	! Ajout PPM
43	!--------------------------------------------------------
44	REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
45	!-------------------------------------------------------------
46	! Variables locales
47	!-------------------------------------------------------------
48
49	REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
50	REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
51	REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
52	REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
53	INTEGER iadvtr
54	INTEGER ij,l,iq,iiq
55	REAL zdpmin, zdpmax
56	EXTERNAL minmax
57	SAVE iadvtr, massem, pbaruc, pbarvc
58	DATA iadvtr/0/
59	!----------------------------------------------------------
60	! Rajouts pour PPM
61	!----------------------------------------------------------
62	INTEGER indice,n
63	REAL dtbon ! Pas de temps adaptatif pour que CFL<1
64	REAL CFLmaxz,aaa,bbb ! CFL maximum
65	REAL psppm(iim,jjp1) ! pression au sol
66	REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
67	REAL qppm(iim*jjp1,llm,nqtot)
68	REAL fluxwppm(iim,jjp1,llm)
69	REAL apppm(llmp1), bpppm(llmp1)
70	LOGICAL dum,fill
71	DATA fill/.true./
72	DATA dum/.true./
73
74	integer,save :: countcfl=0
75	real cflx(ip1jmp1,llm)
76	real cfly(ip1jm,llm)
77	real cflz(ip1jmp1,llm)
78	real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
79
80	IF(iadvtr.EQ.0) THEN
81	CALL initial0(ijp1llm,pbaruc)
82	CALL initial0(ijmllm,pbarvc)
83	ENDIF
84
85	! accumulation des flux de masse horizontaux
86	DO l=1,llm
87	DO ij = 1,ip1jmp1
88	pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
89	ENDDO
90	DO ij = 1,ip1jm
91	pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
92	ENDDO
93	ENDDO
94
95	! selection de la masse instantannee des mailles avant le transport.
96	IF(iadvtr.EQ.0) THEN
97
98	CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
99	!cc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
100	!
101	ENDIF
102
103	iadvtr = iadvtr+1
104	iapptrac = iadvtr
105
106
107	! Test pour savoir si on advecte a ce pas de temps
108	IF ( iadvtr.EQ.iapp_tracvl ) THEN
109
110	!c .. Modif P.Le Van ( 20/12/97 ) ....
111	!c
112
113	! traitement des flux de masse avant advection.
114	! 1. calcul de w
115	! 2. groupement des mailles pres du pole.
116
117	CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
118
119	! ... Flux de masse diaganostiques traceurs
120	flxw = wg / REAL(iapp_tracvl)
121
122	! test sur l'eventuelle creation de valeurs negatives de la masse
123	DO l=1,llm-1
124	DO ij = iip2+1,ip1jm
125	zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) &
126	- pbarvg(ij-iip1,l) + pbarvg(ij,l) &
127	+ wg(ij,l+1) - wg(ij,l)
128	ENDDO
129	CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
130	DO ij = iip2,ip1jm
131	zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
132	ENDDO
133
134
135	CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
136
137	IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
138	PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin, &
139	' MAX:', zdpmax
140	ENDIF
141
142	ENDDO
143
144
145	!-------------------------------------------------------------------
146	! Calcul des criteres CFL en X, Y et Z
147	!-------------------------------------------------------------------
148
149	if (countcfl == 0. ) then
150	cflxmax(:)=0.
151	cflymax(:)=0.
152	cflzmax(:)=0.
153	endif
154
155	countcfl=countcfl+iapp_tracvl
156	cflx(:,:)=0.
157	cfly(:,:)=0.
158	cflz(:,:)=0.
159	do l=1,llm
160	do ij=iip2,ip1jm-1
161	if (pbarug(ij,l)>=0.) then
162	cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
163	else
164	cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
165	endif
166	enddo
167	enddo
168	do l=1,llm
169	do ij=iip2,ip1jm-1,iip1
170	cflx(ij+iip1,l)=cflx(ij,l)
171	enddo
172	enddo
173
174	do l=1,llm
175	do ij=1,ip1jm
176	if (pbarvg(ij,l)>=0.) then
177	cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
178	else
179	cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
180	endif
181	enddo
182	enddo
183
184	do l=2,llm
185	do ij=1,ip1jm
186	if (wg(ij,l)>=0.) then
187	cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
188	else
189	cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
190	endif
191	enddo
192	enddo
193
194	do l=1,llm
195	cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
196	cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
197	cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
198	enddo
199
200	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
201	! Par defaut, on sort le diagnostic des CFL tous les jours.
202	! Si on veut le sortir a chaque pas d'advection en cas de plantage
203	! if (countcfl==iapp_tracvl) then
204	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
205	if (countcfl==day_step) then
206	do l=1,llm
207	write(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), &
208	cflzmax(l)
209	enddo
210	countcfl=0
211	endif
212
213	!-------------------------------------------------------------------
214	! Advection proprement dite (Modification Le Croller (07/2001)
215	!-------------------------------------------------------------------
216
217	!----------------------------------------------------
218	! Calcul des moyennes basées sur la masse
219	!----------------------------------------------------
220	call massbar(massem,massebx,masseby)
221
222	!-----------------------------------------------------------
223	! Appel des sous programmes d'advection
224	!-----------------------------------------------------------
225	do iq=1,nqtot
226	! call clock(t_initial)
227	if(iadv(iq) == 0) cycle
228	! ----------------------------------------------------------------
229	! Schema de Van Leer I MUSCL
230	! ----------------------------------------------------------------
231	if(iadv(iq).eq.10) THEN
232	call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
233
234	! ----------------------------------------------------------------
235	! Schema "pseudo amont" + test sur humidite specifique
236	! pour la vapeur d'eau. F. Codron
237	! ----------------------------------------------------------------
238	else if(iadv(iq).eq.14) then
239	!
240	CALL vlspltqs( q(1,1,1), 2., massem, wg , &
241	pbarug,pbarvg,dtvr,p,pk,teta )
242	! ----------------------------------------------------------------
243	! Schema de Frederic Hourdin
244	! ----------------------------------------------------------------
245	else if(iadv(iq).eq.12) then
246	! Pas de temps adaptatif
247	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
248	if (n.GT.1) then
249	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
250	dtvr,'n=',n
251	endif
252	do indice=1,n
253	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
254	end do
255	else if(iadv(iq).eq.13) then
256	! Pas de temps adaptatif
257	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
258	if (n.GT.1) then
259	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
260	dtvr,'n=',n
261	endif
262	do indice=1,n
263	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
264	end do
265	! ----------------------------------------------------------------
266	! Schema de pente SLOPES
267	! ----------------------------------------------------------------
268	else if (iadv(iq).eq.20) then
269	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
270
271	! ----------------------------------------------------------------
272	! Schema de Prather
273	! ----------------------------------------------------------------
274	else if (iadv(iq).eq.30) then
275	! Pas de temps adaptatif
276	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
277	if (n.GT.1) then
278	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
279	dtvr,'n=',n
280	endif
281	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
282	n,dtbon)
283
284	! ----------------------------------------------------------------
285	! Schemas PPM Lin et Rood
286	! ----------------------------------------------------------------
287	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
288	iadv(iq).LE.18)) then
289
290	! Test sur le flux horizontal
291	! Pas de temps adaptatif
292	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
293	if (n.GT.1) then
294	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
295	dtvr,'n=',n
296	endif
297	! Test sur le flux vertical
298	CFLmaxz=0.
299	do l=2,llm
300	do ij=iip2,ip1jm
301	aaa=wg(ij,l)*dtvr/massem(ij,l)
302	CFLmaxz=max(CFLmaxz,aaa)
303	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
304	CFLmaxz=max(CFLmaxz,bbb)
305	enddo
306	enddo
307	if (CFLmaxz.GE.1) then
308	write(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
309	endif
310
311	!-----------------------------------------------------------
312	! Ss-prg interface LMDZ.4->PPM3d
313	!-----------------------------------------------------------
314
315	call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
316	apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
317	unatppm,vnatppm,psppm)
318
319	do indice=1,n
320	!----------------------------------------------------------------
321	! VL (version PPM) horiz. et PPM vert.
322	!----------------------------------------------------------------
323	if (iadv(iq).eq.11) then
324	! Ss-prg PPM3d de Lin
325	call ppm3d(1,qppm(1,1,iq), &
326	psppm,psppm, &
327	unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
328	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
329	fill,dum,220.)
330
331	!-------------------------------------------------------------
332	! Monotonic PPM
333	!-------------------------------------------------------------
334	else if (iadv(iq).eq.16) then
335	! Ss-prg PPM3d de Lin
336	call ppm3d(1,qppm(1,1,iq), &
337	psppm,psppm, &
338	unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
339	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
340	fill,dum,220.)
341	!-------------------------------------------------------------
342
343	!-------------------------------------------------------------
344	! Semi Monotonic PPM
345	!-------------------------------------------------------------
346	else if (iadv(iq).eq.17) then
347	! Ss-prg PPM3d de Lin
348	call ppm3d(1,qppm(1,1,iq), &
349	psppm,psppm, &
350	unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
351	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
352	fill,dum,220.)
353	!-------------------------------------------------------------
354
355	!-------------------------------------------------------------
356	! Positive Definite PPM
357	!-------------------------------------------------------------
358	else if (iadv(iq).eq.18) then
359	! Ss-prg PPM3d de Lin
360	call ppm3d(1,qppm(1,1,iq), &
361	psppm,psppm, &
362	unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
363	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
364	fill,dum,220.)
365	!-------------------------------------------------------------
366	endif
367	enddo
368	!-----------------------------------------------------------------
369	! Ss-prg interface PPM3d-LMDZ.4
370	!-----------------------------------------------------------------
371	call interpost(q(1,1,iq),qppm(1,1,iq))
372	endif
373	!----------------------------------------------------------------------
374
375	!-----------------------------------------------------------------
376	! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
377	! et Nord j=1
378	!-----------------------------------------------------------------
379
380	! call traceurpole(q(1,1,iq),massem)
381
382	! calcul du temps cpu pour un schema donne
383
384	! call clock(t_final)
385	!ym tps_cpu=t_final-t_initial
386	!ym cpuadv(iq)=cpuadv(iq)+tps_cpu
387
388	end DO
389
390
391	!------------------------------------------------------------------
392	! on reinitialise a zero les flux de masse cumules
393	!---------------------------------------------------
394	iadvtr=0
395
396	ENDIF ! if iadvtr.EQ.iapp_tracvl
397
398	END SUBROUTINE advtrac

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

Download in other formats:

Original Format