Context Navigation

advtrac.F90 @ 2039

Last change on this file since 2039 was 1987, checked in by Ehouarn Millour, 10 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: 14.3 KB

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

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

Download in other formats:

Original Format