Context Navigation

prather.F @ 5456

Last change on this file since 5456 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

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, 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