Context Navigation

advtrac.F @ 1098

Last change on this file since 1098 was 960, checked in by lsce, 17 years ago

Ajoute du parametre config_inca dans conf_gcm.F config_inca='none'(sans INCA, par defaut) config_inca='chem'(avec INCA config chemie) config_inca='aero'(avec INCA config aerosol)
Menage parmis les cles CPP INCA
Enleve le calcul d'omega dans calfis.F et active le calcul correspondant dans advtrac.F(avant uniquement pour INCA).

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

File size: 11.7 KB

Line
1	!
2	! $Header$
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	IMPLICIT NONE
18	c
19	#include "dimensions.h"
20	#include "paramet.h"
21	#include "comconst.h"
22	#include "comvert.h"
23	#include "comdissip.h"
24	#include "comgeom2.h"
25	#include "logic.h"
26	#include "temps.h"
27	#include "control.h"
28	#include "ener.h"
29	#include "description.h"
30	#include "advtrac.h"
31
32	c-------------------------------------------------------------------
33	c Arguments
34	c-------------------------------------------------------------------
35	c Ajout PPM
36	c--------------------------------------------------------
37	REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
38	c--------------------------------------------------------
39	INTEGER iapptrac
40	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
41	REAL q(ip1jmp1,llm,nqmx),masse(ip1jmp1,llm)
42	REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
43	REAL pk(ip1jmp1,llm)
44	REAL flxw(ip1jmp1,llm)
45
46	c-------------------------------------------------------------
47	c Variables locales
48	c-------------------------------------------------------------
49
50	REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
51	REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
52	REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
53	REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
54	real cpuadv(nqmx)
55	common/cpuadv/cpuadv
56
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,nqmx)
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
79	IF(iadvtr.EQ.0) THEN
80	CALL initial0(ijp1llm,pbaruc)
81	CALL initial0(ijmllm,pbarvc)
82	ENDIF
83
84	c 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	c 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	ccc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
99	c
100	ENDIF
101
102	iadvtr = iadvtr+1
103	iapptrac = iadvtr
104
105
106	c Test pour savoir si on advecte a ce pas de temps
107	IF ( iadvtr.EQ.iapp_tracvl ) THEN
108
109	cc .. Modif P.Le Van ( 20/12/97 ) ....
110	cc
111
112	c traitement des flux de masse avant advection.
113	c 1. calcul de w
114	c 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 / FLOAT(iapp_tracvl)
120
121	c 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	s - pbarvg(ij-iip1,l) + pbarvg(ij,l)
126	s + 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	s ' MAX:', zdpmax
139	ENDIF
140
141	ENDDO
142
143	c-------------------------------------------------------------------
144	c Advection proprement dite (Modification Le Croller (07/2001)
145	c-------------------------------------------------------------------
146
147	c----------------------------------------------------
148	c Calcul des moyennes basées sur la masse
149	c----------------------------------------------------
150	call massbar(massem,massebx,masseby)
151
152	c-----------------------------------------------------------
153	c Appel des sous programmes d'advection
154	c-----------------------------------------------------------
155	do iq=1,nqmx
156	c call clock(t_initial)
157	if(iadv(iq) == 0) cycle
158	c ----------------------------------------------------------------
159	c Schema de Van Leer I MUSCL
160	c ----------------------------------------------------------------
161	if(iadv(iq).eq.10) THEN
162	call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
163
164
165	c ----------------------------------------------------------------
166	c Schema "pseudo amont" + test sur humidite specifique
167	C pour la vapeur d'eau. F. Codron
168	c ----------------------------------------------------------------
169	else if(iadv(iq).eq.14) then
170	c
171	CALL vlspltqs( q(1,1,1), 2., massem, wg ,
172	* pbarug,pbarvg,dtvr,p,pk,teta )
173	c ----------------------------------------------------------------
174	c Schema de Frederic Hourdin
175	c ----------------------------------------------------------------
176	else if(iadv(iq).eq.12) then
177	c Pas de temps adaptatif
178	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
179	if (n.GT.1) then
180	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
181	s dtvr,'n=',n
182	endif
183	do indice=1,n
184	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
185	end do
186	else if(iadv(iq).eq.13) then
187	c Pas de temps adaptatif
188	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
189	if (n.GT.1) then
190	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
191	s dtvr,'n=',n
192	endif
193	do indice=1,n
194	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
195	end do
196	c ----------------------------------------------------------------
197	c Schema de pente SLOPES
198	c ----------------------------------------------------------------
199	else if (iadv(iq).eq.20) then
200	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
201
202	c ----------------------------------------------------------------
203	c Schema de Prather
204	c ----------------------------------------------------------------
205	else if (iadv(iq).eq.30) then
206	c Pas de temps adaptatif
207	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
208	if (n.GT.1) then
209	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
210	s dtvr,'n=',n
211	endif
212	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
213	s n,dtbon)
214
215	c ----------------------------------------------------------------
216	c Schemas PPM Lin et Rood
217	c ----------------------------------------------------------------
218	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
219	s iadv(iq).LE.18)) then
220
221	c Test sur le flux horizontal
222	c Pas de temps adaptatif
223	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
224	if (n.GT.1) then
225	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
226	s dtvr,'n=',n
227	endif
228	c Test sur le flux vertical
229	CFLmaxz=0.
230	do l=2,llm
231	do ij=iip2,ip1jm
232	aaa=wg(ij,l)*dtvr/massem(ij,l)
233	CFLmaxz=max(CFLmaxz,aaa)
234	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
235	CFLmaxz=max(CFLmaxz,bbb)
236	enddo
237	enddo
238	if (CFLmaxz.GE.1) then
239	write(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
240	endif
241
242	c-----------------------------------------------------------
243	c Ss-prg interface LMDZ.4->PPM3d
244	c-----------------------------------------------------------
245
246	call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
247	s apppm,bpppm,massebx,masseby,pbarug,pbarvg,
248	s unatppm,vnatppm,psppm)
249
250	do indice=1,n
251	c---------------------------------------------------------------------
252	c VL (version PPM) horiz. et PPM vert.
253	c---------------------------------------------------------------------
254	if (iadv(iq).eq.11) then
255	c Ss-prg PPM3d de Lin
256	call ppm3d(1,qppm(1,1,iq),
257	s psppm,psppm,
258	s unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
259	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
260	s fill,dum,220.)
261
262	c----------------------------------------------------------------------
263	c Monotonic PPM
264	c----------------------------------------------------------------------
265	else if (iadv(iq).eq.16) then
266	c Ss-prg PPM3d de Lin
267	call ppm3d(1,qppm(1,1,iq),
268	s psppm,psppm,
269	s unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
270	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
271	s fill,dum,220.)
272	c---------------------------------------------------------------------
273
274	c---------------------------------------------------------------------
275	c Semi Monotonic PPM
276	c---------------------------------------------------------------------
277	else if (iadv(iq).eq.17) then
278	c Ss-prg PPM3d de Lin
279	call ppm3d(1,qppm(1,1,iq),
280	s psppm,psppm,
281	s unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
282	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
283	s fill,dum,220.)
284	c---------------------------------------------------------------------
285
286	c---------------------------------------------------------------------
287	c Positive Definite PPM
288	c---------------------------------------------------------------------
289	else if (iadv(iq).eq.18) then
290	c Ss-prg PPM3d de Lin
291	call ppm3d(1,qppm(1,1,iq),
292	s psppm,psppm,
293	s unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
294	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
295	s fill,dum,220.)
296	c---------------------------------------------------------------------
297	endif
298	enddo
299	c-----------------------------------------------------------------
300	c Ss-prg interface PPM3d-LMDZ.4
301	c-----------------------------------------------------------------
302	call interpost(q(1,1,iq),qppm(1,1,iq))
303	endif
304	c----------------------------------------------------------------------
305
306	c-----------------------------------------------------------------
307	c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
308	c et Nord j=1
309	c-----------------------------------------------------------------
310
311	c call traceurpole(q(1,1,iq),massem)
312
313	c calcul du temps cpu pour un schema donne
314
315	c call clock(t_final)
316	cym tps_cpu=t_final-t_initial
317	cym cpuadv(iq)=cpuadv(iq)+tps_cpu
318
319	end DO
320
321
322	c------------------------------------------------------------------
323	c on reinitialise a zero les flux de masse cumules
324	c---------------------------------------------------
325	iadvtr=0
326
327	ENDIF ! if iadvtr.EQ.iapp_tracvl
328
329	RETURN
330	END
331

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

Download in other formats:

Original Format