Context Navigation

advtrac.F @ 3400

Last change on this file since 3400 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

tous les FLOAT() sont remplacés par des REAL()
tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

all FLOAT() instructions have been replaced by REAL() instructions
all STOP instructions in phylmd have been replaced by calls to abort_gcm
the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
$Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 13.8 KB

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

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

Download in other formats:

Original Format