/[lmdze]/trunk/dyn3d/Vlsplt/vlx.f
ViewVC logotype

Annotation of /trunk/dyn3d/Vlsplt/vlx.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 31 - (hide annotations)
Thu Apr 1 14:59:19 2010 UTC (14 years, 2 months ago) by guez
Original Path: trunk/libf/dyn3d/Vlsplt/vlx.f90
File size: 8391 byte(s) 
Split "vlsplt.f" in single-procedure files. Gathered the files in
directory "dyn3d/Vlsplt".

Defined "pbarum(:, 1, :)" and "pbarum(:, jjm + 1, :)" in procedure
"groupe".
1 guez 31 SUBROUTINE vlx(q,pente_max,masse,u_m)
2    
3     ! Auteurs: P.Le Van, F.Hourdin, F.Forget
4     !
5     ! *******************************************************
6     ! Shema d'advection " pseudo amont " .
7     ! ************************************************************
8     ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le
9     ! s-pg ....
10     !
11     !
12     ! --------------------------------------------------------------
13     use dimens_m
14     use paramet_m
15     use comconst
16     use comvert
17     use logic
18     IMPLICIT NONE
19     !
20     !
21     !
22     ! Arguments:
23     ! ----------
24     REAL masse(ip1jmp1,llm),pente_max
25     REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
26     REAL q(ip1jmp1,llm)
27     REAL w(ip1jmp1,llm)
28     !
29     ! Local
30     ! ---------
31     !
32     INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
33     INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
34     !
35     REAL new_m,zu_m,zdum(ip1jmp1,llm)
36     REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
37     REAL zz(ip1jmp1)
38     REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
39     REAL u_mq(ip1jmp1,llm)
40    
41     Logical extremum,first,testcpu
42     SAVE first,testcpu
43    
44     REAL SSUM
45     REAL temps0,temps1,temps2,temps3,temps4,temps5,second
46     SAVE temps0,temps1,temps2,temps3,temps4,temps5
47    
48     REAL z1,z2,z3
49    
50     DATA first,testcpu/.true.,.false./
51    
52     IF(first) THEN
53     temps1=0.
54     temps2=0.
55     temps3=0.
56     temps4=0.
57     temps5=0.
58     first=.false.
59     ENDIF
60    
61     ! calcul de la pente a droite et a gauche de la maille
62    
63    
64     IF (pente_max.gt.-1.e-5) THEN
65     ! IF (pente_max.gt.10) THEN
66    
67     ! calcul des pentes avec limitation, Van Leer scheme I:
68     ! -----------------------------------------------------
69    
70     ! calcul de la pente aux points u
71     DO l = 1, llm
72     DO ij=iip2,ip1jm-1
73     dxqu(ij)=q(ij+1,l)-q(ij,l)
74     ! IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
75     ! sigu(ij)=u_m(ij,l)/masse(ij,l)
76     ENDDO
77     DO ij=iip1+iip1,ip1jm,iip1
78     dxqu(ij)=dxqu(ij-iim)
79     ! sigu(ij)=sigu(ij-iim)
80     ENDDO
81    
82     DO ij=iip2,ip1jm
83     adxqu(ij)=abs(dxqu(ij))
84     ENDDO
85    
86     ! calcul de la pente maximum dans la maille en valeur absolue
87    
88     DO ij=iip2+1,ip1jm
89     dxqmax(ij,l)=pente_max* &
90     & min(adxqu(ij-1),adxqu(ij))
91     ! limitation subtile
92     ! , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
93    
94    
95     ENDDO
96    
97     DO ij=iip1+iip1,ip1jm,iip1
98     dxqmax(ij-iim,l)=dxqmax(ij,l)
99     ENDDO
100    
101     DO ij=iip2+1,ip1jm
102     IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
103     dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
104     ELSE
105     ! extremum local
106     dxq(ij,l)=0.
107     ENDIF
108     dxq(ij,l)=0.5*dxq(ij,l)
109     dxq(ij,l)= &
110     & sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
111     ENDDO
112    
113     ! l=1,llm
114     ENDDO
115     !print*,'Ok calcul des pentes'
116    
117     ! (pente_max.lt.-1.e-5)
118     ELSE
119    
120     ! Pentes produits:
121     ! ----------------
122    
123     DO l = 1, llm
124     DO ij=iip2,ip1jm-1
125     dxqu(ij)=q(ij+1,l)-q(ij,l)
126     ENDDO
127     DO ij=iip1+iip1,ip1jm,iip1
128     dxqu(ij)=dxqu(ij-iim)
129     ENDDO
130    
131     DO ij=iip2+1,ip1jm
132     zz(ij)=dxqu(ij-1)*dxqu(ij)
133     zz(ij)=zz(ij)+zz(ij)
134     IF(zz(ij).gt.0) THEN
135     dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
136     ELSE
137     ! extremum local
138     dxq(ij,l)=0.
139     ENDIF
140     ENDDO
141    
142     ENDDO
143    
144     ! (pente_max.lt.-1.e-5)
145     ENDIF
146    
147     ! bouclage de la pente en iip1:
148     ! -----------------------------
149    
150     DO l=1,llm
151     DO ij=iip1+iip1,ip1jm,iip1
152     dxq(ij-iim,l)=dxq(ij,l)
153     ENDDO
154     DO ij=1,ip1jmp1
155     iadvplus(ij,l)=0
156     ENDDO
157    
158     ENDDO
159    
160     ! print*,'Bouclage en iip1'
161    
162     ! calcul des flux a gauche et a droite
163    
164     ! on cumule le flux correspondant a toutes les mailles dont la masse
165     ! au travers de la paroi pENDant le pas de temps.
166     !print*,'Cumule ....'
167    
168     DO l=1,llm
169     DO ij=iip2,ip1jm-1
170     ! print*,'masse(',ij,')=',masse(ij,l)
171     IF (u_m(ij,l).gt.0.) THEN
172     zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
173     u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
174     ELSE
175     zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
176     u_mq(ij,l)=u_m(ij,l) &
177     & *(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
178     ENDIF
179     ENDDO
180     ENDDO
181     ! stop
182    
183     ! go to 9999
184     ! detection des points ou on advecte plus que la masse de la
185     ! maille
186     DO l=1,llm
187     DO ij=iip2,ip1jm-1
188     IF(zdum(ij,l).lt.0) THEN
189     iadvplus(ij,l)=1
190     u_mq(ij,l)=0.
191     ENDIF
192     ENDDO
193     ENDDO
194     !print*,'Ok test 1'
195     DO l=1,llm
196     DO ij=iip1+iip1,ip1jm,iip1
197     iadvplus(ij,l)=iadvplus(ij-iim,l)
198     ENDDO
199     ENDDO
200     ! print*,'Ok test 2'
201    
202    
203     ! traitement special pour le cas ou on advecte en longitude plus
204     ! que le
205     ! contenu de la maille.
206     ! cette partie est mal vectorisee.
207    
208     ! calcul du nombre de maille sur lequel on advecte plus que la maille.
209    
210     n0=0
211     DO l=1,llm
212     nl(l)=0
213     DO ij=iip2,ip1jm
214     nl(l)=nl(l)+iadvplus(ij,l)
215     ENDDO
216     n0=n0+nl(l)
217     ENDDO
218    
219     IF(n0.gt.0) THEN
220     !C PRINT*,'Nombre de points pour lesquels on advect plus que le'
221     !C & ,'contenu de la maille : ',n0
222    
223     DO l=1,llm
224     IF(nl(l).gt.0) THEN
225     iju=0
226     ! indicage des mailles concernees par le traitement special
227     DO ij=iip2,ip1jm
228     IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
229     iju=iju+1
230     indu(iju)=ij
231     ENDIF
232     ENDDO
233     niju=iju
234     ! PRINT*,'niju,nl',niju,nl(l)
235    
236     ! traitement des mailles
237     DO iju=1,niju
238     ij=indu(iju)
239     j=(ij-1)/iip1+1
240     zu_m=u_m(ij,l)
241     u_mq(ij,l)=0.
242     IF(zu_m.gt.0.) THEN
243     ijq=ij
244     i=ijq-(j-1)*iip1
245     ! accumulation pour les mailles completements advectees
246     do while(zu_m.gt.masse(ijq,l))
247     u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
248     zu_m=zu_m-masse(ijq,l)
249     i=mod(i-2+iim,iim)+1
250     ijq=(j-1)*iip1+i
251     ENDDO
252     ! ajout de la maille non completement advectee
253     u_mq(ij,l)=u_mq(ij,l)+zu_m* &
254     & (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
255     ELSE
256     ijq=ij+1
257     i=ijq-(j-1)*iip1
258     ! accumulation pour les mailles completements advectees
259     do while(-zu_m.gt.masse(ijq,l))
260     u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
261     zu_m=zu_m+masse(ijq,l)
262     i=mod(i,iim)+1
263     ijq=(j-1)*iip1+i
264     ENDDO
265     ! ajout de la maille non completement advectee
266     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)- &
267     & 0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
268     ENDIF
269     ENDDO
270     ENDIF
271     ENDDO
272     ! n0.gt.0
273     ENDIF
274     9999 continue
275    
276    
277     ! bouclage en latitude
278     !print*,'cvant bouclage en latitude'
279     DO l=1,llm
280     DO ij=iip1+iip1,ip1jm,iip1
281     u_mq(ij,l)=u_mq(ij-iim,l)
282     ENDDO
283     ENDDO
284    
285    
286     ! calcul des tENDances
287    
288     DO l=1,llm
289     DO ij=iip2+1,ip1jm
290     new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
291     q(ij,l)=(q(ij,l)*masse(ij,l)+ &
292     & u_mq(ij-1,l)-u_mq(ij,l)) &
293     & /new_m
294     masse(ij,l)=new_m
295     ENDDO
296     ! ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
297     DO ij=iip1+iip1,ip1jm,iip1
298     q(ij-iim,l)=q(ij,l)
299     masse(ij-iim,l)=masse(ij,l)
300     ENDDO
301     ENDDO
302    
303    
304     RETURN
305     END
  ViewVC Help
Powered by ViewVC 1.1.21