Context Navigation

advtrac.F90 @ 2598

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

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

Download in other formats:

Original Format