My Project
 All Classes Files Functions Variables Macros
prather.F
Go to the documentation of this file.
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,.5*dt*nt,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,.5*dt*nt,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