Context Navigation

advtrac_p.F @ 1119

Last change on this file since 1119 was 985, checked in by Laurent Fairhead, 16 years ago
Mise a jour de dyn3dpar par rapport a dyn3d, inclusion OpenMP et filtre FFT YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.6 KB

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

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

Download in other formats:

Original Format