Context Navigation

prather.F @ 3000

Last change on this file since 3000 was 1422, checked in by milmd, 10 years ago
In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).
File size: 10.6 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE prather (q,w,masse,pbaru,pbarv,nt,dt)
5
6	USE comconst_mod, ONLY: pi
7
8	IMPLICIT NONE
9
10	c=======================================================================
11	c Adaptation LMDZ: A.Armengaud (LGGE)
12	c ----------------
13	c
14	c ************************************************
15	c Transport des traceurs par la methode de prather
16	c Ref :
17	c
18	c ************************************************
19	c q,w,pext,pbaru et pbarv : arguments d'entree pour le s-pg
20	c
21	c=======================================================================
22
23
24	#include "dimensions.h"
25	#include "paramet.h"
26	#include "comgeom2.h"
27
28	c Arguments:
29	c ----------
30	INTEGER iq,nt
31	REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
32	REAL masse(iip1,jjp1,llm)
33	REAL q( iip1,jjp1,llm,0:9)
34	REAL w( ip1jmp1,llm )
35	integer ordre,ilim
36
37	c Local:
38	c ------
39	LOGICAL limit
40	real zq(iip1,jjp1,llm)
41	REAL sm ( iip1,jjp1, llm )
42	REAL s0( iip1,jjp1,llm ), sx( iip1,jjp1,llm )
43	REAL sy( iip1,jjp1,llm ), sz( iip1,jjp1,llm )
44	REAL sxx( iip1,jjp1,llm)
45	REAL sxy( iip1,jjp1,llm)
46	REAL sxz( iip1,jjp1,llm)
47	REAL syy( iip1,jjp1,llm )
48	REAL syz( iip1,jjp1,llm )
49	REAL szz( iip1,jjp1,llm ),zz
50	INTEGER i,j,l,indice
51	real sxn(iip1),sxs(iip1)
52
53	real sinlon(iip1),sinlondlon(iip1)
54	real coslon(iip1),coslondlon(iip1)
55	real qmin,qmax
56	save qmin,qmax
57	save sinlon,coslon,sinlondlon,coslondlon
58	real dyn1,dyn2,dys1,dys2,qpn,qps,dqzpn,dqzps
59	real masn,mass
60	c
61	REAL SSUM
62	integer ismax,ismin
63	EXTERNAL SSUM, convflu,ismin,ismax
64	logical first
65	save first
66	EXTERNAL advxp,advyp,advzp
67
68
69	data first/.true./
70	data qmin,qmax/-1.e33,1.e33/
71
72
73	c==========================================================================
74	c==========================================================================
75	c MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
76	c==========================================================================
77	c==========================================================================
78	REAL dt
79	c==========================================================================
80	limit = .TRUE.
81
82	if(first) then
83	print*,'SCHEMA PRATHER'
84	first=.false.
85	do i=2,iip1
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)
95
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
111	endif
112	c Fin modif Fred
113
114	c *** On calcule la masse d'air en kg
115
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
123
124	c *** q contient les qqtes de traceur avant l'advection
125
126	c *** Affectation des tableaux S a partir de Q
127
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	c *** Appel des subroutines d'advection en X, en Y et en Z
145	c *** Advection avec "time-splitting"
146
147	c-----------------------------------------------------------
148	do indice =1,nt
149	call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz
150	. ,sxx,sxy,sxz,syy,syz,szz,1 )
151	end do
152	do l=1,llm
153	do i=1,iip1
154	sy(i,1,l)=0.
155	sy(i,jjp1,l)=0.
156	enddo
157	enddo
158	c---------------------------------------------------------
159	call advyp( limit,.5dtnt,pbarv,sm,s0,sx,sy,sz
160	. ,sxx,sxy,sxz,syy,syz,szz,1 )
161	c---------------------------------------------------------
162
163	c---------------------------------------------------------
164	do j=1,jjp1
165	do i=1,iip1
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.
174	enddo
175	enddo
176	call advzp( limit,dt*nt,w,sm,s0,sx,sy,sz
177	. ,sxx,sxy,sxz,syy,syz,szz,1 )
178	do l=1,llm
179	do i=1,iip1
180	sy(i,1,l)=0.
181	sy(i,jjp1,l)=0.
182	enddo
183	enddo
184
185	c---------------------------------------------------------
186
187	c---------------------------------------------------------
188	call advyp( limit,.5dtnt,pbarv,sm,s0,sx,sy,sz
189	. ,sxx,sxy,sxz,syy,syz,szz,1 )
190	c---------------------------------------------------------
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
205	do indice=1,nt
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	c---------------------------------------------------------
210	c---------------------------------------------------------
211	c *** On repasse les S dans la variable qpr
212	c *** 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	c---------------------------------------------------------
232	c go to 777
233	c filtrages aux poles
234
235	c Traitements specifiques au pole
236
237	c filtrages aux poles
238	DO l=1,llm
239	c 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
246	do i=1,iip1
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	c enddo
253	c print*,'qpn',qpn,'qps',qps
254	c print*,'dqzpn',dqzpn,'dqzps',dqzps
255	c enddo
256	dyn1=0.
257	dys1=0.
258	dyn2=0.
259	dys2=0.
260	do i=1,iim
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
268	do i=1,iim
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)
280	do i=1,iim
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)
286	do i=1,iim
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
294	do l=1,llm
295	do i=1,iim
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.
308	enddo
309	ENDDO
310
311	777 continue
312	c
313	c bouclage en longitude
314	do l=1,llm
315	do j=1,jjp1
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).lt.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)
340	c PRINT*,' PBL EN SORTIE D'' ADVZP'
341	q(i,j,l,0)=0.
342	c STOP
343	ENDIF
344	ENDDO
345	ENDDO
346	do j=1,jjp1,jjm
347	do i=1,iip1
348	IF (q(i,j,l,0).lt.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	c STOP
357	ENDIF
358	enddo
359	enddo
360	ENDDO
361	RETURN
362	END

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

Download in other formats:

Original Format