Context Navigation

advtrac.F90 @ 1019

Last change on this file since 1019 was 1019, checked in by emillour, 11 years ago

Common dynamics; keep up with updates (seq and ) in LMDZ5 (up tio rev 1845):

General stuff:
makelmdz_fcm: add options -j # (compile using # threads) and -full, and to keep up

with Earth model, possibility to compile with various versions of orchidee

bld.cfg: adaptations to enable compiling using multiple threads
build_gcm: adaptations to enable compiling using multiple threads
makelmdz: keep up with Earth model: possibility to compile with various versions of orchidee + cosmetic changes + library directory name change

bibio:
wxios.F90 : Added for possible future use of XIOS library

filtrez:
mkl_dft_type.f90 & mkl_dfti.f90 : MKL (for MKL FFT) interface definitions
filtreg_mod : limit use of FFT to parallel mode
mod_filtre_fft.F90 & mod_filtre_fft_lov.F90 : swich to use parallel_lmdz

dyn3d:
abort_gcm.F : add things for xios
advtrac.F90 : minor change in CFL outputs
ce0l.F90 : indicesol.h is now module indice_sol_mod
comvert.h : cosmetic change on comments
gcm.F : add xios and use module indice_sol_mod (for INCA)
inigeom.F : move two computations outside loop

dyn3dpar:
parallel.F90 => parallel_lmdz.F90 : and change all the "use parallel" into "use parallel_lmdz" in all files in dyn3dpar
comvert.h : cosmetic change on comments
gcm.F : add xios and use module indice_sol_mod (for INCA)
leapfrog_p.F : add xios + correction for times in Newtonian case
ce0l.F90 : indicesol.h is now module indice_sol_mod
inigeom.F : move two computations outside loop

File size: 14.1 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	!
	11	USE infotrac
	12	USE control_mod
[1]	13
	14
[270]	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"
[1]	28
[270]	29	!-------------------------------------------------------------------
	30	! Arguments
	31	!-------------------------------------------------------------------
	32	! Ajout PPM
	33	!--------------------------------------------------------
	34	REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
	35	!--------------------------------------------------------
	36	INTEGER iapptrac
	37	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
	38	REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
	39	REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
	40	REAL pk(ip1jmp1,llm)
	41	REAL flxw(ip1jmp1,llm)
[1]	42
[270]	43	!-------------------------------------------------------------
	44	! Variables locales
	45	!-------------------------------------------------------------
[1]	46
[270]	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./
[1]	71
[270]	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)
[1]	77
[270]	78	IF(iadvtr.EQ.0) THEN
	79	CALL initial0(ijp1llm,pbaruc)
	80	CALL initial0(ijmllm,pbarvc)
	81	ENDIF
[1]	82
[270]	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
[1]	92
[270]	93	! selection de la masse instantannee des mailles avant le transport.
	94	IF(iadvtr.EQ.0) THEN
[1]	95
[270]	96	CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
	97	!cc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
	98	!
	99	ENDIF
[1]	100
[270]	101	iadvtr = iadvtr+1
	102	iapptrac = iadvtr
[1]	103
	104
[270]	105	! Test pour savoir si on advecte a ce pas de temps
	106	IF ( iadvtr.EQ.iapp_tracvl ) THEN
[1]	107
[270]	108	!c .. Modif P.Le Van ( 20/12/97 ) ....
	109	!c
[1]	110
[270]	111	! traitement des flux de masse avant advection.
	112	! 1. calcul de w
	113	! 2. groupement des mailles pres du pole.
[1]	114
[270]	115	CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
[1]	116
[270]	117	! ... Flux de masse diaganostiques traceurs
	118	flxw = wg / REAL(iapp_tracvl)
[1]	119
[270]	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
[1]	131
	132
[270]	133	CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
[1]	134
[270]	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
[1]	139
[270]	140	ENDDO
[1]	141
	142
[270]	143	!-------------------------------------------------------------------
	144	! Calcul des criteres CFL en X, Y et Z
	145	!-------------------------------------------------------------------
[1]	146
[270]	147	if (countcfl == 0. ) then
	148	cflxmax(:)=0.
	149	cflymax(:)=0.
	150	cflzmax(:)=0.
	151	endif
[1]	152
[270]	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
[1]	171
[270]	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
[1]	181
[270]	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
[1]	198	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[270]	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
[1]	202	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[270]	203	if (countcfl==day_step) then
	204	do l=1,llm
[1019]	205	write(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), &
	206	cflzmax(l)
[270]	207	enddo
	208	countcfl=0
	209	endif
[1]	210
[270]	211	!-------------------------------------------------------------------
	212	! Advection proprement dite (Modification Le Croller (07/2001)
	213	!-------------------------------------------------------------------
[1]	214
[270]	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	!-----------------------------------------------------------
	223	do iq=1,nqtot
	224	! call clock(t_initial)
[1]	225	if(iadv(iq) == 0) cycle
[270]	226	! ----------------------------------------------------------------
	227	! Schema de Van Leer I MUSCL
	228	! ----------------------------------------------------------------
[1]	229	if(iadv(iq).eq.10) THEN
[270]	230	call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
[1]	231
[270]	232	! ----------------------------------------------------------------
	233	! Schema "pseudo amont" + test sur humidite specifique
	234	! pour la vapeur d'eau. F. Codron
	235	! ----------------------------------------------------------------
[1]	236	else if(iadv(iq).eq.14) then
[270]	237	!
	238	CALL vlspltqs( q(1,1,1), 2., massem, wg , &
	239	pbarug,pbarvg,dtvr,p,pk,teta )
	240	! ----------------------------------------------------------------
	241	! Schema de Frederic Hourdin
	242	! ----------------------------------------------------------------
[1]	243	else if(iadv(iq).eq.12) then
[270]	244	! Pas de temps adaptatif
[1]	245	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	246	if (n.GT.1) then
[270]	247	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
	248	dtvr,'n=',n
[1]	249	endif
	250	do indice=1,n
[270]	251	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
[1]	252	end do
	253	else if(iadv(iq).eq.13) then
[270]	254	! Pas de temps adaptatif
[1]	255	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	256	if (n.GT.1) then
[270]	257	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
	258	dtvr,'n=',n
[1]	259	endif
[270]	260	do indice=1,n
	261	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
	262	end do
	263	! ----------------------------------------------------------------
	264	! Schema de pente SLOPES
	265	! ----------------------------------------------------------------
[1]	266	else if (iadv(iq).eq.20) then
[270]	267	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
[1]	268
[270]	269	! ----------------------------------------------------------------
	270	! Schema de Prather
	271	! ----------------------------------------------------------------
[1]	272	else if (iadv(iq).eq.30) then
[270]	273	! Pas de temps adaptatif
[1]	274	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	275	if (n.GT.1) then
[270]	276	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
	277	dtvr,'n=',n
[1]	278	endif
[270]	279	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
	280	n,dtbon)
[1]	281
[270]	282	! ----------------------------------------------------------------
	283	! Schemas PPM Lin et Rood
	284	! ----------------------------------------------------------------
	285	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
	286	iadv(iq).LE.18)) then
[1]	287
[270]	288	! Test sur le flux horizontal
	289	! Pas de temps adaptatif
	290	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	291	if (n.GT.1) then
	292	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
	293	dtvr,'n=',n
	294	endif
	295	! Test sur le flux vertical
	296	CFLmaxz=0.
	297	do l=2,llm
	298	do ij=iip2,ip1jm
	299	aaa=wg(ij,l)*dtvr/massem(ij,l)
	300	CFLmaxz=max(CFLmaxz,aaa)
	301	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
	302	CFLmaxz=max(CFLmaxz,bbb)
	303	enddo
[1]	304	enddo
[270]	305	if (CFLmaxz.GE.1) then
	306	write(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
	307	endif
[1]	308
[270]	309	!-----------------------------------------------------------
	310	! Ss-prg interface LMDZ.4->PPM3d
	311	!-----------------------------------------------------------
[1]	312
[270]	313	call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
	314	apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
	315	unatppm,vnatppm,psppm)
[1]	316
[270]	317	do indice=1,n
	318	!----------------------------------------------------------------
	319	! VL (version PPM) horiz. et PPM vert.
	320	!----------------------------------------------------------------
	321	if (iadv(iq).eq.11) then
	322	! Ss-prg PPM3d de Lin
	323	call ppm3d(1,qppm(1,1,iq), &
	324	psppm,psppm, &
	325	unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
	326	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
	327	fill,dum,220.)
[1]	328
[270]	329	!-------------------------------------------------------------
	330	! Monotonic PPM
	331	!-------------------------------------------------------------
	332	else if (iadv(iq).eq.16) then
	333	! Ss-prg PPM3d de Lin
	334	call ppm3d(1,qppm(1,1,iq), &
	335	psppm,psppm, &
	336	unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
	337	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
	338	fill,dum,220.)
	339	!-------------------------------------------------------------
[1]	340
[270]	341	!-------------------------------------------------------------
	342	! Semi Monotonic PPM
	343	!-------------------------------------------------------------
	344	else if (iadv(iq).eq.17) then
	345	! Ss-prg PPM3d de Lin
	346	call ppm3d(1,qppm(1,1,iq), &
	347	psppm,psppm, &
	348	unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
	349	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
	350	fill,dum,220.)
	351	!-------------------------------------------------------------
[1]	352
[270]	353	!-------------------------------------------------------------
	354	! Positive Definite PPM
	355	!-------------------------------------------------------------
	356	else if (iadv(iq).eq.18) then
	357	! Ss-prg PPM3d de Lin
	358	call ppm3d(1,qppm(1,1,iq), &
	359	psppm,psppm, &
	360	unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
	361	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
	362	fill,dum,220.)
	363	!-------------------------------------------------------------
	364	endif
	365	enddo
	366	!-----------------------------------------------------------------
	367	! Ss-prg interface PPM3d-LMDZ.4
	368	!-----------------------------------------------------------------
	369	call interpost(q(1,1,iq),qppm(1,1,iq))
	370	endif
	371	!----------------------------------------------------------------------
[1]	372
[270]	373	!-----------------------------------------------------------------
	374	! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
	375	! et Nord j=1
	376	!-----------------------------------------------------------------
[1]	377
[270]	378	! call traceurpole(q(1,1,iq),massem)
[1]	379
[270]	380	! calcul du temps cpu pour un schema donne
[1]	381
[270]	382	! call clock(t_final)
	383	!ym tps_cpu=t_final-t_initial
	384	!ym cpuadv(iq)=cpuadv(iq)+tps_cpu
[1]	385
[270]	386	end DO
[1]	387
	388
[270]	389	!------------------------------------------------------------------
	390	! on reinitialise a zero les flux de masse cumules
	391	!---------------------------------------------------
	392	iadvtr=0
[1]	393
[270]	394	ENDIF ! if iadvtr.EQ.iapp_tracvl
[1]	395
[270]	396	END SUBROUTINE advtrac

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

Download in other formats:

Original Format