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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 39 - (hide annotations)
Tue Jan 25 15:11:05 2011 UTC (13 years, 4 months ago) by guez
File size: 12024 byte(s)
"pi" comes from "nr_util". Removed subroutine "initialize" in module
"comconst".

Copied the content of "fxy_sin.h" into "fxysinus", instead of getting
it from an "include" line. Removed file "fxy_sin.h".

"ps" has rank 2 in "gcm" and "dynetat0".

Assumed-shape for argument "q" of "integrd".

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 guez 39 USE nr_util, ONLY : pi
11 guez 3 IMPLICIT NONE
12    
13     c=======================================================================
14     c Adaptation LMDZ: A.Armengaud (LGGE)
15     c ----------------
16     c
17     c ********************************************************************
18     c Transport des traceurs par la methode des pentes
19     c ********************************************************************
20     c Reference possible : Russel. G.L., Lerner J.A.:
21     c A new Finite-Differencing Scheme for Traceur Transport
22     c Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81
23     c ********************************************************************
24     c q,w,masse,pbaru et pbarv
25     c sont des arguments d'entree pour le s-pg ....
26     c
27     c=======================================================================
28    
29    
30    
31     c Arguments:
32     c ----------
33     integer mode
34 guez 31 REAL, intent(in):: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
35 guez 3 REAL q( iip1,jjp1,llm,0:3)
36     REAL w( ip1jmp1,llm )
37     REAL masse( iip1,jjp1,llm)
38     c Local:
39     c ------
40     LOGICAL limit
41     REAL sm ( iip1,jjp1, llm )
42     REAL s0( iip1,jjp1,llm ), sx( iip1,jjp1,llm )
43     REAL sy( iip1,jjp1,llm ), sz( iip1,jjp1,llm )
44     real masn,mass,zz
45     INTEGER i,j,l,iq
46    
47     c modif Fred 24 03 96
48    
49     real sinlon(iip1),sinlondlon(iip1)
50     real coslon(iip1),coslondlon(iip1)
51     save sinlon,coslon,sinlondlon,coslondlon
52     real dyn1,dyn2,dys1,dys2
53     real qpn,qps,dqzpn,dqzps
54     real smn,sms,s0n,s0s,sxn(iip1),sxs(iip1)
55     real qmin,zq,pente_max
56     c
57     REAL SSUM
58     integer ismax,ismin,lati,latf
59     EXTERNAL SSUM, convflu,ismin,ismax
60     logical first
61     save first
62     c fin modif
63    
64     c EXTERNAL masskg
65     EXTERNAL advx
66     EXTERNAL advy
67     EXTERNAL advz
68    
69     c modif Fred 24 03 96
70     data first/.true./
71    
72     limit = .TRUE.
73     pente_max=2
74     c if (mode.eq.1.or.mode.eq.3) then
75     c if (mode.eq.1) then
76     if (mode.ge.1) then
77     lati=2
78     latf=jjm
79     else
80     lati=1
81     latf=jjp1
82     endif
83    
84     qmin=0.4995
85     qmin=0.
86     if(first) then
87     print*,'SCHEMA AMONT NOUVEAU'
88     first=.false.
89     do i=2,iip1
90     coslon(i)=cos(rlonv(i))
91     sinlon(i)=sin(rlonv(i))
92     coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
93     sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
94     print*,coslondlon(i),sinlondlon(i)
95     enddo
96     coslon(1)=coslon(iip1)
97     coslondlon(1)=coslondlon(iip1)
98     sinlon(1)=sinlon(iip1)
99     sinlondlon(1)=sinlondlon(iip1)
100     print*,'sum sinlondlon ',ssum(iim,sinlondlon,1)/sinlondlon(1)
101     print*,'sum coslondlon ',ssum(iim,coslondlon,1)/coslondlon(1)
102     DO l = 1,llm
103     DO j = 1,jjp1
104     DO i = 1,iip1
105     q ( i,j,l,1 )=0.
106     q ( i,j,l,2 )=0.
107     q ( i,j,l,3 )=0.
108     ENDDO
109     ENDDO
110     ENDDO
111    
112     endif
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 *** On calcule la masse d'air en kg
131    
132     DO l = 1,llm
133     DO j = 1,jjp1
134     DO i = 1,iip1
135     sm ( i,j,llm+1-l)=masse( i,j,l )
136     ENDDO
137     ENDDO
138     ENDDO
139    
140     c *** On converti les champs S en atome (resp. kg)
141     c *** Les routines d'advection traitent les champs
142     c *** a advecter si ces derniers sont en atome (resp. kg)
143     c *** A optimiser !!!
144    
145     DO l = 1,llm
146     DO j = 1,jjp1
147     DO i = 1,iip1
148     s0(i,j,l) = s0(i,j,l) * sm ( i,j,l )
149     sx(i,j,l) = sx(i,j,l) * sm ( i,j,l )
150     sy(i,j,l) = sy(i,j,l) * sm ( i,j,l )
151     sz(i,j,l) = sz(i,j,l) * sm ( i,j,l )
152     ENDDO
153     ENDDO
154     ENDDO
155    
156     c *** Appel des subroutines d'advection en X, en Y et en Z
157     c *** Advection avec "time-splitting"
158    
159     if(mode.eq.2) then
160     do l=1,llm
161     s0s=0.
162     s0n=0.
163     dyn1=0.
164     dys1=0.
165     dyn2=0.
166     dys2=0.
167     smn=0.
168     sms=0.
169     do i=1,iim
170     smn=smn+sm(i,1,l)
171     sms=sms+sm(i,jjp1,l)
172     s0n=s0n+s0(i,1,l)
173     s0s=s0s+s0(i,jjp1,l)
174     zz=sy(i,1,l)/sm(i,1,l)
175     dyn1=dyn1+sinlondlon(i)*zz
176     dyn2=dyn2+coslondlon(i)*zz
177     zz=sy(i,jjp1,l)/sm(i,jjp1,l)
178     dys1=dys1+sinlondlon(i)*zz
179     dys2=dys2+coslondlon(i)*zz
180     enddo
181     do i=1,iim
182     sy(i,1,l)=dyn1*sinlon(i)+dyn2*coslon(i)
183     sy(i,jjp1,l)=dys1*sinlon(i)+dys2*coslon(i)
184     enddo
185     do i=1,iim
186     s0(i,1,l)=s0n/smn+sy(i,1,l)
187     s0(i,jjp1,l)=s0s/sms-sy(i,jjp1,l)
188     enddo
189    
190     s0(iip1,1,l)=s0(1,1,l)
191     s0(iip1,jjp1,l)=s0(1,jjp1,l)
192    
193     do i=1,iim
194     sxn(i)=s0(i+1,1,l)-s0(i,1,l)
195     sxs(i)=s0(i+1,jjp1,l)-s0(i,jjp1,l)
196     c on rerentre les masses
197     enddo
198     do i=1,iim
199     sy(i,1,l)=sy(i,1,l)*sm(i,1,l)
200     sy(i,jjp1,l)=sy(i,jjp1,l)*sm(i,jjp1,l)
201     s0(i,1,l)=s0(i,1,l)*sm(i,1,l)
202     s0(i,jjp1,l)=s0(i,jjp1,l)*sm(i,jjp1,l)
203     enddo
204     sxn(iip1)=sxn(1)
205     sxs(iip1)=sxs(1)
206     do i=1,iim
207     sx(i+1,1,l)=0.25*(sxn(i)+sxn(i+1))*sm(i+1,1,l)
208     sx(i+1,jjp1,l)=0.25*(sxs(i)+sxs(i+1))*sm(i+1,jjp1,l)
209     enddo
210     s0(iip1,1,l)=s0(1,1,l)
211     s0(iip1,jjp1,l)=s0(1,jjp1,l)
212     sy(iip1,1,l)=sy(1,1,l)
213     sy(iip1,jjp1,l)=sy(1,jjp1,l)
214     sx(1,1,l)=sx(iip1,1,l)
215     sx(1,jjp1,l)=sx(iip1,jjp1,l)
216     enddo
217     endif
218    
219     if (mode.eq.4) then
220     do l=1,llm
221     do i=1,iip1
222     sx(i,1,l)=0.
223     sx(i,jjp1,l)=0.
224     sy(i,1,l)=0.
225     sy(i,jjp1,l)=0.
226     enddo
227     enddo
228     endif
229     call limx(s0,sx,sm,pente_max)
230     call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
231     if (mode.eq.4) then
232     do l=1,llm
233     do i=1,iip1
234     sx(i,1,l)=0.
235     sx(i,jjp1,l)=0.
236     sy(i,1,l)=0.
237     sy(i,jjp1,l)=0.
238     enddo
239     enddo
240     endif
241     call limy(s0,sy,sm,pente_max)
242     call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )
243     do j=1,jjp1
244     do i=1,iip1
245     sz(i,j,1)=0.
246     sz(i,j,llm)=0.
247     enddo
248     enddo
249     call limz(s0,sz,sm,pente_max)
250     call advz( limit,dtvr,w,sm,s0,sx,sy,sz )
251     if (mode.eq.4) then
252     do l=1,llm
253     do i=1,iip1
254     sx(i,1,l)=0.
255     sx(i,jjp1,l)=0.
256     sy(i,1,l)=0.
257     sy(i,jjp1,l)=0.
258     enddo
259     enddo
260     endif
261     call limy(s0,sy,sm,pente_max)
262     call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )
263     do l=1,llm
264     do j=1,jjp1
265     sm(iip1,j,l)=sm(1,j,l)
266     s0(iip1,j,l)=s0(1,j,l)
267     sx(iip1,j,l)=sx(1,j,l)
268     sy(iip1,j,l)=sy(1,j,l)
269     sz(iip1,j,l)=sz(1,j,l)
270     enddo
271     enddo
272    
273    
274     if (mode.eq.4) then
275     do l=1,llm
276     do i=1,iip1
277     sx(i,1,l)=0.
278     sx(i,jjp1,l)=0.
279     sy(i,1,l)=0.
280     sy(i,jjp1,l)=0.
281     enddo
282     enddo
283     endif
284     call limx(s0,sx,sm,pente_max)
285     call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
286     c *** On repasse les S dans la variable q directement 14/10/94
287     c On revient a des rapports de melange en divisant par la masse
288    
289     c En dehors des poles:
290    
291     DO l = 1,llm
292     DO j = 1,jjp1
293     DO i = 1,iim
294     q(i,j,llm+1-l,0)=s0(i,j,l)/sm(i,j,l)
295     q(i,j,llm+1-l,1)=sx(i,j,l)/sm(i,j,l)
296     q(i,j,llm+1-l,2)=sy(i,j,l)/sm(i,j,l)
297     q(i,j,llm+1-l,3)=sz(i,j,l)/sm(i,j,l)
298     ENDDO
299     ENDDO
300     ENDDO
301    
302     c Traitements specifiques au pole
303    
304     if(mode.ge.1) then
305     DO l=1,llm
306     c filtrages aux poles
307     masn=ssum(iim,sm(1,1,l),1)
308     mass=ssum(iim,sm(1,jjp1,l),1)
309     qpn=ssum(iim,s0(1,1,l),1)/masn
310     qps=ssum(iim,s0(1,jjp1,l),1)/mass
311     dqzpn=ssum(iim,sz(1,1,l),1)/masn
312     dqzps=ssum(iim,sz(1,jjp1,l),1)/mass
313     do i=1,iip1
314     q( i,1,llm+1-l,3)=dqzpn
315     q( i,jjp1,llm+1-l,3)=dqzps
316     q( i,1,llm+1-l,0)=qpn
317     q( i,jjp1,llm+1-l,0)=qps
318     enddo
319     if(mode.eq.3) then
320     dyn1=0.
321     dys1=0.
322     dyn2=0.
323     dys2=0.
324     do i=1,iim
325     dyn1=dyn1+sinlondlon(i)*sy(i,1,l)/sm(i,1,l)
326     dyn2=dyn2+coslondlon(i)*sy(i,1,l)/sm(i,1,l)
327     dys1=dys1+sinlondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)
328     dys2=dys2+coslondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)
329     enddo
330     do i=1,iim
331     q(i,1,llm+1-l,2)=
332     s (sinlon(i)*dyn1+coslon(i)*dyn2)
333     q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)
334     q(i,jjp1,llm+1-l,2)=
335     s (sinlon(i)*dys1+coslon(i)*dys2)
336     q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
337     s -q(i,jjp1,llm+1-l,2)
338     enddo
339     endif
340     if(mode.eq.1) then
341     c on filtre les valeurs au bord de la "grande maille pole"
342     dyn1=0.
343     dys1=0.
344     dyn2=0.
345     dys2=0.
346     do i=1,iim
347     zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)
348     dyn1=dyn1+sinlondlon(i)*zz
349     dyn2=dyn2+coslondlon(i)*zz
350     zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)
351     dys1=dys1+sinlondlon(i)*zz
352     dys2=dys2+coslondlon(i)*zz
353     enddo
354     do i=1,iim
355     q(i,1,llm+1-l,2)=
356     s (sinlon(i)*dyn1+coslon(i)*dyn2)/2.
357     q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)
358     q(i,jjp1,llm+1-l,2)=
359     s (sinlon(i)*dys1+coslon(i)*dys2)/2.
360     q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
361     s -q(i,jjp1,llm+1-l,2)
362     enddo
363     q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)
364     q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)
365    
366     do i=1,iim
367     sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)
368     sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)
369     enddo
370     sxn(iip1)=sxn(1)
371     sxs(iip1)=sxs(1)
372     do i=1,iim
373     q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))
374     q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))
375     enddo
376     q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)
377     q(1,jjp1,llm+1-l,1)=q(iip1,jjp1,llm+1-l,1)
378    
379     endif
380    
381     ENDDO
382     endif
383    
384     c bouclage en longitude
385     do iq=0,3
386     do l=1,llm
387     do j=1,jjp1
388     q(iip1,j,l,iq)=q(1,j,l,iq)
389     enddo
390     enddo
391     enddo
392    
393     DO l = 1,llm
394     DO j = 1,jjp1
395     DO i = 1,iip1
396     IF (q(i,j,l,0).lt.0.) THEN
397     q(i,j,l,0)=0.
398     ENDIF
399     ENDDO
400     ENDDO
401     ENDDO
402    
403     do l=1,llm
404     do j=1,jjp1
405     do i=1,iip1
406     if(q(i,j,l,0).lt.qmin)
407     , print*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0)
408     enddo
409     enddo
410     enddo
411     RETURN
412     END

  ViewVC Help
Powered by ViewVC 1.1.21