Context Navigation

advtrac.F90 @ 5305

Last change on this file since 5305 was 2622, checked in by Ehouarn Millour, 8 years ago
Some code tidying: turn ener.h into ener_mod.F90 EM
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.7 KB

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

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

Download in other formats:

Original Format