Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

prather.f90

Last change on this file was 5159, checked in by abarral, 3 months ago
Put dimensions.h and paramet.h into modules
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: 9.7 KB

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

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

Download in other formats:

Original Format