Context Navigation

advtrac_loc.F @ 1823

Last change on this file since 1823 was 1823, checked in by Ehouarn Millour, 11 years ago
Remplacement de parallel.F90 (en conflit avec orchidée) par parallel_lmdz.F90. UG ......................................... Renaming parallel.F90 (conflicting with orchidée) into parallel_lmdz.F90. UG
File size: 12.6 KB

Line
1	!
2	! $Id$
3	!
4	c
5	c
6	#define DEBUG_IO
7	#undef DEBUG_IO
8	SUBROUTINE advtrac_loc(pbarug,pbarvg ,wg,
9	* p, massem,q,teta,
10	* pk )
11
12	c Auteur : F. Hourdin
13	c
14	c Modif. P. Le Van (20/12/97)
15	c F. Codron (10/99)
16	c D. Le Croller (07/2001)
17	c M.A Filiberti (04/2002)
18	c
19	USE parallel_lmdz
20	USE Write_Field_loc
21	USE Write_Field
22	USE Bands
23	USE mod_hallo
24	USE Vampir
25	USE times
26	USE infotrac
27	USE control_mod
28	USE advtrac_mod
29	IMPLICIT NONE
30	c
31	#include "dimensions.h"
32	#include "paramet.h"
33	#include "comconst.h"
34	#include "comvert.h"
35	#include "comdissip.h"
36	#include "comgeom2.h"
37	#include "logic.h"
38	#include "temps.h"
39	#include "ener.h"
40	#include "description.h"
41
42	c-------------------------------------------------------------------
43	c Arguments
44	c-------------------------------------------------------------------
45	c Ajout PPM
46	c--------------------------------------------------------
47	REAL massebx(ijb_u:ije_u,llm),masseby(ijb_v:ije_v,llm)
48	c--------------------------------------------------------
49	INTEGER iapptrac
50	REAL pbarug(ijb_u:ije_u,llm),pbarvg(ijb_v:ije_v,llm)
51	REAL wg(ijb_u:ije_u,llm)
52	REAL q(ijb_u:ije_u,llm,nqtot),massem(ijb_u:ije_u,llm)
53	REAL p( ijb_u:ije_u,llmp1 ),teta(ijb_u:ije_u,llm)
54	REAL pk(ijb_u:ije_u,llm)
55
56	c-------------------------------------------------------------
57	c Variables locales
58	c-------------------------------------------------------------
59
60	REAL zdp(ijb_u:ije_u)
61	REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
62	INTEGER,SAVE :: iadvtr=0
63	c$OMP THREADPRIVATE(iadvtr)
64	INTEGER ij,l,iq,iiq
65	REAL zdpmin, zdpmax
66	c----------------------------------------------------------
67	c Rajouts pour PPM
68	c----------------------------------------------------------
69	INTEGER indice,n
70	REAL dtbon ! Pas de temps adaptatif pour que CFL<1
71	REAL CFLmaxz,aaa,bbb ! CFL maximum
72	REAL psppm(iim,jjb_u:jje_u) ! pression au sol
73	REAL unatppm(iim,jjb_u:jje_u,llm),vnatppm(iim,jjb_u:jje_u,llm)
74	REAL qppm(iim*jjnb_u,llm,nqtot)
75	REAL fluxwppm(iim,jjb_u:jje_u,llm)
76	REAL apppm(llmp1), bpppm(llmp1)
77	LOGICAL dum,fill
78	DATA fill/.true./
79	DATA dum/.true./
80	integer ijb,ije,ijbu,ijbv,ijeu,ijev,j
81	type(Request) :: testRequest
82
83	c test sur l'eventuelle creation de valeurs negatives de la masse
84	ijb=ij_begin
85	ije=ij_end
86	if (pole_nord) ijb=ij_begin+iip1
87	if (pole_sud) ije=ij_end-iip1
88
89	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
90	DO l=1,llm-1
91	DO ij = ijb+1,ije
92	zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l)
93	s - pbarvg(ij-iip1,l) + pbarvg(ij,l)
94	s + wg(ij,l+1) - wg(ij,l)
95	ENDDO
96
97	c CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
98	c ym ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
99
100	do ij=ijb,ije-iip1+1,iip1
101	zdp(ij)=zdp(ij+iip1-1)
102	enddo
103
104	DO ij = ijb,ije
105	zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
106	ENDDO
107
108
109	c CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
110	c ym ---> eventuellement a revoir
111	CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
112
113	IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
114	PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin,
115	s ' MAX:', zdpmax
116	ENDIF
117
118	ENDDO
119	c$OMP END DO NOWAIT
120
121	c-------------------------------------------------------------------
122	c Advection proprement dite (Modification Le Croller (07/2001)
123	c-------------------------------------------------------------------
124
125	c----------------------------------------------------
126	c Calcul des moyennes basées sur la masse
127	c----------------------------------------------------
128
129	cym ----> Normalement, inutile pour les schémas classiques
130	cym ----> Revérifier lors de la parallélisation des autres schemas
131
132	cym call massbar_p(massem,massebx,masseby)
133
134	#ifdef DEBUG_IO
135	CALL WriteField_u('massem',massem)
136	CALL WriteField_u('wg',wg)
137	CALL WriteField_u('pbarug',pbarug)
138	CALL WriteField_v('pbarvg',pbarvg)
139	CALL WriteField_u('p_tmp',p)
140	CALL WriteField_u('pk_tmp',pk)
141	CALL WriteField_u('teta_tmp',teta)
142	do j=1,nqtot
143	call WriteField_u('q_adv'//trim(int2str(j)),
144	. q(:,:,j))
145	enddo
146	#endif
147
148	!
149	! call Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
150	!
151	! call SendRequest(TestRequest)
152	!c$OMP BARRIER
153	! call WaitRequest(TestRequest)
154	c$OMP BARRIER
155
156	call vlspltgen_loc( q,iadv, 2., massem, wg ,
157	* pbarug,pbarvg,dtvr,p,
158	* pk,teta )
159
160	#ifdef DEBUG_IO
161	do j=1,nqtot
162	call WriteField_u('q_adv'//trim(int2str(j)),
163	. q(:,:,j))
164	enddo
165	#endif
166
167	GOTO 1234
168	c-----------------------------------------------------------
169	c Appel des sous programmes d'advection
170	c-----------------------------------------------------------
171	do iq=1,nqtot
172	c call clock(t_initial)
173	if(iadv(iq) == 0) cycle
174	c ----------------------------------------------------------------
175	c Schema de Van Leer I MUSCL
176	c ----------------------------------------------------------------
177	if(iadv(iq).eq.10) THEN
178
179	!LF call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
180
181	c ----------------------------------------------------------------
182	c Schema "pseudo amont" + test sur humidite specifique
183	C pour la vapeur d'eau. F. Codron
184	c ----------------------------------------------------------------
185	else if(iadv(iq).eq.14) then
186	c
187	cym stop 'advtrac : appel à vlspltqs :schema non parallelise'
188	!LF CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
189	!LF * pbarug,pbarvg,dtvr,p,pk,teta )
190	c ----------------------------------------------------------------
191	c Schema de Frederic Hourdin
192	c ----------------------------------------------------------------
193	else if(iadv(iq).eq.12) then
194	stop 'advtrac : schema non parallelise'
195	c Pas de temps adaptatif
196	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
197	if (n.GT.1) then
198	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
199	s dtvr,'n=',n
200	endif
201	do indice=1,n
202	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
203	end do
204	else if(iadv(iq).eq.13) then
205	stop 'advtrac : schema non parallelise'
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	do indice=1,n
213	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
214	end do
215	c ----------------------------------------------------------------
216	c Schema de pente SLOPES
217	c ----------------------------------------------------------------
218	else if (iadv(iq).eq.20) then
219	stop 'advtrac : schema non parallelise'
220
221	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
222
223	c ----------------------------------------------------------------
224	c Schema de Prather
225	c ----------------------------------------------------------------
226	else if (iadv(iq).eq.30) then
227	stop 'advtrac : schema non parallelise'
228	c Pas de temps adaptatif
229	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
230	if (n.GT.1) then
231	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
232	s dtvr,'n=',n
233	endif
234	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
235	s n,dtbon)
236	c ----------------------------------------------------------------
237	c Schemas PPM Lin et Rood
238	c ----------------------------------------------------------------
239	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
240	s iadv(iq).LE.18)) then
241
242	stop 'advtrac : schema non parallelise'
243
244	c Test sur le flux horizontal
245	c Pas de temps adaptatif
246	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
247	if (n.GT.1) then
248	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
249	s dtvr,'n=',n
250	endif
251	c Test sur le flux vertical
252	CFLmaxz=0.
253	do l=2,llm
254	do ij=iip2,ip1jm
255	aaa=wg(ij,l)*dtvr/massem(ij,l)
256	CFLmaxz=max(CFLmaxz,aaa)
257	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
258	CFLmaxz=max(CFLmaxz,bbb)
259	enddo
260	enddo
261	if (CFLmaxz.GE.1) then
262	write(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
263	endif
264
265	c-----------------------------------------------------------
266	c Ss-prg interface LMDZ.4->PPM3d
267	c-----------------------------------------------------------
268
269	call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
270	s apppm,bpppm,massebx,masseby,pbarug,pbarvg,
271	s unatppm,vnatppm,psppm)
272
273	do indice=1,n
274	c---------------------------------------------------------------------
275	c VL (version PPM) horiz. et PPM vert.
276	c---------------------------------------------------------------------
277	if (iadv(iq).eq.11) 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,2,2,2,1,
282	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
283	s fill,dum,220.)
284
285	c----------------------------------------------------------------------
286	c Monotonic PPM
287	c----------------------------------------------------------------------
288	else if (iadv(iq).eq.16) then
289	c Ss-prg PPM3d de Lin
290	call ppm3d(1,qppm(1,1,iq),
291	s psppm,psppm,
292	s unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
293	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
294	s fill,dum,220.)
295	c---------------------------------------------------------------------
296
297	c---------------------------------------------------------------------
298	c Semi Monotonic PPM
299	c---------------------------------------------------------------------
300	else if (iadv(iq).eq.17) then
301	c Ss-prg PPM3d de Lin
302	call ppm3d(1,qppm(1,1,iq),
303	s psppm,psppm,
304	s unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
305	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
306	s fill,dum,220.)
307	c---------------------------------------------------------------------
308
309	c---------------------------------------------------------------------
310	c Positive Definite PPM
311	c---------------------------------------------------------------------
312	else if (iadv(iq).eq.18) then
313	c Ss-prg PPM3d de Lin
314	call ppm3d(1,qppm(1,1,iq),
315	s psppm,psppm,
316	s unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
317	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
318	s fill,dum,220.)
319	c---------------------------------------------------------------------
320	endif
321	enddo
322	c-----------------------------------------------------------------
323	c Ss-prg interface PPM3d-LMDZ.4
324	c-----------------------------------------------------------------
325	call interpost(q(1,1,iq),qppm(1,1,iq))
326	endif
327	c----------------------------------------------------------------------
328
329	c-----------------------------------------------------------------
330	c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
331	c et Nord j=1
332	c-----------------------------------------------------------------
333
334	c call traceurpole(q(1,1,iq),massem)
335
336	c calcul du temps cpu pour un schema donne
337
338
339	end DO
340
341	1234 CONTINUE
342	c$OMP BARRIER
343
344	if (planet_type=="earth") then
345
346	ijb=ij_begin
347	ije=ij_end
348
349	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
350	DO l = 1, llm
351	DO ij = ijb, ije
352	finmasse(ij,l) = p(ij,l) - p(ij,l+1)
353	ENDDO
354	ENDDO
355	c$OMP END DO
356
357	CALL qminimum_loc( q, 2, finmasse )
358
359	endif ! of if (planet_type=="earth")
360
361	RETURN
362	END
363

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

Download in other formats:

Original Format