Context Navigation

advtrac.F90 @ 2256

Last change on this file since 2256 was 2239, checked in by Ehouarn Millour, 10 years ago

Reorganizing physics/dynamics interface:

what is related to dynamics-physics interface is now in a seperate directory: dynlmdz_phy* for physics in phy*
1d model and related dependencies (including a couple from "dynamics", set up as symbolic links) is now in subdirectory "dyn1d" of phy*.
"bibio" directory is now "misc" and should only contain autonomous utilities.
"cosp" is now a subdirectory of phylmd.

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.2 KB

Rev	Line
[1279]	1	! $Id: advtrac.F90 2239 2015-03-23 07:27:30Z fairhead $
[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
[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	!-----------------------------------------------------------
	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