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, 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,.5*dt*nt,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,.5*dt*nt,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