Context Navigation

prather.F @ 1481

Last change on this file since 1481 was 764, checked in by Laurent Fairhead, 17 years ago
Merge entre la version V3_conv et le HEAD YM, JG, LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 10.6 KB

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

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

Download in other formats:

Original Format