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

Annotation of /trunk/dyn3d/pentes_ini.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3 - (hide annotations)
Wed Feb 27 13:16:39 2008 UTC (16 years, 3 months ago) by guez
Original Path: trunk/libf/dyn3d/pentes_ini.f
File size: 13825 byte(s)
Initial import
1 guez 3 !
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