Context Navigation

advtrac.F90 @ 5360

Last change on this file since 5360 was 2603, checked in by Ehouarn Millour, 8 years ago
Cleanup in the dynamics: turn logic.h into module logic_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.8 KB

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

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

Download in other formats:

Original Format