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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 32 - (show annotations)
Tue Apr 6 17:52:58 2010 UTC (14 years, 1 month ago) by guez
Original Path: trunk/libf/dyn3d/Vlsplt/vlx.f90
File size: 8330 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 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 ! (pente_max.lt.-1.e-5)
116 ELSE
117
118 ! Pentes produits:
119 ! ----------------
120
121 DO l = 1, llm
122 DO ij=iip2,ip1jm-1
123 dxqu(ij)=q(ij+1,l)-q(ij,l)
124 ENDDO
125 DO ij=iip1+iip1,ip1jm,iip1
126 dxqu(ij)=dxqu(ij-iim)
127 ENDDO
128
129 DO ij=iip2+1,ip1jm
130 zz(ij)=dxqu(ij-1)*dxqu(ij)
131 zz(ij)=zz(ij)+zz(ij)
132 IF(zz(ij).gt.0) THEN
133 dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
134 ELSE
135 ! extremum local
136 dxq(ij,l)=0.
137 ENDIF
138 ENDDO
139
140 ENDDO
141
142 ! (pente_max.lt.-1.e-5)
143 ENDIF
144
145 ! bouclage de la pente en iip1:
146 ! -----------------------------
147
148 DO l=1,llm
149 DO ij=iip1+iip1,ip1jm,iip1
150 dxq(ij-iim,l)=dxq(ij,l)
151 ENDDO
152 DO ij=1,ip1jmp1
153 iadvplus(ij,l)=0
154 ENDDO
155
156 ENDDO
157
158 ! calcul des flux a gauche et a droite
159
160 ! on cumule le flux correspondant a toutes les mailles dont la masse
161 ! au travers de la paroi pENDant le pas de temps.
162 !print*,'Cumule ....'
163
164 DO l=1,llm
165 DO ij=iip2,ip1jm-1
166 ! print*,'masse(',ij,')=',masse(ij,l)
167 IF (u_m(ij,l).gt.0.) THEN
168 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
169 u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
170 ELSE
171 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
172 u_mq(ij,l)=u_m(ij,l) &
173 & *(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
174 ENDIF
175 ENDDO
176 ENDDO
177 ! stop
178
179 ! go to 9999
180 ! detection des points ou on advecte plus que la masse de la
181 ! maille
182 DO l=1,llm
183 DO ij=iip2,ip1jm-1
184 IF(zdum(ij,l).lt.0) THEN
185 iadvplus(ij,l)=1
186 u_mq(ij,l)=0.
187 ENDIF
188 ENDDO
189 ENDDO
190 !print*,'Ok test 1'
191 DO l=1,llm
192 DO ij=iip1+iip1,ip1jm,iip1
193 iadvplus(ij,l)=iadvplus(ij-iim,l)
194 ENDDO
195 ENDDO
196 ! print*,'Ok test 2'
197
198
199 ! traitement special pour le cas ou on advecte en longitude plus
200 ! que le
201 ! contenu de la maille.
202 ! cette partie est mal vectorisee.
203
204 ! calcul du nombre de maille sur lequel on advecte plus que la maille.
205
206 n0=0
207 DO l=1,llm
208 nl(l)=0
209 DO ij=iip2,ip1jm
210 nl(l)=nl(l)+iadvplus(ij,l)
211 ENDDO
212 n0=n0+nl(l)
213 ENDDO
214
215 IF(n0.gt.0) THEN
216 !C PRINT*,'Nombre de points pour lesquels on advect plus que le'
217 !C & ,'contenu de la maille : ',n0
218
219 DO l=1,llm
220 IF(nl(l).gt.0) THEN
221 iju=0
222 ! indicage des mailles concernees par le traitement special
223 DO ij=iip2,ip1jm
224 IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
225 iju=iju+1
226 indu(iju)=ij
227 ENDIF
228 ENDDO
229 niju=iju
230 ! PRINT*,'niju,nl',niju,nl(l)
231
232 ! traitement des mailles
233 DO iju=1,niju
234 ij=indu(iju)
235 j=(ij-1)/iip1+1
236 zu_m=u_m(ij,l)
237 u_mq(ij,l)=0.
238 IF(zu_m.gt.0.) THEN
239 ijq=ij
240 i=ijq-(j-1)*iip1
241 ! accumulation pour les mailles completements advectees
242 do while(zu_m.gt.masse(ijq,l))
243 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
244 zu_m=zu_m-masse(ijq,l)
245 i=mod(i-2+iim,iim)+1
246 ijq=(j-1)*iip1+i
247 ENDDO
248 ! ajout de la maille non completement advectee
249 u_mq(ij,l)=u_mq(ij,l)+zu_m* &
250 & (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
251 ELSE
252 ijq=ij+1
253 i=ijq-(j-1)*iip1
254 ! accumulation pour les mailles completements advectees
255 do while(-zu_m.gt.masse(ijq,l))
256 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
257 zu_m=zu_m+masse(ijq,l)
258 i=mod(i,iim)+1
259 ijq=(j-1)*iip1+i
260 ENDDO
261 ! ajout de la maille non completement advectee
262 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)- &
263 & 0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
264 ENDIF
265 ENDDO
266 ENDIF
267 ENDDO
268 ! n0.gt.0
269 ENDIF
270 9999 continue
271
272
273 ! bouclage en latitude
274 !print*,'cvant bouclage en latitude'
275 DO l=1,llm
276 DO ij=iip1+iip1,ip1jm,iip1
277 u_mq(ij,l)=u_mq(ij-iim,l)
278 ENDDO
279 ENDDO
280
281
282 ! calcul des tENDances
283
284 DO l=1,llm
285 DO ij=iip2+1,ip1jm
286 new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
287 q(ij,l)=(q(ij,l)*masse(ij,l)+ &
288 & u_mq(ij-1,l)-u_mq(ij,l)) &
289 & /new_m
290 masse(ij,l)=new_m
291 ENDDO
292 ! ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
293 DO ij=iip1+iip1,ip1jm,iip1
294 q(ij-iim,l)=q(ij,l)
295 masse(ij-iim,l)=masse(ij,l)
296 ENDDO
297 ENDDO
298
299
300 RETURN
301 END
  ViewVC Help
Powered by ViewVC 1.1.21