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

Annotation of /trunk/dyn3d/vlxqs.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 254 - (hide annotations)
Mon Feb 5 10:39:38 2018 UTC (6 years, 3 months ago) by guez
File size: 7550 byte(s) 
Move Sources/* to root directory.
1 guez 43 SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat)
2     !
3     ! Auteurs: P.Le Van, F.Hourdin, F.Forget
4     !
5     ! ********************************************************************
6     ! Shema d'advection " pseudo amont " .
7     ! ********************************************************************
8     !
9     ! --------------------------------------------------------------------
10     use dimens_m
11     use paramet_m
12     use comconst
13 guez 66 use disvert_m
14 guez 57 use conf_gcm_m
15 guez 43 IMPLICIT NONE
16     !
17     !
18     !
19     ! Arguments:
20     ! ----------
21     REAL masse(ip1jmp1,llm),pente_max
22     REAL u_m( ip1jmp1,llm )
23     REAL q(ip1jmp1,llm)
24     REAL qsat(ip1jmp1,llm)
25     !
26     ! Local
27     ! ---------
28     !
29     INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
30     INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
31     !
32     REAL new_m,zu_m,zdum(ip1jmp1,llm)
33     REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
34     REAL zz(ip1jmp1)
35     REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
36     REAL u_mq(ip1jmp1,llm)
37    
38     ! calcul de la pente a droite et a gauche de la maille
39    
40    
41     IF (pente_max.gt.-1.e-5) THEN
42     ! IF (pente_max.gt.10) THEN
43    
44     ! calcul des pentes avec limitation, Van Leer scheme I:
45     ! -----------------------------------------------------
46    
47     ! calcul de la pente aux points u
48     DO l = 1, llm
49     DO ij=iip2,ip1jm-1
50     dxqu(ij)=q(ij+1,l)-q(ij,l)
51     ! IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
52     ! sigu(ij)=u_m(ij,l)/masse(ij,l)
53     ENDDO
54     DO ij=iip1+iip1,ip1jm,iip1
55     dxqu(ij)=dxqu(ij-iim)
56     ! sigu(ij)=sigu(ij-iim)
57     ENDDO
58    
59     DO ij=iip2,ip1jm
60     adxqu(ij)=abs(dxqu(ij))
61     ENDDO
62    
63     ! calcul de la pente maximum dans la maille en valeur absolue
64    
65     DO ij=iip2+1,ip1jm
66     dxqmax(ij,l)=pente_max* &
67     min(adxqu(ij-1),adxqu(ij))
68     ! limitation subtile
69     ! , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
70    
71    
72     ENDDO
73    
74     DO ij=iip1+iip1,ip1jm,iip1
75     dxqmax(ij-iim,l)=dxqmax(ij,l)
76     ENDDO
77    
78     DO ij=iip2+1,ip1jm
79     IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
80     dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
81     ELSE
82     ! extremum local
83     dxq(ij,l)=0.
84     ENDIF
85     dxq(ij,l)=0.5*dxq(ij,l)
86     dxq(ij,l)= &
87     sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
88     ENDDO
89    
90     ENDDO ! l=1,llm
91    
92     ELSE ! (pente_max.lt.-1.e-5)
93    
94     ! Pentes produits:
95     ! ----------------
96    
97     DO l = 1, llm
98     DO ij=iip2,ip1jm-1
99     dxqu(ij)=q(ij+1,l)-q(ij,l)
100     ENDDO
101     DO ij=iip1+iip1,ip1jm,iip1
102     dxqu(ij)=dxqu(ij-iim)
103     ENDDO
104    
105     DO ij=iip2+1,ip1jm
106     zz(ij)=dxqu(ij-1)*dxqu(ij)
107     zz(ij)=zz(ij)+zz(ij)
108     IF(zz(ij).gt.0) THEN
109     dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
110     ELSE
111     ! extremum local
112     dxq(ij,l)=0.
113     ENDIF
114     ENDDO
115    
116     ENDDO
117    
118     ENDIF ! (pente_max.lt.-1.e-5)
119    
120     ! bouclage de la pente en iip1:
121     ! -----------------------------
122    
123     DO l=1,llm
124     DO ij=iip1+iip1,ip1jm,iip1
125     dxq(ij-iim,l)=dxq(ij,l)
126     ENDDO
127    
128     DO ij=1,ip1jmp1
129     iadvplus(ij,l)=0
130     ENDDO
131    
132     ENDDO
133    
134    
135     ! calcul des flux a gauche et a droite
136    
137     ! on cumule le flux correspondant a toutes les mailles dont la masse
138     ! au travers de la paroi pENDant le pas de temps.
139     ! le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
140     DO l=1,llm
141     DO ij=iip2,ip1jm-1
142     IF (u_m(ij,l).gt.0.) THEN
143     zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
144     u_mq(ij,l)=u_m(ij,l)* &
145     min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
146     ELSE
147     zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
148     u_mq(ij,l)=u_m(ij,l)* &
149     min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
150     ENDIF
151     ENDDO
152     ENDDO
153    
154    
155     ! detection des points ou on advecte plus que la masse de la
156     ! maille
157     DO l=1,llm
158     DO ij=iip2,ip1jm-1
159     IF(zdum(ij,l).lt.0) THEN
160     iadvplus(ij,l)=1
161     u_mq(ij,l)=0.
162     ENDIF
163     ENDDO
164     ENDDO
165     DO l=1,llm
166     DO ij=iip1+iip1,ip1jm,iip1
167     iadvplus(ij,l)=iadvplus(ij-iim,l)
168     ENDDO
169     ENDDO
170    
171    
172    
173     ! traitement special pour le cas ou on advecte en longitude plus que le
174     ! contenu de la maille.
175     ! cette partie est mal vectorisee.
176    
177     ! pas d'influence de la pression saturante (pour l'instant)
178    
179     ! calcul du nombre de maille sur lequel on advecte plus que la maille.
180    
181     n0=0
182     DO l=1,llm
183     nl(l)=0
184     DO ij=iip2,ip1jm
185     nl(l)=nl(l)+iadvplus(ij,l)
186     ENDDO
187     n0=n0+nl(l)
188     ENDDO
189    
190     IF(n0.gt.0) THEN
191     DO l=1,llm
192     IF(nl(l).gt.0) THEN
193     iju=0
194     ! indicage des mailles concernees par le traitement special
195     DO ij=iip2,ip1jm
196     IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
197     iju=iju+1
198     indu(iju)=ij
199     ENDIF
200     ENDDO
201     niju=iju
202    
203     ! traitement des mailles
204     DO iju=1,niju
205     ij=indu(iju)
206     j=(ij-1)/iip1+1
207     zu_m=u_m(ij,l)
208     u_mq(ij,l)=0.
209     IF(zu_m.gt.0.) THEN
210     ijq=ij
211     i=ijq-(j-1)*iip1
212     ! accumulation pour les mailles completements advectees
213     do while(zu_m.gt.masse(ijq,l))
214     u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
215     zu_m=zu_m-masse(ijq,l)
216     i=mod(i-2+iim,iim)+1
217     ijq=(j-1)*iip1+i
218     ENDDO
219     ! ajout de la maille non completement advectee
220     u_mq(ij,l)=u_mq(ij,l)+zu_m* &
221     (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
222     ELSE
223     ijq=ij+1
224     i=ijq-(j-1)*iip1
225     ! accumulation pour les mailles completements advectees
226     do while(-zu_m.gt.masse(ijq,l))
227     u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
228     zu_m=zu_m+masse(ijq,l)
229     i=mod(i,iim)+1
230     ijq=(j-1)*iip1+i
231     ENDDO
232     ! ajout de la maille non completement advectee
233     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)- &
234     0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
235     ENDIF
236     ENDDO
237     ENDIF
238     ENDDO
239     ENDIF ! n0.gt.0
240    
241    
242    
243     ! bouclage en latitude
244    
245     DO l=1,llm
246     DO ij=iip1+iip1,ip1jm,iip1
247     u_mq(ij,l)=u_mq(ij-iim,l)
248     ENDDO
249     ENDDO
250    
251    
252     ! calcul des tendances
253    
254     DO l=1,llm
255     DO ij=iip2+1,ip1jm
256     new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
257     q(ij,l)=(q(ij,l)*masse(ij,l)+ &
258     u_mq(ij-1,l)-u_mq(ij,l)) &
259     /new_m
260     masse(ij,l)=new_m
261     ENDDO
262     ! Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
263     DO ij=iip1+iip1,ip1jm,iip1
264     q(ij-iim,l)=q(ij,l)
265     masse(ij-iim,l)=masse(ij,l)
266     ENDDO
267     ENDDO
268    
269     ! CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
270     ! CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
271    
272    
273     RETURN
274     END
  ViewVC Help
Powered by ViewVC 1.1.21