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

Annotation of /trunk/dyn3d/pentes_ini.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 32 - (hide annotations)
Tue Apr 6 17:52:58 2010 UTC (14 years, 2 months ago) by guez
Original Path: trunk/libf/dyn3d/pentes_ini.f
File size: 11995 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 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 guez 31 REAL, intent(in):: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
34 guez 3 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    
113     c *** q contient les qqtes de traceur avant l'advection
114    
115     c *** Affectation des tableaux S a partir de Q
116     c *** Rem : utilisation de SCOPY ulterieurement
117    
118     DO l = 1,llm
119     DO j = 1,jjp1
120     DO i = 1,iip1
121     s0( i,j,llm+1-l ) = q ( i,j,l,0 )
122     sx( i,j,llm+1-l ) = q ( i,j,l,1 )
123     sy( i,j,llm+1-l ) = q ( i,j,l,2 )
124     sz( i,j,llm+1-l ) = q ( i,j,l,3 )
125     ENDDO
126     ENDDO
127     ENDDO
128    
129     c *** On calcule la masse d'air en kg
130    
131     DO l = 1,llm
132     DO j = 1,jjp1
133     DO i = 1,iip1
134     sm ( i,j,llm+1-l)=masse( i,j,l )
135     ENDDO
136     ENDDO
137     ENDDO
138    
139     c *** On converti les champs S en atome (resp. kg)
140     c *** Les routines d'advection traitent les champs
141     c *** a advecter si ces derniers sont en atome (resp. kg)
142     c *** A optimiser !!!
143    
144     DO l = 1,llm
145     DO j = 1,jjp1
146     DO i = 1,iip1
147     s0(i,j,l) = s0(i,j,l) * sm ( i,j,l )
148     sx(i,j,l) = sx(i,j,l) * sm ( i,j,l )
149     sy(i,j,l) = sy(i,j,l) * sm ( i,j,l )
150     sz(i,j,l) = sz(i,j,l) * sm ( i,j,l )
151     ENDDO
152     ENDDO
153     ENDDO
154    
155     c *** Appel des subroutines d'advection en X, en Y et en Z
156     c *** Advection avec "time-splitting"
157    
158     if(mode.eq.2) then
159     do l=1,llm
160     s0s=0.
161     s0n=0.
162     dyn1=0.
163     dys1=0.
164     dyn2=0.
165     dys2=0.
166     smn=0.
167     sms=0.
168     do i=1,iim
169     smn=smn+sm(i,1,l)
170     sms=sms+sm(i,jjp1,l)
171     s0n=s0n+s0(i,1,l)
172     s0s=s0s+s0(i,jjp1,l)
173     zz=sy(i,1,l)/sm(i,1,l)
174     dyn1=dyn1+sinlondlon(i)*zz
175     dyn2=dyn2+coslondlon(i)*zz
176     zz=sy(i,jjp1,l)/sm(i,jjp1,l)
177     dys1=dys1+sinlondlon(i)*zz
178     dys2=dys2+coslondlon(i)*zz
179     enddo
180     do i=1,iim
181     sy(i,1,l)=dyn1*sinlon(i)+dyn2*coslon(i)
182     sy(i,jjp1,l)=dys1*sinlon(i)+dys2*coslon(i)
183     enddo
184     do i=1,iim
185     s0(i,1,l)=s0n/smn+sy(i,1,l)
186     s0(i,jjp1,l)=s0s/sms-sy(i,jjp1,l)
187     enddo
188    
189     s0(iip1,1,l)=s0(1,1,l)
190     s0(iip1,jjp1,l)=s0(1,jjp1,l)
191    
192     do i=1,iim
193     sxn(i)=s0(i+1,1,l)-s0(i,1,l)
194     sxs(i)=s0(i+1,jjp1,l)-s0(i,jjp1,l)
195     c on rerentre les masses
196     enddo
197     do i=1,iim
198     sy(i,1,l)=sy(i,1,l)*sm(i,1,l)
199     sy(i,jjp1,l)=sy(i,jjp1,l)*sm(i,jjp1,l)
200     s0(i,1,l)=s0(i,1,l)*sm(i,1,l)
201     s0(i,jjp1,l)=s0(i,jjp1,l)*sm(i,jjp1,l)
202     enddo
203     sxn(iip1)=sxn(1)
204     sxs(iip1)=sxs(1)
205     do i=1,iim
206     sx(i+1,1,l)=0.25*(sxn(i)+sxn(i+1))*sm(i+1,1,l)
207     sx(i+1,jjp1,l)=0.25*(sxs(i)+sxs(i+1))*sm(i+1,jjp1,l)
208     enddo
209     s0(iip1,1,l)=s0(1,1,l)
210     s0(iip1,jjp1,l)=s0(1,jjp1,l)
211     sy(iip1,1,l)=sy(1,1,l)
212     sy(iip1,jjp1,l)=sy(1,jjp1,l)
213     sx(1,1,l)=sx(iip1,1,l)
214     sx(1,jjp1,l)=sx(iip1,jjp1,l)
215     enddo
216     endif
217    
218     if (mode.eq.4) then
219     do l=1,llm
220     do i=1,iip1
221     sx(i,1,l)=0.
222     sx(i,jjp1,l)=0.
223     sy(i,1,l)=0.
224     sy(i,jjp1,l)=0.
225     enddo
226     enddo
227     endif
228     call limx(s0,sx,sm,pente_max)
229     call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
230     if (mode.eq.4) then
231     do l=1,llm
232     do i=1,iip1
233     sx(i,1,l)=0.
234     sx(i,jjp1,l)=0.
235     sy(i,1,l)=0.
236     sy(i,jjp1,l)=0.
237     enddo
238     enddo
239     endif
240     call limy(s0,sy,sm,pente_max)
241     call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )
242     do j=1,jjp1
243     do i=1,iip1
244     sz(i,j,1)=0.
245     sz(i,j,llm)=0.
246     enddo
247     enddo
248     call limz(s0,sz,sm,pente_max)
249     call advz( limit,dtvr,w,sm,s0,sx,sy,sz )
250     if (mode.eq.4) then
251     do l=1,llm
252     do i=1,iip1
253     sx(i,1,l)=0.
254     sx(i,jjp1,l)=0.
255     sy(i,1,l)=0.
256     sy(i,jjp1,l)=0.
257     enddo
258     enddo
259     endif
260     call limy(s0,sy,sm,pente_max)
261     call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )
262     do l=1,llm
263     do j=1,jjp1
264     sm(iip1,j,l)=sm(1,j,l)
265     s0(iip1,j,l)=s0(1,j,l)
266     sx(iip1,j,l)=sx(1,j,l)
267     sy(iip1,j,l)=sy(1,j,l)
268     sz(iip1,j,l)=sz(1,j,l)
269     enddo
270     enddo
271    
272    
273     if (mode.eq.4) then
274     do l=1,llm
275     do i=1,iip1
276     sx(i,1,l)=0.
277     sx(i,jjp1,l)=0.
278     sy(i,1,l)=0.
279     sy(i,jjp1,l)=0.
280     enddo
281     enddo
282     endif
283     call limx(s0,sx,sm,pente_max)
284     call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
285     c *** On repasse les S dans la variable q directement 14/10/94
286     c On revient a des rapports de melange en divisant par la masse
287    
288     c En dehors des poles:
289    
290     DO l = 1,llm
291     DO j = 1,jjp1
292     DO i = 1,iim
293     q(i,j,llm+1-l,0)=s0(i,j,l)/sm(i,j,l)
294     q(i,j,llm+1-l,1)=sx(i,j,l)/sm(i,j,l)
295     q(i,j,llm+1-l,2)=sy(i,j,l)/sm(i,j,l)
296     q(i,j,llm+1-l,3)=sz(i,j,l)/sm(i,j,l)
297     ENDDO
298     ENDDO
299     ENDDO
300    
301     c Traitements specifiques au pole
302    
303     if(mode.ge.1) then
304     DO l=1,llm
305     c filtrages aux poles
306     masn=ssum(iim,sm(1,1,l),1)
307     mass=ssum(iim,sm(1,jjp1,l),1)
308     qpn=ssum(iim,s0(1,1,l),1)/masn
309     qps=ssum(iim,s0(1,jjp1,l),1)/mass
310     dqzpn=ssum(iim,sz(1,1,l),1)/masn
311     dqzps=ssum(iim,sz(1,jjp1,l),1)/mass
312     do i=1,iip1
313     q( i,1,llm+1-l,3)=dqzpn
314     q( i,jjp1,llm+1-l,3)=dqzps
315     q( i,1,llm+1-l,0)=qpn
316     q( i,jjp1,llm+1-l,0)=qps
317     enddo
318     if(mode.eq.3) then
319     dyn1=0.
320     dys1=0.
321     dyn2=0.
322     dys2=0.
323     do i=1,iim
324     dyn1=dyn1+sinlondlon(i)*sy(i,1,l)/sm(i,1,l)
325     dyn2=dyn2+coslondlon(i)*sy(i,1,l)/sm(i,1,l)
326     dys1=dys1+sinlondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)
327     dys2=dys2+coslondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)
328     enddo
329     do i=1,iim
330     q(i,1,llm+1-l,2)=
331     s (sinlon(i)*dyn1+coslon(i)*dyn2)
332     q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)
333     q(i,jjp1,llm+1-l,2)=
334     s (sinlon(i)*dys1+coslon(i)*dys2)
335     q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
336     s -q(i,jjp1,llm+1-l,2)
337     enddo
338     endif
339     if(mode.eq.1) then
340     c on filtre les valeurs au bord de la "grande maille pole"
341     dyn1=0.
342     dys1=0.
343     dyn2=0.
344     dys2=0.
345     do i=1,iim
346     zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)
347     dyn1=dyn1+sinlondlon(i)*zz
348     dyn2=dyn2+coslondlon(i)*zz
349     zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)
350     dys1=dys1+sinlondlon(i)*zz
351     dys2=dys2+coslondlon(i)*zz
352     enddo
353     do i=1,iim
354     q(i,1,llm+1-l,2)=
355     s (sinlon(i)*dyn1+coslon(i)*dyn2)/2.
356     q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)
357     q(i,jjp1,llm+1-l,2)=
358     s (sinlon(i)*dys1+coslon(i)*dys2)/2.
359     q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
360     s -q(i,jjp1,llm+1-l,2)
361     enddo
362     q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)
363     q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)
364    
365     do i=1,iim
366     sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)
367     sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)
368     enddo
369     sxn(iip1)=sxn(1)
370     sxs(iip1)=sxs(1)
371     do i=1,iim
372     q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))
373     q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))
374     enddo
375     q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)
376     q(1,jjp1,llm+1-l,1)=q(iip1,jjp1,llm+1-l,1)
377    
378     endif
379    
380     ENDDO
381     endif
382    
383     c bouclage en longitude
384     do iq=0,3
385     do l=1,llm
386     do j=1,jjp1
387     q(iip1,j,l,iq)=q(1,j,l,iq)
388     enddo
389     enddo
390     enddo
391    
392     DO l = 1,llm
393     DO j = 1,jjp1
394     DO i = 1,iip1
395     IF (q(i,j,l,0).lt.0.) THEN
396     q(i,j,l,0)=0.
397     ENDIF
398     ENDDO
399     ENDDO
400     ENDDO
401    
402     do l=1,llm
403     do j=1,jjp1
404     do i=1,iip1
405     if(q(i,j,l,0).lt.qmin)
406     , print*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0)
407     enddo
408     enddo
409     enddo
410     RETURN
411     END

  ViewVC Help
Powered by ViewVC 1.1.21