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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 32 - (hide 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 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     ! (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