Context Navigation

← Previous Revision
Next Revision →
Blame
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

Line
1
2	! $Header$
3
4	SUBROUTINE prather(q,w,masse,pbaru,pbarv,nt,dt)
5
6	USE comconst_mod, ONLY: pi
7	USE lmdz_ssum_scopy, ONLY: ssum
8	USE lmdz_comgeom2
9
10	USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
11	USE lmdz_paramet
12	IMPLICIT NONE
13
14	!=======================================================================
15	! Adaptation LMDZ: A.Armengaud (LGGE)
16	! ----------------
17
18	! ************************************************
19	! Transport des traceurs par la methode de prather
20	! Ref :
21
22	! ************************************************
23	! q,w,pext,pbaru et pbarv : arguments d'entree pour le s-pg
24
25	!=======================================================================
26
27
28
29
30
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 )
38	INTEGER :: ordre,ilim
39
40	! Local:
41	! ------
42	LOGICAL :: limit
43	REAL :: zq(iip1,jjp1,llm)
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
54	REAL :: sxn(iip1),sxs(iip1)
55
56	REAL :: sinlon(iip1),sinlondlon(iip1)
57	REAL :: coslon(iip1),coslondlon(iip1)
58	REAL :: qmin,qmax
59	save qmin,qmax
60	save sinlon,coslon,sinlondlon,coslondlon
61	REAL :: dyn1,dyn2,dys1,dys2,qpn,qps,dqzpn,dqzps
62	REAL :: masn,mass
63
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	!==========================================================================
74	!==========================================================================
75	! MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
76	!==========================================================================
77	!==========================================================================
78	REAL :: dt
79	!==========================================================================
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	! Fin modif Fred
113
114	! *** 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	! *** q contient les qqtes de traceur avant l'advection
125
126	! *** 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	! *** Appel des subroutines d'advection en X, en Y et en Z
145	! *** Advection avec "time-splitting"
146
147	!-----------------------------------------------------------
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	!---------------------------------------------------------
159	CALL advyp( limit,.5dtnt,pbarv,sm,s0,sx,sy,sz &
160	,sxx,sxy,sxz,syy,syz,szz,1 )
161	!---------------------------------------------------------
162
163	!---------------------------------------------------------
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	!---------------------------------------------------------
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
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	!---------------------------------------------------------
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
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	! 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.
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
313	! 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)<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	! PRINT*,' PBL EN SORTIE D'' ADVZP'
341	q(i,j,l,0)=0.
342	! 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)<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