Context Navigation

advtrac_loc.F @ 5360

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

fix in dyn3d the array out of bound issue (in qminimum) that was only corrected in dyn3dmem in rev 2285.
comment out many invasive debug writes.

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

File size: 13.0 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, ONLY: nqtot, iadv, ok_iso_verif
27	USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
28	USE advtrac_mod, ONLY: finmasse
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),SAVE :: testRequest
82	!$OMP THREADPRIVATE(testRequest)
83
84	c test sur l''eventuelle creation de valeurs negatives de la masse
85	ijb=ij_begin
86	ije=ij_end
87	if (pole_nord) ijb=ij_begin+iip1
88	if (pole_sud) ije=ij_end-iip1
89
90	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
91	DO l=1,llm-1
92	DO ij = ijb+1,ije
93	zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l)
94	s - pbarvg(ij-iip1,l) + pbarvg(ij,l)
95	s + wg(ij,l+1) - wg(ij,l)
96	ENDDO
97
98	c CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
99	c ym ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
100
101	do ij=ijb,ije-iip1+1,iip1
102	zdp(ij)=zdp(ij+iip1-1)
103	enddo
104
105	DO ij = ijb,ije
106	zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
107	ENDDO
108
109
110	c CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
111	c ym ---> eventuellement a revoir
112	CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
113
114	IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
115	PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin,
116	s ' MAX:', zdpmax
117	ENDIF
118
119	ENDDO
120	c$OMP END DO NOWAIT
121
122	c-------------------------------------------------------------------
123	c Advection proprement dite (Modification Le Croller (07/2001)
124	c-------------------------------------------------------------------
125
126	c----------------------------------------------------
127	c Calcul des moyennes basées sur la masse
128	c----------------------------------------------------
129
130	cym ----> Normalement, inutile pour les schémas classiques
131	cym ----> Revérifier lors de la parallélisation des autres schemas
132
133	cym call massbar_p(massem,massebx,masseby)
134
135	#ifdef DEBUG_IO
136	CALL WriteField_u('massem',massem)
137	CALL WriteField_u('wg',wg)
138	CALL WriteField_u('pbarug',pbarug)
139	CALL WriteField_v('pbarvg',pbarvg)
140	CALL WriteField_u('p_tmp',p)
141	CALL WriteField_u('pk_tmp',pk)
142	CALL WriteField_u('teta_tmp',teta)
143	do j=1,nqtot
144	call WriteField_u('q_adv'//trim(int2str(j)),
145	. q(:,:,j))
146	enddo
147	#endif
148
149	!
150	! call Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
151	!
152	! call SendRequest(TestRequest)
153	!c$OMP BARRIER
154	! call WaitRequest(TestRequest)
155	c$OMP BARRIER
156
157	!write(,) 'advtrac 157: appel de vlspltgen_loc'
158	call vlspltgen_loc( q,iadv, 2., massem, wg ,
159	* pbarug,pbarvg,dtvr,p,
160	* pk,teta )
161
162	!write(,) 'advtrac 162: apres appel vlspltgen_loc'
163	if (ok_iso_verif) then
164	call check_isotopes(q,ijb_u,ije_u,'advtrac 162')
165	endif !if (ok_iso_verif) then
166
167	#ifdef DEBUG_IO
168	do j=1,nqtot
169	call WriteField_u('q_adv'//trim(int2str(j)),
170	. q(:,:,j))
171	enddo
172	#endif
173
174	GOTO 1234
175	c-----------------------------------------------------------
176	c Appel des sous programmes d'advection
177	c-----------------------------------------------------------
178	do iq=1,nqtot
179	c call clock(t_initial)
180	if(iadv(iq) == 0) cycle
181	c ----------------------------------------------------------------
182	c Schema de Van Leer I MUSCL
183	c ----------------------------------------------------------------
184	if(iadv(iq).eq.10) THEN
185
186	!LF call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
187
188	c ----------------------------------------------------------------
189	c Schema "pseudo amont" + test sur humidite specifique
190	C pour la vapeur d'eau. F. Codron
191	c ----------------------------------------------------------------
192	else if(iadv(iq).eq.14) then
193	c
194	cym stop 'advtrac : appel à vlspltqs :schema non parallelise'
195	!LF CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
196	!LF * pbarug,pbarvg,dtvr,p,pk,teta )
197	c ----------------------------------------------------------------
198	c Schema de Frederic Hourdin
199	c ----------------------------------------------------------------
200	else if(iadv(iq).eq.12) then
201	stop 'advtrac : schema non parallelise'
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,1)
210	end do
211	else if(iadv(iq).eq.13) then
212	stop 'advtrac : schema non parallelise'
213	c Pas de temps adaptatif
214	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
215	if (n.GT.1) then
216	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
217	s dtvr,'n=',n
218	endif
219	do indice=1,n
220	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
221	end do
222	c ----------------------------------------------------------------
223	c Schema de pente SLOPES
224	c ----------------------------------------------------------------
225	else if (iadv(iq).eq.20) then
226	stop 'advtrac : schema non parallelise'
227
228	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
229
230	c ----------------------------------------------------------------
231	c Schema de Prather
232	c ----------------------------------------------------------------
233	else if (iadv(iq).eq.30) then
234	stop 'advtrac : schema non parallelise'
235	c Pas de temps adaptatif
236	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
237	if (n.GT.1) then
238	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
239	s dtvr,'n=',n
240	endif
241	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
242	s n,dtbon)
243	c ----------------------------------------------------------------
244	c Schemas PPM Lin et Rood
245	c ----------------------------------------------------------------
246	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
247	s iadv(iq).LE.18)) then
248
249	stop 'advtrac : schema non parallelise'
250
251	c Test sur le flux horizontal
252	c Pas de temps adaptatif
253	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
254	if (n.GT.1) then
255	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
256	s dtvr,'n=',n
257	endif
258	c Test sur le flux vertical
259	CFLmaxz=0.
260	do l=2,llm
261	do ij=iip2,ip1jm
262	aaa=wg(ij,l)*dtvr/massem(ij,l)
263	CFLmaxz=max(CFLmaxz,aaa)
264	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
265	CFLmaxz=max(CFLmaxz,bbb)
266	enddo
267	enddo
268	if (CFLmaxz.GE.1) then
269	write(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
270	endif
271
272	c-----------------------------------------------------------
273	c Ss-prg interface LMDZ.4->PPM3d
274	c-----------------------------------------------------------
275
276	call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
277	s apppm,bpppm,massebx,masseby,pbarug,pbarvg,
278	s unatppm,vnatppm,psppm)
279
280	do indice=1,n
281	c---------------------------------------------------------------------
282	c VL (version PPM) horiz. et PPM vert.
283	c---------------------------------------------------------------------
284	if (iadv(iq).eq.11) then
285	c Ss-prg PPM3d de Lin
286	call ppm3d(1,qppm(1,1,iq),
287	s psppm,psppm,
288	s unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
289	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
290	s fill,dum,220.)
291
292	c----------------------------------------------------------------------
293	c Monotonic PPM
294	c----------------------------------------------------------------------
295	else if (iadv(iq).eq.16) then
296	c Ss-prg PPM3d de Lin
297	call ppm3d(1,qppm(1,1,iq),
298	s psppm,psppm,
299	s unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
300	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
301	s fill,dum,220.)
302	c---------------------------------------------------------------------
303
304	c---------------------------------------------------------------------
305	c Semi Monotonic PPM
306	c---------------------------------------------------------------------
307	else if (iadv(iq).eq.17) then
308	c Ss-prg PPM3d de Lin
309	call ppm3d(1,qppm(1,1,iq),
310	s psppm,psppm,
311	s unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
312	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
313	s fill,dum,220.)
314	c---------------------------------------------------------------------
315
316	c---------------------------------------------------------------------
317	c Positive Definite PPM
318	c---------------------------------------------------------------------
319	else if (iadv(iq).eq.18) then
320	c Ss-prg PPM3d de Lin
321	call ppm3d(1,qppm(1,1,iq),
322	s psppm,psppm,
323	s unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
324	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
325	s fill,dum,220.)
326	c---------------------------------------------------------------------
327	endif
328	enddo
329	c-----------------------------------------------------------------
330	c Ss-prg interface PPM3d-LMDZ.4
331	c-----------------------------------------------------------------
332	call interpost(q(1,1,iq),qppm(1,1,iq))
333	endif
334	c----------------------------------------------------------------------
335
336	c-----------------------------------------------------------------
337	c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
338	c et Nord j=1
339	c-----------------------------------------------------------------
340
341	c call traceurpole(q(1,1,iq),massem)
342
343	c calcul du temps cpu pour un schema donne
344
345
346	end DO
347
348	1234 CONTINUE
349	c$OMP BARRIER
350
351	if (planet_type=="earth") then
352
353	ijb=ij_begin
354	ije=ij_end
355
356	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
357	DO l = 1, llm
358	DO ij = ijb, ije
359	finmasse(ij,l) = p(ij,l) - p(ij,l+1)
360	ENDDO
361	ENDDO
362	c$OMP END DO
363
364	! CRisi: on passe nqtot et non nq
365	CALL qminimum_loc( q, nqtot, finmasse )
366
367	endif ! of if (planet_type=="earth")
368
369	RETURN
370	END
371

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

Download in other formats:

Original Format