Context Navigation

advtrac.F @ 1035

Last change on this file since 1035 was 524, checked in by lmdzadmin, 20 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 12.2 KB

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

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

Download in other formats:

Original Format