Context Navigation

advtrac.F90 @ 3553

Last change on this file since 3553 was 1508, checked in by emillour, 9 years ago

Common dynamics:
Updates in the dynamics (seq and ) to keep up with updates
in LMDZ5 (up to LMDZ5 trunk, rev 2325):
IMPORTANT: Modifications for isotopes are only done in dyn3d, not in dyn3dpar

as in LMDZ5 these modifications were done in dyn3dmem.
Related LMDZ5 revisions are r2270 and r2281

in dynlonlat_phylonlat:
add module "grid_atob_m.F90" (a regridding utility so far only used by phylmd/ce0l.F90, used to be dyn3d_common/grid_atob.F)

in misc:
follow up updates on wxios.F (add missing_val module variable)

in dyn3d_common:
pression.F => pression.F90
misc_mod.F90: moved from misc to dyn3d_common
added new iso_verif_dyn.F
covcont.F => covcont.F90
infotrac.F90 : add handling of isotopes (reading of corresponding traceur.def for planets not implemented)
dynetat0.F => dynetat0.F90 with some code factorization
dynredem.F => dynredem.F90 with some code factorization
added dynredem_mod.F90: routines used by dynredem
iniacademic.F90 : added isotopes-related initialization for Earth case

in dyn3d:
added check_isotopes.F
modified (isotopes) advtrac.F90, caladvtrac.F
guide_mod.F90: ported updates
leapfrog.F : (isotopes) updates (NB: call integrd with nqtot tracers)
qminimium.F : adaptations for isotopes (copied over, except that #include comvert.h is not needed).
vlsplt.F: adaptations for isotopes (copied over, except than #include logic.h, comvert.h not needed, and replace "include comconst.h" with use comconst_mod, ONLY: pi)
vlspltqs.F : same as vlsplt.F, but also keeping added modification for CP(T)

in dyn3dpar:
leapfrog_p.F: remove unecessary #ifdef CPP_EARTH cpp flag. and call integrd_p with nqtot tracers (only important for Earth)
dynredem_p.F => dynredem_p.F90 and some code factorization
and no isotopes-relates changes in dyn3dpar (since these changes have been made in LMDZ5 dyn3dmem).

File size: 14.7 KB

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

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

Download in other formats:

Original Format