Context Navigation

advtrac.F90 @ 2519

Last change on this file since 2519 was 2286, checked in by Ehouarn Millour, 9 years ago

fix in dyn3d the array out of bound issue (in qminimum) that was only corrected in dyn3dmem in rev 2285.
comment out many invasive debug writes.

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 2286 2015-05-20 13:27:07Z 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
[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
[2239]	81	pbaruc(:,:)=0
	82	pbarvc(:,:)=0
[1549]	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	!-----------------------------------------------------------
[2270]	225
	226	if (ok_iso_verif) then
[2281]	227	write(,) 'advtrac 227'
	228	call check_isotopes_seq(q,ip1jmp1,'advtrac 162')
[2270]	229	endif !if (ok_iso_verif) then
	230
	231	do iq=1,nqperes
[1549]	232	! call clock(t_initial)
[524]	233	if(iadv(iq) == 0) cycle
[1549]	234	! ----------------------------------------------------------------
	235	! Schema de Van Leer I MUSCL
	236	! ----------------------------------------------------------------
[524]	237	if(iadv(iq).eq.10) THEN
[2270]	238	! CRisi: on fait passer tout q pour avoir acces aux fils
	239
[2286]	240	!write(,) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)
[2270]	241	call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq)
	242
[1549]	243	! ----------------------------------------------------------------
	244	! Schema "pseudo amont" + test sur humidite specifique
	245	! pour la vapeur d'eau. F. Codron
	246	! ----------------------------------------------------------------
[524]	247	else if(iadv(iq).eq.14) then
[1549]	248	!
[2286]	249	!write(,) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
[2270]	250	CALL vlspltqs( q, 2., massem, wg , &
	251	pbarug,pbarvg,dtvr,p,pk,teta,iq)
	252
[1549]	253	! ----------------------------------------------------------------
	254	! Schema de Frederic Hourdin
	255	! ----------------------------------------------------------------
[524]	256	else if(iadv(iq).eq.12) then
[1549]	257	! Pas de temps adaptatif
[524]	258	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	259	if (n.GT.1) then
[1549]	260	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
	261	dtvr,'n=',n
[524]	262	endif
	263	do indice=1,n
[1549]	264	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
[524]	265	end do
	266	else if(iadv(iq).eq.13) then
[1549]	267	! Pas de temps adaptatif
[524]	268	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	269	if (n.GT.1) then
[1549]	270	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
	271	dtvr,'n=',n
[524]	272	endif
[1549]	273	do indice=1,n
	274	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
	275	end do
	276	! ----------------------------------------------------------------
	277	! Schema de pente SLOPES
	278	! ----------------------------------------------------------------
[524]	279	else if (iadv(iq).eq.20) then
[1549]	280	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
[524]	281
[1549]	282	! ----------------------------------------------------------------
	283	! Schema de Prather
	284	! ----------------------------------------------------------------
[524]	285	else if (iadv(iq).eq.30) then
[1549]	286	! Pas de temps adaptatif
[524]	287	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	288	if (n.GT.1) then
[1549]	289	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
	290	dtvr,'n=',n
[524]	291	endif
[1549]	292	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
	293	n,dtbon)
[960]	294
[1549]	295	! ----------------------------------------------------------------
	296	! Schemas PPM Lin et Rood
	297	! ----------------------------------------------------------------
	298	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
	299	iadv(iq).LE.18)) then
[524]	300
[1549]	301	! Test sur le flux horizontal
	302	! Pas de temps adaptatif
	303	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	304	if (n.GT.1) then
	305	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
	306	dtvr,'n=',n
	307	endif
	308	! Test sur le flux vertical
	309	CFLmaxz=0.
	310	do l=2,llm
	311	do ij=iip2,ip1jm
	312	aaa=wg(ij,l)*dtvr/massem(ij,l)
	313	CFLmaxz=max(CFLmaxz,aaa)
	314	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
	315	CFLmaxz=max(CFLmaxz,bbb)
	316	enddo
[524]	317	enddo
[1549]	318	if (CFLmaxz.GE.1) then
	319	write(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
	320	endif
[524]	321
[1549]	322	!-----------------------------------------------------------
	323	! Ss-prg interface LMDZ.4->PPM3d
	324	!-----------------------------------------------------------
[524]	325
[1549]	326	call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
	327	apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
	328	unatppm,vnatppm,psppm)
[524]	329
[1549]	330	do indice=1,n
	331	!----------------------------------------------------------------
	332	! VL (version PPM) horiz. et PPM vert.
	333	!----------------------------------------------------------------
	334	if (iadv(iq).eq.11) then
	335	! Ss-prg PPM3d de Lin
	336	call ppm3d(1,qppm(1,1,iq), &
	337	psppm,psppm, &
	338	unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
	339	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
	340	fill,dum,220.)
[524]	341
[1549]	342	!-------------------------------------------------------------
	343	! Monotonic PPM
	344	!-------------------------------------------------------------
	345	else if (iadv(iq).eq.16) then
	346	! Ss-prg PPM3d de Lin
	347	call ppm3d(1,qppm(1,1,iq), &
	348	psppm,psppm, &
	349	unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
	350	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
	351	fill,dum,220.)
	352	!-------------------------------------------------------------
[524]	353
[1549]	354	!-------------------------------------------------------------
	355	! Semi Monotonic PPM
	356	!-------------------------------------------------------------
	357	else if (iadv(iq).eq.17) then
	358	! Ss-prg PPM3d de Lin
	359	call ppm3d(1,qppm(1,1,iq), &
	360	psppm,psppm, &
	361	unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
	362	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
	363	fill,dum,220.)
	364	!-------------------------------------------------------------
[524]	365
[1549]	366	!-------------------------------------------------------------
	367	! Positive Definite PPM
	368	!-------------------------------------------------------------
	369	else if (iadv(iq).eq.18) then
	370	! Ss-prg PPM3d de Lin
	371	call ppm3d(1,qppm(1,1,iq), &
	372	psppm,psppm, &
	373	unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
	374	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
	375	fill,dum,220.)
	376	!-------------------------------------------------------------
	377	endif
	378	enddo
	379	!-----------------------------------------------------------------
	380	! Ss-prg interface PPM3d-LMDZ.4
	381	!-----------------------------------------------------------------
	382	call interpost(q(1,1,iq),qppm(1,1,iq))
	383	endif
	384	!----------------------------------------------------------------------
[524]	385
[1549]	386	!-----------------------------------------------------------------
	387	! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
	388	! et Nord j=1
	389	!-----------------------------------------------------------------
[524]	390
[1549]	391	! call traceurpole(q(1,1,iq),massem)
[524]	392
[1549]	393	! calcul du temps cpu pour un schema donne
[524]	394
[1549]	395	! call clock(t_final)
	396	!ym tps_cpu=t_final-t_initial
	397	!ym cpuadv(iq)=cpuadv(iq)+tps_cpu
[524]	398
[1549]	399	end DO
[524]	400
[2270]	401	if (ok_iso_verif) then
[2281]	402	write(,) 'advtrac 402'
[2270]	403	call check_isotopes_seq(q,ip1jmp1,'advtrac 397')
	404	endif !if (ok_iso_verif) then
[524]	405
[1549]	406	!------------------------------------------------------------------
	407	! on reinitialise a zero les flux de masse cumules
	408	!---------------------------------------------------
	409	iadvtr=0
[524]	410
[1549]	411	ENDIF ! if iadvtr.EQ.iapp_tracvl
[524]	412
[1549]	413	END SUBROUTINE advtrac

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

Download in other formats:

Original Format