Context Navigation

advtrac.F @ 1340

Last change on this file since 1340 was 1279, checked in by Laurent Fairhead, 15 years ago
Merged LMDZ4-dev branch changes r1241:1278 into the trunk Running trunk and LMDZ4-dev in LMDZOR configuration on local machine (sequential) and SX8 (4-proc) yields identical results (restart and restartphy are identical binarily) Log history from r1241 to r1278 is available by switching to source:LMDZ4/branches/LMDZ4-dev-20091210
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 13.7 KB

Rev	Line
[524]	1	!
[1279]	2	! $Id: advtrac.F 1279 2009-12-10 09:02:56Z fairhead $
[524]	3	!
	4	c
	5	c
	6	SUBROUTINE advtrac(pbaru,pbarv ,
	7	* p, masse,q,iapptrac,teta,
	8	* flxw,
[762]	9	* pk)
[524]	10	c Auteur : F. Hourdin
	11	c
	12	c Modif. P. Le Van (20/12/97)
	13	c F. Codron (10/99)
	14	c D. Le Croller (07/2001)
	15	c M.A Filiberti (04/2002)
	16	c
[1146]	17	USE infotrac
	18
[524]	19	IMPLICIT NONE
	20	c
	21	#include "dimensions.h"
	22	#include "paramet.h"
	23	#include "comconst.h"
	24	#include "comvert.h"
	25	#include "comdissip.h"
	26	#include "comgeom2.h"
	27	#include "logic.h"
	28	#include "temps.h"
	29	#include "control.h"
	30	#include "ener.h"
	31	#include "description.h"
[1279]	32	#include "iniprint.h"
[524]	33
	34	c-------------------------------------------------------------------
	35	c Arguments
	36	c-------------------------------------------------------------------
	37	c Ajout PPM
	38	c--------------------------------------------------------
	39	REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
	40	c--------------------------------------------------------
	41	INTEGER iapptrac
	42	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
[1146]	43	REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
[524]	44	REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
	45	REAL pk(ip1jmp1,llm)
[960]	46	REAL flxw(ip1jmp1,llm)
[524]	47
	48	c-------------------------------------------------------------
	49	c Variables locales
	50	c-------------------------------------------------------------
	51
	52	REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
	53	REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
	54	REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
	55	REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
	56	INTEGER iadvtr
	57	INTEGER ij,l,iq,iiq
	58	REAL zdpmin, zdpmax
	59	EXTERNAL minmax
	60	SAVE iadvtr, massem, pbaruc, pbarvc
	61	DATA iadvtr/0/
	62	c----------------------------------------------------------
	63	c Rajouts pour PPM
	64	c----------------------------------------------------------
	65	INTEGER indice,n
	66	REAL dtbon ! Pas de temps adaptatif pour que CFL<1
	67	REAL CFLmaxz,aaa,bbb ! CFL maximum
	68	REAL psppm(iim,jjp1) ! pression au sol
	69	REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
[1146]	70	REAL qppm(iim*jjp1,llm,nqtot)
[524]	71	REAL fluxwppm(iim,jjp1,llm)
	72	REAL apppm(llmp1), bpppm(llmp1)
	73	LOGICAL dum,fill
	74	DATA fill/.true./
	75	DATA dum/.true./
	76
[1279]	77	integer,save :: countcfl=0
	78	real cflx(ip1jmp1,llm)
	79	real cfly(ip1jm,llm)
	80	real cflz(ip1jmp1,llm)
	81	real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
[524]	82
	83	IF(iadvtr.EQ.0) THEN
	84	CALL initial0(ijp1llm,pbaruc)
	85	CALL initial0(ijmllm,pbarvc)
	86	ENDIF
	87
	88	c accumulation des flux de masse horizontaux
	89	DO l=1,llm
	90	DO ij = 1,ip1jmp1
	91	pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
	92	ENDDO
	93	DO ij = 1,ip1jm
	94	pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
	95	ENDDO
	96	ENDDO
	97
	98	c selection de la masse instantannee des mailles avant le transport.
	99	IF(iadvtr.EQ.0) THEN
	100
	101	CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
	102	ccc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
	103	c
	104	ENDIF
	105
	106	iadvtr = iadvtr+1
	107	iapptrac = iadvtr
	108
	109
	110	c Test pour savoir si on advecte a ce pas de temps
	111	IF ( iadvtr.EQ.iapp_tracvl ) THEN
	112
	113	cc .. Modif P.Le Van ( 20/12/97 ) ....
	114	cc
	115
	116	c traitement des flux de masse avant advection.
	117	c 1. calcul de w
	118	c 2. groupement des mailles pres du pole.
	119
	120	CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
	121
	122	! ... Flux de masse diaganostiques traceurs
	123	flxw = wg / FLOAT(iapp_tracvl)
	124
	125	c test sur l'eventuelle creation de valeurs negatives de la masse
	126	DO l=1,llm-1
	127	DO ij = iip2+1,ip1jm
	128	zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l)
	129	s - pbarvg(ij-iip1,l) + pbarvg(ij,l)
	130	s + wg(ij,l+1) - wg(ij,l)
	131	ENDDO
	132	CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
	133	DO ij = iip2,ip1jm
	134	zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
	135	ENDDO
	136
	137
	138	CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
	139
	140	IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
	141	PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin,
	142	s ' MAX:', zdpmax
	143	ENDIF
	144
	145	ENDDO
	146
[1279]	147
[524]	148	c-------------------------------------------------------------------
[1279]	149	! Calcul des criteres CFL en X, Y et Z
	150	c-------------------------------------------------------------------
	151
	152	if (countcfl == 0. ) then
	153	cflxmax(:)=0.
	154	cflymax(:)=0.
	155	cflzmax(:)=0.
	156	endif
	157
	158	countcfl=countcfl+iapp_tracvl
	159	cflx(:,:)=0.
	160	cfly(:,:)=0.
	161	cflz(:,:)=0.
	162	do l=1,llm
	163	do ij=iip2,ip1jm-1
	164	if (pbarug(ij,l)>=0.) then
	165	cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
	166	else
	167	cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
	168	endif
	169	enddo
	170	enddo
	171	do l=1,llm
	172	do ij=iip2,ip1jm-1,iip1
	173	cflx(ij+iip1,l)=cflx(ij,l)
	174	enddo
	175	enddo
	176
	177	do l=1,llm
	178	do ij=1,ip1jm
	179	if (pbarvg(ij,l)>=0.) then
	180	cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
	181	else
	182	cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
	183	endif
	184	enddo
	185	enddo
	186
	187	do l=2,llm
	188	do ij=1,ip1jm
	189	if (wg(ij,l)>=0.) then
	190	cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
	191	else
	192	cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
	193	endif
	194	enddo
	195	enddo
	196
	197	do l=1,llm
	198	cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
	199	cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
	200	cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
	201	enddo
	202
	203	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	204	! Par defaut, on sort le diagnostic des CFL tous les jours.
	205	! Si on veut le sortir a chaque pas d'advection en cas de plantage
	206	! if (countcfl==iapp_tracvl) then
	207	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	208	if (countcfl==day_step) then
	209	do l=1,llm
	210	write(lunout,*) 'L, CFLmax '
	211	s ,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l))
	212	enddo
	213	countcfl=0
	214	endif
	215
	216	c-------------------------------------------------------------------
[524]	217	c Advection proprement dite (Modification Le Croller (07/2001)
	218	c-------------------------------------------------------------------
	219
	220	c----------------------------------------------------
	221	c Calcul des moyennes basées sur la masse
	222	c----------------------------------------------------
	223	call massbar(massem,massebx,masseby)
	224
	225	c-----------------------------------------------------------
	226	c Appel des sous programmes d'advection
	227	c-----------------------------------------------------------
[1146]	228	do iq=1,nqtot
[524]	229	c call clock(t_initial)
	230	if(iadv(iq) == 0) cycle
	231	c ----------------------------------------------------------------
	232	c Schema de Van Leer I MUSCL
	233	c ----------------------------------------------------------------
	234	if(iadv(iq).eq.10) THEN
	235	call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
	236
	237
	238	c ----------------------------------------------------------------
	239	c Schema "pseudo amont" + test sur humidite specifique
	240	C pour la vapeur d'eau. F. Codron
	241	c ----------------------------------------------------------------
	242	else if(iadv(iq).eq.14) then
	243	c
	244	CALL vlspltqs( q(1,1,1), 2., massem, wg ,
	245	* pbarug,pbarvg,dtvr,p,pk,teta )
	246	c ----------------------------------------------------------------
	247	c Schema de Frederic Hourdin
	248	c ----------------------------------------------------------------
	249	else if(iadv(iq).eq.12) then
	250	c Pas de temps adaptatif
	251	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	252	if (n.GT.1) then
	253	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
	254	s dtvr,'n=',n
	255	endif
	256	do indice=1,n
	257	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
	258	end do
	259	else if(iadv(iq).eq.13) then
	260	c Pas de temps adaptatif
	261	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	262	if (n.GT.1) then
	263	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
	264	s dtvr,'n=',n
	265	endif
	266	do indice=1,n
	267	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
	268	end do
	269	c ----------------------------------------------------------------
	270	c Schema de pente SLOPES
	271	c ----------------------------------------------------------------
	272	else if (iadv(iq).eq.20) then
	273	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
	274
	275	c ----------------------------------------------------------------
	276	c Schema de Prather
	277	c ----------------------------------------------------------------
	278	else if (iadv(iq).eq.30) then
	279	c Pas de temps adaptatif
	280	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	281	if (n.GT.1) then
	282	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
	283	s dtvr,'n=',n
	284	endif
	285	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
	286	s n,dtbon)
[960]	287
[524]	288	c ----------------------------------------------------------------
	289	c Schemas PPM Lin et Rood
	290	c ----------------------------------------------------------------
	291	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
	292	s iadv(iq).LE.18)) then
	293
	294	c Test sur le flux horizontal
	295	c Pas de temps adaptatif
	296	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	297	if (n.GT.1) then
	298	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
	299	s dtvr,'n=',n
	300	endif
	301	c Test sur le flux vertical
	302	CFLmaxz=0.
	303	do l=2,llm
	304	do ij=iip2,ip1jm
	305	aaa=wg(ij,l)*dtvr/massem(ij,l)
	306	CFLmaxz=max(CFLmaxz,aaa)
	307	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
	308	CFLmaxz=max(CFLmaxz,bbb)
	309	enddo
	310	enddo
	311	if (CFLmaxz.GE.1) then
	312	write(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
	313	endif
	314
	315	c-----------------------------------------------------------
	316	c Ss-prg interface LMDZ.4->PPM3d
	317	c-----------------------------------------------------------
	318
	319	call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
	320	s apppm,bpppm,massebx,masseby,pbarug,pbarvg,
	321	s unatppm,vnatppm,psppm)
	322
	323	do indice=1,n
	324	c---------------------------------------------------------------------
	325	c VL (version PPM) horiz. et PPM vert.
	326	c---------------------------------------------------------------------
	327	if (iadv(iq).eq.11) then
	328	c Ss-prg PPM3d de Lin
	329	call ppm3d(1,qppm(1,1,iq),
	330	s psppm,psppm,
	331	s unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
	332	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
	333	s fill,dum,220.)
	334
	335	c----------------------------------------------------------------------
	336	c Monotonic PPM
	337	c----------------------------------------------------------------------
	338	else if (iadv(iq).eq.16) then
	339	c Ss-prg PPM3d de Lin
	340	call ppm3d(1,qppm(1,1,iq),
	341	s psppm,psppm,
	342	s unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
	343	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
	344	s fill,dum,220.)
	345	c---------------------------------------------------------------------
	346
	347	c---------------------------------------------------------------------
	348	c Semi Monotonic PPM
	349	c---------------------------------------------------------------------
	350	else if (iadv(iq).eq.17) then
	351	c Ss-prg PPM3d de Lin
	352	call ppm3d(1,qppm(1,1,iq),
	353	s psppm,psppm,
	354	s unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
	355	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
	356	s fill,dum,220.)
	357	c---------------------------------------------------------------------
	358
	359	c---------------------------------------------------------------------
	360	c Positive Definite PPM
	361	c---------------------------------------------------------------------
	362	else if (iadv(iq).eq.18) then
	363	c Ss-prg PPM3d de Lin
	364	call ppm3d(1,qppm(1,1,iq),
	365	s psppm,psppm,
	366	s unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
	367	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
	368	s fill,dum,220.)
	369	c---------------------------------------------------------------------
	370	endif
	371	enddo
	372	c-----------------------------------------------------------------
	373	c Ss-prg interface PPM3d-LMDZ.4
	374	c-----------------------------------------------------------------
	375	call interpost(q(1,1,iq),qppm(1,1,iq))
	376	endif
	377	c----------------------------------------------------------------------
	378
	379	c-----------------------------------------------------------------
	380	c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
	381	c et Nord j=1
	382	c-----------------------------------------------------------------
	383
	384	c call traceurpole(q(1,1,iq),massem)
	385
	386	c calcul du temps cpu pour un schema donne
	387
	388	c call clock(t_final)
[566]	389	cym tps_cpu=t_final-t_initial
	390	cym cpuadv(iq)=cpuadv(iq)+tps_cpu
[524]	391
	392	end DO
	393
	394
	395	c------------------------------------------------------------------
	396	c on reinitialise a zero les flux de masse cumules
	397	c---------------------------------------------------
	398	iadvtr=0
	399
	400	ENDIF ! if iadvtr.EQ.iapp_tracvl
	401
	402	RETURN
	403	END
	404

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

Download in other formats:

Original Format