Context Navigation

advtrac_p.F90 @ 2404

Last change on this file since 2404 was 1987, checked in by Ehouarn Millour, 11 years ago
Add updating pressure, mass and Exner function (ie: all variables which depend on surface pressure) after adding physics tendencies (which include a surface pressure tendency). Note that this change induces slight changes in GCM results with respect to previous svn version of the code, even if surface pressure tendency is zero (because of recomputation of polar values as an average over polar points on the dynamics grid). EM
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 Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 17.0 KB

Line
1	! $Id: advtrac_p.F90 1987 2014-02-24 15:05:47Z fhourdin $
2
3	SUBROUTINE advtrac_p(pbaru,pbarv , p, masse,q,iapptrac,teta, flxw, pk)
4
5	! Auteur : F. Hourdin
6	!
7	! Modif. P. Le Van (20/12/97)
8	! F. Codron (10/99)
9	! D. Le Croller (07/2001)
10	! M.A Filiberti (04/2002)
11	!
12	USE parallel_lmdz
13	USE Write_Field_p
14	USE Bands
15	USE mod_hallo
16	USE Vampir
17	USE times
18	USE infotrac, ONLY: nqtot, iadv
19	USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
20	IMPLICIT NONE
21	!
22	include "dimensions.h"
23	include "paramet.h"
24	include "comconst.h"
25	include "comvert.h"
26	include "comdissip.h"
27	include "comgeom2.h"
28	include "logic.h"
29	include "temps.h"
30	include "ener.h"
31	include "description.h"
32
33	!-------------------------------------------------------------------
34	! Arguments
35	!-------------------------------------------------------------------
36	INTEGER,INTENT(OUT) :: iapptrac
37	REAL,INTENT(IN) :: pbaru(ip1jmp1,llm)
38	REAL,INTENT(IN) :: pbarv(ip1jm,llm)
39	REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
40	REAL,INTENT(IN) :: masse(ip1jmp1,llm)
41	REAL,INTENT(IN) :: p( ip1jmp1,llmp1 )
42	REAL,INTENT(IN) :: teta(ip1jmp1,llm)
43	REAL,INTENT(IN) :: pk(ip1jmp1,llm)
44	REAL,INTENT(OUT) :: flxw(ip1jmp1,llm)
45	!-------------------------------------------------------------------
46	! Ajout PPM
47	!--------------------------------------------------------
48	REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
49	!-------------------------------------------------------------
50	! Variables locales
51	!-------------------------------------------------------------
52
53	REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
54	REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
55	REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
56	REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
57	INTEGER iadvtr
58	INTEGER ij,l,iq,iiq
59	REAL zdpmin, zdpmax
60	SAVE iadvtr, massem, pbaruc, pbarvc
61	DATA iadvtr/0/
62	!$OMP THREADPRIVATE(iadvtr)
63	!----------------------------------------------------------
64	! Rajouts pour PPM
65	!----------------------------------------------------------
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,nqtot)
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	REAL,SAVE :: finmasse(ip1jmp1,llm)
78	integer ijb,ije,ijb_u,ijb_v,ije_u,ije_v,j
79	type(Request) :: Request_vanleer
80	REAL,SAVE :: p_tmp( ip1jmp1,llmp1 )
81	REAL,SAVE :: teta_tmp(ip1jmp1,llm)
82	REAL,SAVE :: pk_tmp(ip1jmp1,llm)
83
84	ijb_u=ij_begin
85	ije_u=ij_end
86
87	ijb_v=ij_begin-iip1
88	ije_v=ij_end
89	if (pole_nord) ijb_v=ij_begin
90	if (pole_sud) ije_v=ij_end-iip1
91
92	IF(iadvtr.EQ.0) THEN
93	! CALL initial0(ijp1llm,pbaruc)
94	! CALL initial0(ijmllm,pbarvc)
95	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
96	DO l=1,llm
97	pbaruc(ijb_u:ije_u,l)=0.
98	pbarvc(ijb_v:ije_v,l)=0.
99	ENDDO
100	!$OMP END DO NOWAIT
101	ENDIF
102
103	! accumulation des flux de masse horizontaux
104	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
105	DO l=1,llm
106	DO ij = ijb_u,ije_u
107	pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
108	ENDDO
109	DO ij = ijb_v,ije_v
110	pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
111	ENDDO
112	ENDDO
113	!$OMP END DO NOWAIT
114
115	! selection de la masse instantannee des mailles avant le transport.
116	IF(iadvtr.EQ.0) THEN
117
118	! CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
119	ijb=ij_begin
120	ije=ij_end
121
122	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
123	DO l=1,llm
124	massem(ijb:ije,l)=masse(ijb:ije,l)
125	ENDDO
126	!$OMP END DO NOWAIT
127
128	!cc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
129	!
130	ENDIF ! of IF(iadvtr.EQ.0)
131
132	iadvtr = iadvtr+1
133
134	!$OMP MASTER
135	iapptrac = iadvtr
136	!$OMP END MASTER
137
138	! Test pour savoir si on advecte a ce pas de temps
139
140	IF ( iadvtr.EQ.iapp_tracvl ) THEN
141	!$OMP MASTER
142	call suspend_timer(timer_caldyn)
143	!$OMP END MASTER
144
145	ijb=ij_begin
146	ije=ij_end
147
148
149	!c .. Modif P.Le Van ( 20/12/97 ) ....
150	!c
151
152	! traitement des flux de masse avant advection.
153	! 1. calcul de w
154	! 2. groupement des mailles pres du pole.
155
156	CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
157
158	!$OMP BARRIER
159
160	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
161	DO l=1,llmp1
162	p_tmp(ijb:ije,l)=p(ijb:ije,l)
163	ENDDO
164	!$OMP END DO NOWAIT
165
166	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
167	DO l=1,llm
168	pk_tmp(ijb:ije,l)=pk(ijb:ije,l)
169	teta_tmp(ijb:ije,l)=teta(ijb:ije,l)
170	ENDDO
171	!$OMP END DO NOWAIT
172
173	!$OMP MASTER
174	call VTb(VTHallo)
175	!$OMP END MASTER
176
177	call Register_SwapFieldHallo(pbarug,pbarug,ip1jmp1,llm, &
178	jj_Nb_vanleer,0,0,Request_vanleer)
179	call Register_SwapFieldHallo(pbarvg,pbarvg,ip1jm,llm, &
180	jj_Nb_vanleer,1,0,Request_vanleer)
181	call Register_SwapFieldHallo(massem,massem,ip1jmp1,llm, &
182	jj_Nb_vanleer,0,0,Request_vanleer)
183	call Register_SwapFieldHallo(wg,wg,ip1jmp1,llm, &
184	jj_Nb_vanleer,0,0,Request_vanleer)
185	call Register_SwapFieldHallo(teta_tmp,teta_tmp,ip1jmp1,llm, &
186	jj_Nb_vanleer,1,1,Request_vanleer)
187	call Register_SwapFieldHallo(p_tmp,p_tmp,ip1jmp1,llmp1, &
188	jj_Nb_vanleer,1,1,Request_vanleer)
189	call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm, &
190	jj_Nb_vanleer,1,1,Request_vanleer)
191	do j=1,nqtot
192	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, &
193	jj_nb_vanleer,0,0,Request_vanleer)
194	enddo
195
196	call SendRequest(Request_vanleer)
197	!$OMP BARRIER
198	call WaitRequest(Request_vanleer)
199
200
201	!$OMP BARRIER
202	!$OMP MASTER
203	call SetDistrib(jj_nb_vanleer)
204	call VTe(VTHallo)
205	call VTb(VTadvection)
206	call start_timer(timer_vanleer)
207	!$OMP END MASTER
208	!$OMP BARRIER
209
210	! ... Flux de masse diaganostiques traceurs
211	ijb=ij_begin
212	ije=ij_end
213	flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/REAL(iapp_tracvl)
214
215	! test sur l'eventuelle creation de valeurs negatives de la masse
216	ijb=ij_begin
217	ije=ij_end
218	if (pole_nord) ijb=ij_begin+iip1
219	if (pole_sud) ije=ij_end-iip1
220
221	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
222	DO l=1,llm-1
223	DO ij = ijb+1,ije
224	zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) &
225	- pbarvg(ij-iip1,l) + pbarvg(ij,l) &
226	+ wg(ij,l+1) - wg(ij,l)
227	ENDDO
228
229	! CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
230	! ym ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
231
232	do ij=ijb,ije-iip1+1,iip1
233	zdp(ij)=zdp(ij+iip1-1)
234	enddo
235
236	DO ij = ijb,ije
237	zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
238	ENDDO
239
240
241	! CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
242	! ym ---> eventuellement a revoir
243	CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
244
245	IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
246	PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin, &
247	' MAX:', zdpmax
248	ENDIF
249
250	ENDDO
251	!$OMP END DO NOWAIT
252
253	!-------------------------------------------------------------------
254	! Advection proprement dite (Modification Le Croller (07/2001)
255	!-------------------------------------------------------------------
256
257	!----------------------------------------------------
258	! Calcul des moyennes basées sur la masse
259	!----------------------------------------------------
260
261	!ym ----> Normalement, inutile pour les schémas classiques
262	!ym ----> Revérifier lors de la parallélisation des autres schemas
263
264	!ym call massbar_p(massem,massebx,masseby)
265
266	call vlspltgen_p( q,iadv, 2., massem, wg , &
267	pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
268
269
270	GOTO 1234
271	!-----------------------------------------------------------
272	! Appel des sous programmes d'advection
273	!-----------------------------------------------------------
274	do iq=1,nqtot
275	! call clock(t_initial)
276	if(iadv(iq) == 0) cycle
277	! ----------------------------------------------------------------
278	! Schema de Van Leer I MUSCL
279	! ----------------------------------------------------------------
280	if(iadv(iq).eq.10) THEN
281
282	call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
283
284	! ----------------------------------------------------------------
285	! Schema "pseudo amont" + test sur humidite specifique
286	! pour la vapeur d'eau. F. Codron
287	! ----------------------------------------------------------------
288	else if(iadv(iq).eq.14) then
289	!
290	!ym stop 'advtrac : appel à vlspltqs :schema non parallelise'
291	CALL vlspltqs_p( q(1,1,1), 2., massem, wg , &
292	pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
293	! ----------------------------------------------------------------
294	! Schema de Frederic Hourdin
295	! ----------------------------------------------------------------
296	else if(iadv(iq).eq.12) then
297	stop 'advtrac : schema non parallelise'
298	! Pas de temps adaptatif
299	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
300	if (n.GT.1) then
301	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
302	dtvr,'n=',n
303	endif
304	do indice=1,n
305	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
306	end do
307	else if(iadv(iq).eq.13) then
308	stop 'advtrac : schema non parallelise'
309	! Pas de temps adaptatif
310	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
311	if (n.GT.1) then
312	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
313	dtvr,'n=',n
314	endif
315	do indice=1,n
316	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
317	end do
318	! ----------------------------------------------------------------
319	! Schema de pente SLOPES
320	! ----------------------------------------------------------------
321	else if (iadv(iq).eq.20) then
322	stop 'advtrac : schema non parallelise'
323
324	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
325
326	! ----------------------------------------------------------------
327	! Schema de Prather
328	! ----------------------------------------------------------------
329	else if (iadv(iq).eq.30) then
330	stop 'advtrac : schema non parallelise'
331	! Pas de temps adaptatif
332	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
333	if (n.GT.1) then
334	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
335	dtvr,'n=',n
336	endif
337	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
338	n,dtbon)
339	! ----------------------------------------------------------------
340	! Schemas PPM Lin et Rood
341	! ----------------------------------------------------------------
342	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
343	iadv(iq).LE.18)) then
344
345	stop 'advtrac : schema non parallelise'
346
347	! Test sur le flux horizontal
348	! Pas de temps adaptatif
349	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
350	if (n.GT.1) then
351	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
352	dtvr,'n=',n
353	endif
354	! Test sur le flux vertical
355	CFLmaxz=0.
356	do l=2,llm
357	do ij=iip2,ip1jm
358	aaa=wg(ij,l)*dtvr/massem(ij,l)
359	CFLmaxz=max(CFLmaxz,aaa)
360	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
361	CFLmaxz=max(CFLmaxz,bbb)
362	enddo
363	enddo
364	if (CFLmaxz.GE.1) then
365	write(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
366	endif
367
368	!-----------------------------------------------------------
369	! Ss-prg interface LMDZ.4->PPM3d
370	!-----------------------------------------------------------
371
372	call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
373	apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
374	unatppm,vnatppm,psppm)
375
376	do indice=1,n
377	!----------------------------------------------------------------
378	! VL (version PPM) horiz. et PPM vert.
379	!----------------------------------------------------------------
380	if (iadv(iq).eq.11) then
381	! Ss-prg PPM3d de Lin
382	call ppm3d(1,qppm(1,1,iq), &
383	psppm,psppm, &
384	unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
385	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
386	fill,dum,220.)
387
388	!-------------------------------------------------------------
389	! Monotonic PPM
390	!-------------------------------------------------------------
391	else if (iadv(iq).eq.16) then
392	! Ss-prg PPM3d de Lin
393	call ppm3d(1,qppm(1,1,iq), &
394	psppm,psppm, &
395	unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
396	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
397	fill,dum,220.)
398	!-------------------------------------------------------------
399
400	!-------------------------------------------------------------
401	! Semi Monotonic PPM
402	!-------------------------------------------------------------
403	else if (iadv(iq).eq.17) then
404	! Ss-prg PPM3d de Lin
405	call ppm3d(1,qppm(1,1,iq), &
406	psppm,psppm, &
407	unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
408	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
409	fill,dum,220.)
410	!-------------------------------------------------------------
411
412	!-------------------------------------------------------------
413	! Positive Definite PPM
414	!-------------------------------------------------------------
415	else if (iadv(iq).eq.18) then
416	! Ss-prg PPM3d de Lin
417	call ppm3d(1,qppm(1,1,iq), &
418	psppm,psppm, &
419	unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
420	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
421	fill,dum,220.)
422	!-------------------------------------------------------------
423	endif
424	enddo
425	!-----------------------------------------------------------------
426	! Ss-prg interface PPM3d-LMDZ.4
427	!-----------------------------------------------------------------
428	call interpost(q(1,1,iq),qppm(1,1,iq))
429	endif
430	!----------------------------------------------------------------------
431
432	!-----------------------------------------------------------------
433	! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
434	! et Nord j=1
435	!-----------------------------------------------------------------
436
437	! call traceurpole(q(1,1,iq),massem)
438
439	! calcul du temps cpu pour un schema donne
440
441	! call clock(t_final)
442	!ym tps_cpu=t_final-t_initial
443	!ym cpuadv(iq)=cpuadv(iq)+tps_cpu
444
445	end DO
446
447	1234 CONTINUE
448	!$OMP BARRIER
449
450	if (planet_type=="earth") then
451
452	ijb=ij_begin
453	ije=ij_end
454
455	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
456	DO l = 1, llm
457	DO ij = ijb, ije
458	finmasse(ij,l) = p(ij,l) - p(ij,l+1)
459	ENDDO
460	ENDDO
461	!$OMP END DO
462
463	CALL qminimum_p( q, 2, finmasse )
464	endif ! of if (planet_type=="earth")
465
466	!------------------------------------------------------------------
467	! on reinitialise a zero les flux de masse cumules
468	!---------------------------------------------------
469	! iadvtr=0
470
471	!$OMP MASTER
472	call VTe(VTadvection)
473	call stop_timer(timer_vanleer)
474	call VTb(VThallo)
475	!$OMP END MASTER
476
477	do j=1,nqtot
478	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, &
479	jj_nb_caldyn,0,0,Request_vanleer)
480	enddo
481
482	call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm, &
483	jj_nb_caldyn,0,0,Request_vanleer)
484
485	call SendRequest(Request_vanleer)
486	!$OMP BARRIER
487	call WaitRequest(Request_vanleer)
488
489	!$OMP BARRIER
490	!$OMP MASTER
491	call SetDistrib(jj_nb_caldyn)
492	call VTe(VThallo)
493	call resume_timer(timer_caldyn)
494	!$OMP END MASTER
495	!$OMP BARRIER
496	iadvtr=0
497	ENDIF ! if iadvtr.EQ.iapp_tracvl
498
499	END SUBROUTINE advtrac_p

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

Download in other formats:

Original Format