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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3 - (show annotations)
Wed Feb 27 13:16:39 2008 UTC (16 years, 3 months ago) by guez
File size: 13825 byte(s)
Initial import
1 !
2 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/pentes_ini.F,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $
3 !
4 SUBROUTINE pentes_ini (q,w,masse,pbaru,pbarv,mode)
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 des pentes
18 c ********************************************************************
19 c Reference possible : Russel. G.L., Lerner J.A.:
20 c A new Finite-Differencing Scheme for Traceur Transport
21 c Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81
22 c ********************************************************************
23 c q,w,masse,pbaru et pbarv
24 c sont des arguments d'entree pour le s-pg ....
25 c
26 c=======================================================================
27
28
29
30 c Arguments:
31 c ----------
32 integer mode
33 REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
34 REAL q( iip1,jjp1,llm,0:3)
35 REAL w( ip1jmp1,llm )
36 REAL masse( iip1,jjp1,llm)
37 c Local:
38 c ------
39 LOGICAL limit
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 masn,mass,zz
44 INTEGER i,j,l,iq
45
46 c modif Fred 24 03 96
47
48 real sinlon(iip1),sinlondlon(iip1)
49 real coslon(iip1),coslondlon(iip1)
50 save sinlon,coslon,sinlondlon,coslondlon
51 real dyn1,dyn2,dys1,dys2
52 real qpn,qps,dqzpn,dqzps
53 real smn,sms,s0n,s0s,sxn(iip1),sxs(iip1)
54 real qmin,zq,pente_max
55 c
56 REAL SSUM
57 integer ismax,ismin,lati,latf
58 EXTERNAL SSUM, convflu,ismin,ismax
59 logical first
60 save first
61 c fin modif
62
63 c EXTERNAL masskg
64 EXTERNAL advx
65 EXTERNAL advy
66 EXTERNAL advz
67
68 c modif Fred 24 03 96
69 data first/.true./
70
71 limit = .TRUE.
72 pente_max=2
73 c if (mode.eq.1.or.mode.eq.3) then
74 c if (mode.eq.1) then
75 if (mode.ge.1) then
76 lati=2
77 latf=jjm
78 else
79 lati=1
80 latf=jjp1
81 endif
82
83 qmin=0.4995
84 qmin=0.
85 if(first) then
86 print*,'SCHEMA AMONT NOUVEAU'
87 first=.false.
88 do i=2,iip1
89 coslon(i)=cos(rlonv(i))
90 sinlon(i)=sin(rlonv(i))
91 coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
92 sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
93 print*,coslondlon(i),sinlondlon(i)
94 enddo
95 coslon(1)=coslon(iip1)
96 coslondlon(1)=coslondlon(iip1)
97 sinlon(1)=sinlon(iip1)
98 sinlondlon(1)=sinlondlon(iip1)
99 print*,'sum sinlondlon ',ssum(iim,sinlondlon,1)/sinlondlon(1)
100 print*,'sum coslondlon ',ssum(iim,coslondlon,1)/coslondlon(1)
101 DO l = 1,llm
102 DO j = 1,jjp1
103 DO i = 1,iip1
104 q ( i,j,l,1 )=0.
105 q ( i,j,l,2 )=0.
106 q ( i,j,l,3 )=0.
107 ENDDO
108 ENDDO
109 ENDDO
110
111 endif
112 c Fin modif Fred
113
114 c *** q contient les qqtes de traceur avant l'advection
115
116 c *** Affectation des tableaux S a partir de Q
117 c *** Rem : utilisation de SCOPY ulterieurement
118
119 DO l = 1,llm
120 DO j = 1,jjp1
121 DO i = 1,iip1
122 s0( i,j,llm+1-l ) = q ( i,j,l,0 )
123 sx( i,j,llm+1-l ) = q ( i,j,l,1 )
124 sy( i,j,llm+1-l ) = q ( i,j,l,2 )
125 sz( i,j,llm+1-l ) = q ( i,j,l,3 )
126 ENDDO
127 ENDDO
128 ENDDO
129
130 c PRINT*,'----- S0 just before conversion -------'
131 c PRINT*,'S0(16,12,1)=',s0(16,12,1)
132 c PRINT*,'Q(16,12,1,4)=',q(16,12,1,4)
133
134 c *** On calcule la masse d'air en kg
135
136 DO l = 1,llm
137 DO j = 1,jjp1
138 DO i = 1,iip1
139 sm ( i,j,llm+1-l)=masse( i,j,l )
140 ENDDO
141 ENDDO
142 ENDDO
143
144 c *** On converti les champs S en atome (resp. kg)
145 c *** Les routines d'advection traitent les champs
146 c *** a advecter si ces derniers sont en atome (resp. kg)
147 c *** A optimiser !!!
148
149 DO l = 1,llm
150 DO j = 1,jjp1
151 DO i = 1,iip1
152 s0(i,j,l) = s0(i,j,l) * sm ( i,j,l )
153 sx(i,j,l) = sx(i,j,l) * sm ( i,j,l )
154 sy(i,j,l) = sy(i,j,l) * sm ( i,j,l )
155 sz(i,j,l) = sz(i,j,l) * sm ( i,j,l )
156 ENDDO
157 ENDDO
158 ENDDO
159
160 c ss0 = 0.
161 c DO l = 1,llm
162 c DO j = 1,jjp1
163 c DO i = 1,iim
164 c ss0 = ss0 + s0 ( i,j,l )
165 c ENDDO
166 c ENDDO
167 c ENDDO
168 c PRINT*, 'valeur tot s0 avant advection=',ss0
169
170 c *** Appel des subroutines d'advection en X, en Y et en Z
171 c *** Advection avec "time-splitting"
172
173 c-----------------------------------------------------------
174 c PRINT*,'----- S0 just before ADVX -------'
175 c PRINT*,'S0(16,12,1)=',s0(16,12,1)
176
177 c-----------------------------------------------------------
178 c do l=1,llm
179 c do j=1,jjp1
180 c do i=1,iip1
181 c zq=s0(i,j,l)/sm(i,j,l)
182 c if(zq.lt.qmin)
183 c , print*,'avant advx1, s0(',i,',',j,',',l,')=',zq
184 c enddo
185 c enddo
186 c enddo
187 CCC
188 if(mode.eq.2) then
189 do l=1,llm
190 s0s=0.
191 s0n=0.
192 dyn1=0.
193 dys1=0.
194 dyn2=0.
195 dys2=0.
196 smn=0.
197 sms=0.
198 do i=1,iim
199 smn=smn+sm(i,1,l)
200 sms=sms+sm(i,jjp1,l)
201 s0n=s0n+s0(i,1,l)
202 s0s=s0s+s0(i,jjp1,l)
203 zz=sy(i,1,l)/sm(i,1,l)
204 dyn1=dyn1+sinlondlon(i)*zz
205 dyn2=dyn2+coslondlon(i)*zz
206 zz=sy(i,jjp1,l)/sm(i,jjp1,l)
207 dys1=dys1+sinlondlon(i)*zz
208 dys2=dys2+coslondlon(i)*zz
209 enddo
210 do i=1,iim
211 sy(i,1,l)=dyn1*sinlon(i)+dyn2*coslon(i)
212 sy(i,jjp1,l)=dys1*sinlon(i)+dys2*coslon(i)
213 enddo
214 do i=1,iim
215 s0(i,1,l)=s0n/smn+sy(i,1,l)
216 s0(i,jjp1,l)=s0s/sms-sy(i,jjp1,l)
217 enddo
218
219 s0(iip1,1,l)=s0(1,1,l)
220 s0(iip1,jjp1,l)=s0(1,jjp1,l)
221
222 do i=1,iim
223 sxn(i)=s0(i+1,1,l)-s0(i,1,l)
224 sxs(i)=s0(i+1,jjp1,l)-s0(i,jjp1,l)
225 c on rerentre les masses
226 enddo
227 do i=1,iim
228 sy(i,1,l)=sy(i,1,l)*sm(i,1,l)
229 sy(i,jjp1,l)=sy(i,jjp1,l)*sm(i,jjp1,l)
230 s0(i,1,l)=s0(i,1,l)*sm(i,1,l)
231 s0(i,jjp1,l)=s0(i,jjp1,l)*sm(i,jjp1,l)
232 enddo
233 sxn(iip1)=sxn(1)
234 sxs(iip1)=sxs(1)
235 do i=1,iim
236 sx(i+1,1,l)=0.25*(sxn(i)+sxn(i+1))*sm(i+1,1,l)
237 sx(i+1,jjp1,l)=0.25*(sxs(i)+sxs(i+1))*sm(i+1,jjp1,l)
238 enddo
239 s0(iip1,1,l)=s0(1,1,l)
240 s0(iip1,jjp1,l)=s0(1,jjp1,l)
241 sy(iip1,1,l)=sy(1,1,l)
242 sy(iip1,jjp1,l)=sy(1,jjp1,l)
243 sx(1,1,l)=sx(iip1,1,l)
244 sx(1,jjp1,l)=sx(iip1,jjp1,l)
245 enddo
246 endif
247
248 if (mode.eq.4) then
249 do l=1,llm
250 do i=1,iip1
251 sx(i,1,l)=0.
252 sx(i,jjp1,l)=0.
253 sy(i,1,l)=0.
254 sy(i,jjp1,l)=0.
255 enddo
256 enddo
257 endif
258 call limx(s0,sx,sm,pente_max)
259 c call minmaxq(zq,1.e33,-1.e33,'avant advx ')
260 call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
261 c call minmaxq(zq,1.e33,-1.e33,'avant advy ')
262 if (mode.eq.4) then
263 do l=1,llm
264 do i=1,iip1
265 sx(i,1,l)=0.
266 sx(i,jjp1,l)=0.
267 sy(i,1,l)=0.
268 sy(i,jjp1,l)=0.
269 enddo
270 enddo
271 endif
272 call limy(s0,sy,sm,pente_max)
273 call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )
274 c call minmaxq(zq,1.e33,-1.e33,'avant advz ')
275 do j=1,jjp1
276 do i=1,iip1
277 sz(i,j,1)=0.
278 sz(i,j,llm)=0.
279 enddo
280 enddo
281 call limz(s0,sz,sm,pente_max)
282 call advz( limit,dtvr,w,sm,s0,sx,sy,sz )
283 if (mode.eq.4) then
284 do l=1,llm
285 do i=1,iip1
286 sx(i,1,l)=0.
287 sx(i,jjp1,l)=0.
288 sy(i,1,l)=0.
289 sy(i,jjp1,l)=0.
290 enddo
291 enddo
292 endif
293 call limy(s0,sy,sm,pente_max)
294 call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )
295 do l=1,llm
296 do j=1,jjp1
297 sm(iip1,j,l)=sm(1,j,l)
298 s0(iip1,j,l)=s0(1,j,l)
299 sx(iip1,j,l)=sx(1,j,l)
300 sy(iip1,j,l)=sy(1,j,l)
301 sz(iip1,j,l)=sz(1,j,l)
302 enddo
303 enddo
304
305
306 c call minmaxq(zq,1.e33,-1.e33,'avant advx ')
307 if (mode.eq.4) then
308 do l=1,llm
309 do i=1,iip1
310 sx(i,1,l)=0.
311 sx(i,jjp1,l)=0.
312 sy(i,1,l)=0.
313 sy(i,jjp1,l)=0.
314 enddo
315 enddo
316 endif
317 call limx(s0,sx,sm,pente_max)
318 call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
319 c call minmaxq(zq,1.e33,-1.e33,'apres advx ')
320 c do l=1,llm
321 c do j=1,jjp1
322 c do i=1,iip1
323 c zq=s0(i,j,l)/sm(i,j,l)
324 c if(zq.lt.qmin)
325 c , print*,'apres advx2, s0(',i,',',j,',',l,')=',zq
326 c enddo
327 c enddo
328 c enddo
329 c *** On repasse les S dans la variable q directement 14/10/94
330 c On revient a des rapports de melange en divisant par la masse
331
332 c En dehors des poles:
333
334 DO l = 1,llm
335 DO j = 1,jjp1
336 DO i = 1,iim
337 q(i,j,llm+1-l,0)=s0(i,j,l)/sm(i,j,l)
338 q(i,j,llm+1-l,1)=sx(i,j,l)/sm(i,j,l)
339 q(i,j,llm+1-l,2)=sy(i,j,l)/sm(i,j,l)
340 q(i,j,llm+1-l,3)=sz(i,j,l)/sm(i,j,l)
341 ENDDO
342 ENDDO
343 ENDDO
344
345 c Traitements specifiques au pole
346
347 if(mode.ge.1) then
348 DO l=1,llm
349 c filtrages aux poles
350 masn=ssum(iim,sm(1,1,l),1)
351 mass=ssum(iim,sm(1,jjp1,l),1)
352 qpn=ssum(iim,s0(1,1,l),1)/masn
353 qps=ssum(iim,s0(1,jjp1,l),1)/mass
354 dqzpn=ssum(iim,sz(1,1,l),1)/masn
355 dqzps=ssum(iim,sz(1,jjp1,l),1)/mass
356 do i=1,iip1
357 q( i,1,llm+1-l,3)=dqzpn
358 q( i,jjp1,llm+1-l,3)=dqzps
359 q( i,1,llm+1-l,0)=qpn
360 q( i,jjp1,llm+1-l,0)=qps
361 enddo
362 if(mode.eq.3) then
363 dyn1=0.
364 dys1=0.
365 dyn2=0.
366 dys2=0.
367 do i=1,iim
368 dyn1=dyn1+sinlondlon(i)*sy(i,1,l)/sm(i,1,l)
369 dyn2=dyn2+coslondlon(i)*sy(i,1,l)/sm(i,1,l)
370 dys1=dys1+sinlondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)
371 dys2=dys2+coslondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)
372 enddo
373 do i=1,iim
374 q(i,1,llm+1-l,2)=
375 s (sinlon(i)*dyn1+coslon(i)*dyn2)
376 q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)
377 q(i,jjp1,llm+1-l,2)=
378 s (sinlon(i)*dys1+coslon(i)*dys2)
379 q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
380 s -q(i,jjp1,llm+1-l,2)
381 enddo
382 endif
383 if(mode.eq.1) then
384 c on filtre les valeurs au bord de la "grande maille pole"
385 dyn1=0.
386 dys1=0.
387 dyn2=0.
388 dys2=0.
389 do i=1,iim
390 zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)
391 dyn1=dyn1+sinlondlon(i)*zz
392 dyn2=dyn2+coslondlon(i)*zz
393 zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)
394 dys1=dys1+sinlondlon(i)*zz
395 dys2=dys2+coslondlon(i)*zz
396 enddo
397 do i=1,iim
398 q(i,1,llm+1-l,2)=
399 s (sinlon(i)*dyn1+coslon(i)*dyn2)/2.
400 q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)
401 q(i,jjp1,llm+1-l,2)=
402 s (sinlon(i)*dys1+coslon(i)*dys2)/2.
403 q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
404 s -q(i,jjp1,llm+1-l,2)
405 enddo
406 q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)
407 q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)
408
409 do i=1,iim
410 sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)
411 sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)
412 enddo
413 sxn(iip1)=sxn(1)
414 sxs(iip1)=sxs(1)
415 do i=1,iim
416 q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))
417 q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))
418 enddo
419 q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)
420 q(1,jjp1,llm+1-l,1)=q(iip1,jjp1,llm+1-l,1)
421
422 endif
423
424 ENDDO
425 endif
426
427 c bouclage en longitude
428 do iq=0,3
429 do l=1,llm
430 do j=1,jjp1
431 q(iip1,j,l,iq)=q(1,j,l,iq)
432 enddo
433 enddo
434 enddo
435
436 c PRINT*, ' SORTIE DE PENTES --- ca peut glisser ....'
437
438 DO l = 1,llm
439 DO j = 1,jjp1
440 DO i = 1,iip1
441 IF (q(i,j,l,0).lt.0.) THEN
442 c PRINT*,'------------ BIP-----------'
443 c PRINT*,'Q0(',i,j,l,')=',q(i,j,l,0)
444 c PRINT*,'QX(',i,j,l,')=',q(i,j,l,1)
445 c PRINT*,'QY(',i,j,l,')=',q(i,j,l,2)
446 c PRINT*,'QZ(',i,j,l,')=',q(i,j,l,3)
447 c PRINT*,' PBL EN SORTIE DE PENTES'
448 q(i,j,l,0)=0.
449 c STOP
450 ENDIF
451 ENDDO
452 ENDDO
453 ENDDO
454
455 c PRINT*, '-------------------------------------------'
456
457 do l=1,llm
458 do j=1,jjp1
459 do i=1,iip1
460 if(q(i,j,l,0).lt.qmin)
461 , print*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0)
462 enddo
463 enddo
464 enddo
465 RETURN
466 END
467
468
469
470
471
472
473
474
475
476
477
478

  ViewVC Help
Powered by ViewVC 1.1.21