Context Navigation

prather.F @ 653

Last change on this file since 653 was 630, checked in by Laurent Fairhead, 20 years ago
Import d'une version parallele de la dynamique YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 10.6 KB

Rev	Line
[630]	1	!
	2	! $Header$
	3	!
	4	SUBROUTINE prather (q,w,masse,pbaru,pbarv,nt,dt)
	5	IMPLICIT NONE
	6
	7	c=======================================================================
	8	c Adaptation LMDZ: A.Armengaud (LGGE)
	9	c ----------------
	10	c
	11	c ************************************************
	12	c Transport des traceurs par la methode de prather
	13	c Ref :
	14	c
	15	c ************************************************
	16	c q,w,pext,pbaru et pbarv : arguments d'entree pour le s-pg
	17	c
	18	c=======================================================================
	19
	20
	21	#include "dimensions.h"
	22	#include "paramet.h"
	23	#include "comconst.h"
	24	#include "comvert.h"
	25	#include "comgeom2.h"
	26
	27	c Arguments:
	28	c ----------
	29	INTEGER iq,nt
	30	REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
	31	REAL masse(iip1,jjp1,llm)
	32	REAL q( iip1,jjp1,llm,0:9)
	33	REAL w( ip1jmp1,llm )
	34	integer ordre,ilim
	35
	36	c Local:
	37	c ------
	38	LOGICAL limit
	39	real zq(iip1,jjp1,llm)
	40	REAL sm ( iip1,jjp1, llm )
	41	REAL s0( iip1,jjp1,llm ), sx( iip1,jjp1,llm )
	42	REAL sy( iip1,jjp1,llm ), sz( iip1,jjp1,llm )
	43	REAL sxx( iip1,jjp1,llm)
	44	REAL sxy( iip1,jjp1,llm)
	45	REAL sxz( iip1,jjp1,llm)
	46	REAL syy( iip1,jjp1,llm )
	47	REAL syz( iip1,jjp1,llm )
	48	REAL szz( iip1,jjp1,llm ),zz
	49	INTEGER i,j,l,indice
	50	real sxn(iip1),sxs(iip1)
	51
	52	real sinlon(iip1),sinlondlon(iip1)
	53	real coslon(iip1),coslondlon(iip1)
	54	real qmin,qmax
	55	save qmin,qmax
	56	save sinlon,coslon,sinlondlon,coslondlon
	57	real dyn1,dyn2,dys1,dys2,qpn,qps,dqzpn,dqzps
	58	real masn,mass
	59	c
	60	REAL SSUM
	61	integer ismax,ismin
	62	EXTERNAL SSUM, convflu,ismin,ismax
	63	logical first
	64	save first
	65	EXTERNAL advxp,advyp,advzp
	66
	67
	68	data first/.true./
	69	data qmin,qmax/-1.e33,1.e33/
	70
	71
	72	c==========================================================================
	73	c==========================================================================
	74	c MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
	75	c==========================================================================
	76	c==========================================================================
	77	REAL dt
	78	c==========================================================================
	79	limit = .TRUE.
	80
	81	if(first) then
	82	print*,'SCHEMA PRATHER'
	83	first=.false.
	84	do i=2,iip1
	85	coslon(i)=cos(rlonv(i))
	86	sinlon(i)=sin(rlonv(i))
	87	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
	88	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
	89	enddo
	90	coslon(1)=coslon(iip1)
	91	coslondlon(1)=coslondlon(iip1)
	92	sinlon(1)=sinlon(iip1)
	93	sinlondlon(1)=sinlondlon(iip1)
	94
	95	DO l = 1,llm
	96	DO j = 1,jjp1
	97	DO i = 1,iip1
	98	q( i,j,l,1 )=0.
	99	q( i,j,l,2)=0.
	100	q( i,j,l,3)=0.
	101	q( i,j,l,4)=0.
	102	q( i,j,l,5)=0.
	103	q( i,j,l,6)=0.
	104	q( i,j,l,7)=0.
	105	q( i,j,l,8)=0.
	106	q( i,j,l,9)=0.
	107	ENDDO
	108	ENDDO
	109	ENDDO
	110	endif
	111	c Fin modif Fred
	112
	113	c *** On calcule la masse d'air en kg
	114
	115	DO l = 1,llm
	116	DO j = 1,jjp1
	117	DO i = 1,iip1
	118	sm( i,j,llm+1-l ) =masse(i,j,l)
	119	ENDDO
	120	ENDDO
	121	ENDDO
	122
	123	c *** q contient les qqtes de traceur avant l'advection
	124
	125	c *** Affectation des tableaux S a partir de Q
	126
	127	DO l = 1,llm
	128	DO j = 1,jjp1
	129	DO i = 1,iip1
	130	s0( i,j,l) = q ( i,j,llm+1-l,0 )*sm(i,j,l)
	131	sx( i,j,l) = q( i,j,llm+1-l,1 )*sm(i,j,l)
	132	sy( i,j,l) = q( i,j,llm+1-l,2)*sm(i,j,l)
	133	sz( i,j,l) = q( i,j,llm+1-l,3)*sm(i,j,l)
	134	sxx( i,j,l) = q( i,j,llm+1-l,4)*sm(i,j,l)
	135	sxy( i,j,l) = q( i,j,llm+1-l,5)*sm(i,j,l)
	136	sxz( i,j,l) = q( i,j,llm+1-l,6)*sm(i,j,l)
	137	syy( i,j,l) = q( i,j,llm+1-l,7)*sm(i,j,l)
	138	syz( i,j,l) = q( i,j,llm+1-l,8)*sm(i,j,l)
	139	szz( i,j,l) = q( i,j,llm+1-l,9)*sm(i,j,l)
	140	ENDDO
	141	ENDDO
	142	ENDDO
	143	c *** Appel des subroutines d'advection en X, en Y et en Z
	144	c *** Advection avec "time-splitting"
	145
	146	c-----------------------------------------------------------
	147	do indice =1,nt
	148	call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz
	149	. ,sxx,sxy,sxz,syy,syz,szz,1 )
	150	end do
	151	do l=1,llm
	152	do i=1,iip1
	153	sy(i,1,l)=0.
	154	sy(i,jjp1,l)=0.
	155	enddo
	156	enddo
	157	c---------------------------------------------------------
	158	call advyp( limit,.5dtnt,pbarv,sm,s0,sx,sy,sz
	159	. ,sxx,sxy,sxz,syy,syz,szz,1 )
	160	c---------------------------------------------------------
	161
	162	c---------------------------------------------------------
	163	do j=1,jjp1
	164	do i=1,iip1
	165	sz(i,j,1)=0.
	166	sz(i,j,llm)=0.
	167	sxz(i,j,1)=0.
	168	sxz(i,j,llm)=0.
	169	syz(i,j,1)=0.
	170	syz(i,j,llm)=0.
	171	szz(i,j,1)=0.
	172	szz(i,j,llm)=0.
	173	enddo
	174	enddo
	175	call advzp( limit,dt*nt,w,sm,s0,sx,sy,sz
	176	. ,sxx,sxy,sxz,syy,syz,szz,1 )
	177	do l=1,llm
	178	do i=1,iip1
	179	sy(i,1,l)=0.
	180	sy(i,jjp1,l)=0.
	181	enddo
	182	enddo
	183
	184	c---------------------------------------------------------
	185
	186	c---------------------------------------------------------
	187	call advyp( limit,.5dtnt,pbarv,sm,s0,sx,sy,sz
	188	. ,sxx,sxy,sxz,syy,syz,szz,1 )
	189	c---------------------------------------------------------
	190	DO l = 1,llm
	191	DO j = 1,jjp1
	192	s0( iip1,j,l)=s0( 1,j,l )
	193	sx( iip1,j,l)=sx( 1,j,l )
	194	sy( iip1,j,l)=sy( 1,j,l )
	195	sz( iip1,j,l)=sz( 1,j,l )
	196	sxx( iip1,j,l)=sxx( 1,j,l )
	197	sxy( iip1,j,l)=sxy( 1,j,l)
	198	sxz( iip1,j,l)=sxz( 1,j,l )
	199	syy( iip1,j,l)=syy( 1,j,l )
	200	syz( iip1,j,l)=syz( 1,j,l)
	201	szz( iip1,j,l)=szz( 1,j,l )
	202	ENDDO
	203	ENDDO
	204	do indice=1,nt
	205	call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz
	206	. ,sxx,sxy,sxz,syy,syz,szz,1 )
	207	end do
	208	c---------------------------------------------------------
	209	c---------------------------------------------------------
	210	c *** On repasse les S dans la variable qpr
	211	c *** On repasse les S dans la variable q directement 14/10/94
	212
	213	DO l = 1,llm
	214	DO j = 1,jjp1
	215	DO i = 1,iip1
	216	q( i,j,llm+1-l,0 )=s0( i,j,l )/sm(i,j,l)
	217	q( i,j,llm+1-l,1 ) = sx( i,j,l )/sm(i,j,l)
	218	q( i,j,llm+1-l,2 ) = sy( i,j,l )/sm(i,j,l)
	219	q( i,j,llm+1-l,3 ) = sz( i,j,l )/sm(i,j,l)
	220	q( i,j,llm+1-l,4 ) = sxx( i,j,l )/sm(i,j,l)
	221	q( i,j,llm+1-l,5 ) = sxy( i,j,l )/sm(i,j,l)
	222	q( i,j,llm+1-l,6 ) = sxz( i,j,l )/sm(i,j,l)
	223	q( i,j,llm+1-l,7 ) = syy( i,j,l )/sm(i,j,l)
	224	q( i,j,llm+1-l,8 ) = syz( i,j,l )/sm(i,j,l)
	225	q( i,j,llm+1-l,9 ) = szz( i,j,l )/sm(i,j,l)
	226	ENDDO
	227	ENDDO
	228	ENDDO
	229
	230	c---------------------------------------------------------
	231	c go to 777
	232	c filtrages aux poles
	233
	234	c Traitements specifiques au pole
	235
	236	c filtrages aux poles
	237	DO l=1,llm
	238	c filtrages aux poles
	239	masn=ssum(iim,sm(1,1,l),1)
	240	mass=ssum(iim,sm(1,jjp1,l),1)
	241	qpn=ssum(iim,s0(1,1,l),1)/masn
	242	qps=ssum(iim,s0(1,jjp1,l),1)/mass
	243	dqzpn=ssum(iim,sz(1,1,l),1)/masn
	244	dqzps=ssum(iim,sz(1,jjp1,l),1)/mass
	245	do i=1,iip1
	246	q( i,1,llm+1-l,3)=dqzpn
	247	q( i,jjp1,llm+1-l,3)=dqzps
	248	q( i,1,llm+1-l,0)=qpn
	249	q( i,jjp1,llm+1-l,0)=qps
	250	enddo
	251	c enddo
	252	c print*,'qpn',qpn,'qps',qps
	253	c print*,'dqzpn',dqzpn,'dqzps',dqzps
	254	c enddo
	255	dyn1=0.
	256	dys1=0.
	257	dyn2=0.
	258	dys2=0.
	259	do i=1,iim
	260	zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)
	261	dyn1=dyn1+sinlondlon(i)*zz
	262	dyn2=dyn2+coslondlon(i)*zz
	263	zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)
	264	dys1=dys1+sinlondlon(i)*zz
	265	dys2=dys2+coslondlon(i)*zz
	266	enddo
	267	do i=1,iim
	268	q(i,1,llm+1-l,2)=
	269	$ (sinlon(i)dyn1+coslon(i)dyn2)/2.
	270	q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)
	271	$ +q(i,1,llm+1-l,2)
	272	q(i,jjp1,llm+1-l,2)=
	273	$ (sinlon(i)dys1+coslon(i)dys2)/2.
	274	q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
	275	$ -q(i,jjp1,llm+1-l,2)
	276	enddo
	277	q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)
	278	q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)
	279	do i=1,iim
	280	sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)
	281	sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)
	282	enddo
	283	sxn(iip1)=sxn(1)
	284	sxs(iip1)=sxs(1)
	285	do i=1,iim
	286	q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))
	287	q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))
	288	END DO
	289	q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)
	290	q(1,jjp1,llm+1-l,1)=
	291	$ q(iip1,jjp1,llm+1-l,1)
	292	enddo
	293	do l=1,llm
	294	do i=1,iim
	295	q( i,1,llm+1-l,4)=0.
	296	q( i,jjp1,llm+1-l,4)=0.
	297	q( i,1,llm+1-l,5)=0.
	298	q( i,jjp1,llm+1-l,5)=0.
	299	q( i,1,llm+1-l,6)=0.
	300	q( i,jjp1,llm+1-l,6)=0.
	301	q( i,1,llm+1-l,7)=0.
	302	q( i,jjp1,llm+1-l,7)=0.
	303	q( i,1,llm+1-l,8)=0.
	304	q( i,jjp1,llm+1-l,8)=0.
	305	q( i,1,llm+1-l,9)=0.
	306	q( i,jjp1,llm+1-l,9)=0.
	307	enddo
	308	ENDDO
	309
	310	777 continue
	311	c
	312	c bouclage en longitude
	313	do l=1,llm
	314	do j=1,jjp1
	315	q(iip1,j,l,0)=q(1,j,l,0)
	316	q(iip1,j,llm+1-l,0)=q(1,j,llm+1-l,0)
	317	q(iip1,j,llm+1-l,1)=q(1,j,llm+1-l,1)
	318	q(iip1,j,llm+1-l,2)=q(1,j,llm+1-l,2)
	319	q(iip1,j,llm+1-l,3)=q(1,j,llm+1-l,3)
	320	q(iip1,j,llm+1-l,4)=q(1,j,llm+1-l,4)
	321	q(iip1,j,llm+1-l,5)=q(1,j,llm+1-l,5)
	322	q(iip1,j,llm+1-l,6)=q(1,j,llm+1-l,6)
	323	q(iip1,j,llm+1-l,7)=q(1,j,llm+1-l,7)
	324	q(iip1,j,llm+1-l,8)=q(1,j,llm+1-l,8)
	325	q(iip1,j,llm+1-l,9)=q(1,j,llm+1-l,9)
	326	enddo
	327	enddo
	328	DO l = 1,llm
	329	DO j = 2,jjm
	330	DO i = 1,iip1
	331	IF (q(i,j,l,0).lt.0.) THEN
	332	PRINT*,'------------ BIP-----------'
	333	PRINT*,'S0(',i,j,l,')=',q(i,j,l,0),
	334	$ q(i,j-1,l,0)
	335	PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)
	336	PRINT*,'SY(',i,j,l,')=',q(i,j,l,2),
	337	$ q(i,j-1,l,2)
	338	PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)
	339	c PRINT*,' PBL EN SORTIE D'' ADVZP'
	340	q(i,j,l,0)=0.
	341	c STOP
	342	ENDIF
	343	ENDDO
	344	ENDDO
	345	do j=1,jjp1,jjm
	346	do i=1,iip1
	347	IF (q(i,j,l,0).lt.0.) THEN
	348	PRINT*,'------------ BIP 2-----------'
	349	PRINT*,'S0(',i,j,l,')=',q(i,j,l,0)
	350	PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)
	351	PRINT*,'SY(',i,j,l,')=',q(i,j,l,2)
	352	PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)
	353
	354	q(i,j,l,0)=0.
	355	c STOP
	356	ENDIF
	357	enddo
	358	enddo
	359	ENDDO
	360	RETURN
	361	END

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

Download in other formats:

Original Format