Context Navigation

advtrac_loc.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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: 12.7 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),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	call vlspltgen_loc( q,iadv, 2., massem, wg ,
158	* pbarug,pbarvg,dtvr,p,
159	* pk,teta )
160
161	#ifdef DEBUG_IO
162	do j=1,nqtot
163	call WriteField_u('q_adv'//trim(int2str(j)),
164	. q(:,:,j))
165	enddo
166	#endif
167
168	GOTO 1234
169	c-----------------------------------------------------------
170	c Appel des sous programmes d'advection
171	c-----------------------------------------------------------
172	do iq=1,nqtot
173	c call clock(t_initial)
174	if(iadv(iq) == 0) cycle
175	c ----------------------------------------------------------------
176	c Schema de Van Leer I MUSCL
177	c ----------------------------------------------------------------
178	if(iadv(iq).eq.10) THEN
179
180	!LF call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
181
182	c ----------------------------------------------------------------
183	c Schema "pseudo amont" + test sur humidite specifique
184	C pour la vapeur d'eau. F. Codron
185	c ----------------------------------------------------------------
186	else if(iadv(iq).eq.14) then
187	c
188	cym stop 'advtrac : appel à vlspltqs :schema non parallelise'
189	!LF CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
190	!LF * pbarug,pbarvg,dtvr,p,pk,teta )
191	c ----------------------------------------------------------------
192	c Schema de Frederic Hourdin
193	c ----------------------------------------------------------------
194	else if(iadv(iq).eq.12) then
195	stop 'advtrac : schema non parallelise'
196	c Pas de temps adaptatif
197	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
198	if (n.GT.1) then
199	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
200	s dtvr,'n=',n
201	endif
202	do indice=1,n
203	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
204	end do
205	else if(iadv(iq).eq.13) then
206	stop 'advtrac : schema non parallelise'
207	c Pas de temps adaptatif
208	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
209	if (n.GT.1) then
210	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
211	s dtvr,'n=',n
212	endif
213	do indice=1,n
214	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
215	end do
216	c ----------------------------------------------------------------
217	c Schema de pente SLOPES
218	c ----------------------------------------------------------------
219	else if (iadv(iq).eq.20) then
220	stop 'advtrac : schema non parallelise'
221
222	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
223
224	c ----------------------------------------------------------------
225	c Schema de Prather
226	c ----------------------------------------------------------------
227	else if (iadv(iq).eq.30) then
228	stop 'advtrac : schema non parallelise'
229	c Pas de temps adaptatif
230	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
231	if (n.GT.1) then
232	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
233	s dtvr,'n=',n
234	endif
235	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
236	s n,dtbon)
237	c ----------------------------------------------------------------
238	c Schemas PPM Lin et Rood
239	c ----------------------------------------------------------------
240	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
241	s iadv(iq).LE.18)) then
242
243	stop 'advtrac : schema non parallelise'
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
340	end DO
341
342	1234 CONTINUE
343	c$OMP BARRIER
344
345	if (planet_type=="earth") then
346
347	ijb=ij_begin
348	ije=ij_end
349
350	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
351	DO l = 1, llm
352	DO ij = ijb, ije
353	finmasse(ij,l) = p(ij,l) - p(ij,l+1)
354	ENDDO
355	ENDDO
356	c$OMP END DO
357
358	CALL qminimum_loc( q, 2, finmasse )
359
360	endif ! of if (planet_type=="earth")
361
362	RETURN
363	END
364

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

Download in other formats:

Original Format