/[lmdze]/trunk/libf/dyn3d/Vlsplt/vlx.f90
ViewVC logotype

Contents of /trunk/libf/dyn3d/Vlsplt/vlx.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 31 - (show annotations)
Thu Apr 1 14:59:19 2010 UTC (14 years, 1 month ago) by guez
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 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