Context Navigation

advtrac.F90 @ 2600

Last change on this file since 2600 was 2600, checked in by Ehouarn Millour, 9 years ago
Cleanup in the dynamics: turn comvert.h into module comvert_mod.F90 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.8 KB

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

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

Download in other formats:

Original Format