Context Navigation

advtrac.F90 @ 4052

Last change on this file since 4052 was 4050, checked in by dcugnet, 3 years ago

Second commit for new tracers.

include most of the keys in the tracers descriptor vector "tracers(:)".
fix in phylmdiso/cv3_routines: fq_* variables were used where their fxt_* counterparts were expected.
multiple IF(nqdesc(iq)>0) and IF(nqfils(iq)>0) tests suppressed, because they are not needed: "do ... enddo" loops with 0 upper bound are not executed.
remove French accents from comments (encoding problem) in phylmdiso/cv3_routines and phylmdiso/cv30_routines.
modifications in "isotopes_verif_mod", where the call to function "iso_verif_tag17_q_deltad_chn" in "iso_verif_tag17_q_deltad_chn" was not detected at linking stage, although defined in the same module (?).

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

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

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

Download in other formats:

Original Format