Context Navigation

advtrac_p.F90 @ 3000

Last change on this file since 3000 was 2296, checked in by mvals, 5 years ago

Mars GCM:
dyn3dpar: Implementation in the parallel version of the dynamics of the dynamical transport of the isotopes using the same scheme as the one
implemented by Camille Risi in the Earth GCM (see LMDZ6/libf/dyn3dmem, and Risi et al. 2009: genealogy of the tracers+transport of the isotopic
ratio). This is added to prepare the future implementation of the HDO cycle in the GCM. These changes had been already added in the sequential part.
dyn3d: corrections to prevent from dividing by zero in the case of the use of the isotopes scheme.
traceur.def.isotopes: example of how to write the traceur.def if you want to use the new dynamics with the genealogy of the isotopes: air is the
father of all, H2O is the father of HDO.
MV

File size: 18.4 KB

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

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

Download in other formats:

Original Format