Context Navigation

advtrac.F @ 5198

Last change on this file since 5198 was 1446, checked in by Ehouarn Millour, 14 years ago

Implemented modifications to enable running with only one tracer for planet types different from "earth". Rem: If flag 'planet_type' is set to "earth" (default behaviour) then there must be at least 2 tracers for the dynamics to function properly.

These updates do not induce any changes in model outputs with respect to previous revisions.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 13.8 KB

Rev	Line
[524]	1	!
[1279]	2	! $Id: advtrac.F 1446 2010-10-22 09:27:25Z idelkadi $
[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
[1403]	18	USE control_mod
	19
[1146]	20
[524]	21	IMPLICIT NONE
	22	c
	23	#include "dimensions.h"
	24	#include "paramet.h"
	25	#include "comconst.h"
	26	#include "comvert.h"
	27	#include "comdissip.h"
	28	#include "comgeom2.h"
	29	#include "logic.h"
	30	#include "temps.h"
	31	#include "ener.h"
	32	#include "description.h"
[1279]	33	#include "iniprint.h"
[524]	34
	35	c-------------------------------------------------------------------
	36	c Arguments
	37	c-------------------------------------------------------------------
	38	c Ajout PPM
	39	c--------------------------------------------------------
	40	REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
	41	c--------------------------------------------------------
	42	INTEGER iapptrac
	43	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
[1146]	44	REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
[524]	45	REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
	46	REAL pk(ip1jmp1,llm)
[960]	47	REAL flxw(ip1jmp1,llm)
[524]	48
	49	c-------------------------------------------------------------
	50	c Variables locales
	51	c-------------------------------------------------------------
	52
	53	REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
	54	REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
	55	REAL 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	EXTERNAL minmax
	61	SAVE iadvtr, massem, pbaruc, pbarvc
	62	DATA iadvtr/0/
	63	c----------------------------------------------------------
	64	c Rajouts pour PPM
	65	c----------------------------------------------------------
	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)
[1146]	71	REAL qppm(iim*jjp1,llm,nqtot)
[524]	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
[1279]	78	integer,save :: countcfl=0
	79	real cflx(ip1jmp1,llm)
	80	real cfly(ip1jm,llm)
	81	real cflz(ip1jmp1,llm)
	82	real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
[524]	83
	84	IF(iadvtr.EQ.0) THEN
	85	CALL initial0(ijp1llm,pbaruc)
	86	CALL initial0(ijmllm,pbarvc)
	87	ENDIF
	88
	89	c accumulation des flux de masse horizontaux
	90	DO l=1,llm
	91	DO ij = 1,ip1jmp1
	92	pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
	93	ENDDO
	94	DO ij = 1,ip1jm
	95	pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
	96	ENDDO
	97	ENDDO
	98
	99	c selection de la masse instantannee des mailles avant le transport.
	100	IF(iadvtr.EQ.0) THEN
	101
	102	CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
	103	ccc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
	104	c
	105	ENDIF
	106
	107	iadvtr = iadvtr+1
	108	iapptrac = iadvtr
	109
	110
	111	c Test pour savoir si on advecte a ce pas de temps
	112	IF ( iadvtr.EQ.iapp_tracvl ) THEN
	113
	114	cc .. Modif P.Le Van ( 20/12/97 ) ....
	115	cc
	116
	117	c traitement des flux de masse avant advection.
	118	c 1. calcul de w
	119	c 2. groupement des mailles pres du pole.
	120
	121	CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
	122
	123	! ... Flux de masse diaganostiques traceurs
[1403]	124	flxw = wg / REAL(iapp_tracvl)
[524]	125
	126	c test sur l'eventuelle creation de valeurs negatives de la masse
	127	DO l=1,llm-1
	128	DO ij = iip2+1,ip1jm
	129	zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l)
	130	s - pbarvg(ij-iip1,l) + pbarvg(ij,l)
	131	s + wg(ij,l+1) - wg(ij,l)
	132	ENDDO
	133	CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
	134	DO ij = iip2,ip1jm
	135	zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
	136	ENDDO
	137
	138
	139	CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
	140
	141	IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
	142	PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin,
	143	s ' MAX:', zdpmax
	144	ENDIF
	145
	146	ENDDO
	147
[1279]	148
[524]	149	c-------------------------------------------------------------------
[1279]	150	! Calcul des criteres CFL en X, Y et Z
	151	c-------------------------------------------------------------------
	152
	153	if (countcfl == 0. ) then
	154	cflxmax(:)=0.
	155	cflymax(:)=0.
	156	cflzmax(:)=0.
	157	endif
	158
	159	countcfl=countcfl+iapp_tracvl
	160	cflx(:,:)=0.
	161	cfly(:,:)=0.
	162	cflz(:,:)=0.
	163	do l=1,llm
	164	do ij=iip2,ip1jm-1
	165	if (pbarug(ij,l)>=0.) then
	166	cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
	167	else
	168	cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
	169	endif
	170	enddo
	171	enddo
	172	do l=1,llm
	173	do ij=iip2,ip1jm-1,iip1
	174	cflx(ij+iip1,l)=cflx(ij,l)
	175	enddo
	176	enddo
	177
	178	do l=1,llm
	179	do ij=1,ip1jm
	180	if (pbarvg(ij,l)>=0.) then
	181	cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
	182	else
	183	cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
	184	endif
	185	enddo
	186	enddo
	187
	188	do l=2,llm
	189	do ij=1,ip1jm
	190	if (wg(ij,l)>=0.) then
	191	cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
	192	else
	193	cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
	194	endif
	195	enddo
	196	enddo
	197
	198	do l=1,llm
	199	cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
	200	cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
	201	cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
	202	enddo
	203
	204	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	205	! Par defaut, on sort le diagnostic des CFL tous les jours.
	206	! Si on veut le sortir a chaque pas d'advection en cas de plantage
	207	! if (countcfl==iapp_tracvl) then
	208	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	209	if (countcfl==day_step) then
	210	do l=1,llm
	211	write(lunout,*) 'L, CFLmax '
	212	s ,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l))
	213	enddo
	214	countcfl=0
	215	endif
	216
	217	c-------------------------------------------------------------------
[524]	218	c Advection proprement dite (Modification Le Croller (07/2001)
	219	c-------------------------------------------------------------------
	220
	221	c----------------------------------------------------
	222	c Calcul des moyennes basées sur la masse
	223	c----------------------------------------------------
	224	call massbar(massem,massebx,masseby)
	225
	226	c-----------------------------------------------------------
	227	c Appel des sous programmes d'advection
	228	c-----------------------------------------------------------
[1146]	229	do iq=1,nqtot
[524]	230	c call clock(t_initial)
	231	if(iadv(iq) == 0) cycle
	232	c ----------------------------------------------------------------
	233	c Schema de Van Leer I MUSCL
	234	c ----------------------------------------------------------------
	235	if(iadv(iq).eq.10) THEN
	236	call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
	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