Context Navigation

prather.F @ 2597

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

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

Download in other formats:

Original Format