Context Navigation

advtrac.F @ 5360

Last change on this file since 5360 was 1446, checked in by Ehouarn Millour, 14 years ago

Implemented modifications to enable running with only one tracer for planet types different from "earth". Rem: If flag 'planet_type' is set to "earth" (default behaviour) then there must be at least 2 tracers for the dynamics to function properly.

These updates do not induce any changes in model outputs with respect to previous revisions.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 13.8 KB

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

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

Download in other formats:

Original Format