/[lmdze]/trunk/libf/dyn3d/prather.f
ViewVC logotype

Contents of /trunk/libf/dyn3d/prather.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 32 - (show annotations)
Tue Apr 6 17:52:58 2010 UTC (14 years, 1 month ago) by guez
File size: 10800 byte(s) 
Split "stringop.f90" into single-procedure files. Gathered files in directory
"IOIPSL/Stringop".

Split "flincom.f90" into "flincom.f90" and "flinget.f90". Removed
unused procedures from module "flincom". Removed unused argument
"filename" of procedure "flinopen_nozoom".

Removed unused files.

Split "grid_change.f90" into "grid_change.f90" and
"gr_phy_write_3d.f90".

Removed unused procedures from modules "calendar", "ioipslmpp",
"grid_atob", "gath_cpl" and "getincom". Removed unused procedures in
files "ppm3d.f" and "thermcell.f".

Split "mathelp.f90" into "mathelp.f90" and "mathop.f90".

Removed unused variable "dpres" of module "comvert".

Use argument "itau" instead of local variables "iadvtr" and "first" to
control algorithm in procedure "fluxstokenc".

Removed unused arguments of procedure "integrd".

Removed useless computations at the end of procedure "leapfrog".

Merged common block "matrfil" into module "parafilt".
1 !
2 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/prather.F,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $
3 !
4 SUBROUTINE prather (q,w,masse,pbaru,pbarv,nt,dt)
5 use dimens_m
6 use paramet_m
7 use comconst
8 use comvert
9 use comgeom
10 IMPLICIT NONE
11
12 c=======================================================================
13 c Adaptation LMDZ: A.Armengaud (LGGE)
14 c ----------------
15 c
16 c ************************************************
17 c Transport des traceurs par la methode de prather
18 c Ref :
19 c
20 c ************************************************
21 c q,w,pext,pbaru et pbarv : arguments d'entree pour le s-pg
22 c
23 c=======================================================================
24
25
26
27 c Arguments:
28 c ----------
29 INTEGER iq,nt
30 REAL, intent(in):: 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 dyn1=0.
252 dys1=0.
253 dyn2=0.
254 dys2=0.
255 do i=1,iim
256 zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)
257 dyn1=dyn1+sinlondlon(i)*zz
258 dyn2=dyn2+coslondlon(i)*zz
259 zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)
260 dys1=dys1+sinlondlon(i)*zz
261 dys2=dys2+coslondlon(i)*zz
262 enddo
263 do i=1,iim
264 q(i,1,llm+1-l,2)=
265 $ (sinlon(i)*dyn1+coslon(i)*dyn2)/2.
266 q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)
267 $ +q(i,1,llm+1-l,2)
268 q(i,jjp1,llm+1-l,2)=
269 $ (sinlon(i)*dys1+coslon(i)*dys2)/2.
270 q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
271 $ -q(i,jjp1,llm+1-l,2)
272 enddo
273 q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)
274 q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)
275 do i=1,iim
276 sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)
277 sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)
278 enddo
279 sxn(iip1)=sxn(1)
280 sxs(iip1)=sxs(1)
281 do i=1,iim
282 q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))
283 q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))
284 END DO
285 q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)
286 q(1,jjp1,llm+1-l,1)=
287 $ q(iip1,jjp1,llm+1-l,1)
288 enddo
289 do l=1,llm
290 do i=1,iim
291 q( i,1,llm+1-l,4)=0.
292 q( i,jjp1,llm+1-l,4)=0.
293 q( i,1,llm+1-l,5)=0.
294 q( i,jjp1,llm+1-l,5)=0.
295 q( i,1,llm+1-l,6)=0.
296 q( i,jjp1,llm+1-l,6)=0.
297 q( i,1,llm+1-l,7)=0.
298 q( i,jjp1,llm+1-l,7)=0.
299 q( i,1,llm+1-l,8)=0.
300 q( i,jjp1,llm+1-l,8)=0.
301 q( i,1,llm+1-l,9)=0.
302 q( i,jjp1,llm+1-l,9)=0.
303 enddo
304 ENDDO
305
306 777 continue
307 c
308 c bouclage en longitude
309 do l=1,llm
310 do j=1,jjp1
311 q(iip1,j,l,0)=q(1,j,l,0)
312 q(iip1,j,llm+1-l,0)=q(1,j,llm+1-l,0)
313 q(iip1,j,llm+1-l,1)=q(1,j,llm+1-l,1)
314 q(iip1,j,llm+1-l,2)=q(1,j,llm+1-l,2)
315 q(iip1,j,llm+1-l,3)=q(1,j,llm+1-l,3)
316 q(iip1,j,llm+1-l,4)=q(1,j,llm+1-l,4)
317 q(iip1,j,llm+1-l,5)=q(1,j,llm+1-l,5)
318 q(iip1,j,llm+1-l,6)=q(1,j,llm+1-l,6)
319 q(iip1,j,llm+1-l,7)=q(1,j,llm+1-l,7)
320 q(iip1,j,llm+1-l,8)=q(1,j,llm+1-l,8)
321 q(iip1,j,llm+1-l,9)=q(1,j,llm+1-l,9)
322 enddo
323 enddo
324 DO l = 1,llm
325 DO j = 2,jjm
326 DO i = 1,iip1
327 IF (q(i,j,l,0).lt.0.) THEN
328 PRINT*,'------------ BIP-----------'
329 PRINT*,'S0(',i,j,l,')=',q(i,j,l,0),
330 $ q(i,j-1,l,0)
331 PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)
332 PRINT*,'SY(',i,j,l,')=',q(i,j,l,2),
333 $ q(i,j-1,l,2)
334 PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)
335 q(i,j,l,0)=0.
336 ENDIF
337 ENDDO
338 ENDDO
339 do j=1,jjp1,jjm
340 do i=1,iip1
341 IF (q(i,j,l,0).lt.0.) THEN
342 PRINT*,'------------ BIP 2-----------'
343 PRINT*,'S0(',i,j,l,')=',q(i,j,l,0)
344 PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)
345 PRINT*,'SY(',i,j,l,')=',q(i,j,l,2)
346 PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)
347
348 q(i,j,l,0)=0.
349 c STOP
350 ENDIF
351 enddo
352 enddo
353 ENDDO
354 RETURN
355 END
  ViewVC Help
Powered by ViewVC 1.1.21