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

Contents of /trunk/dyn3d/vlxqs.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 157 - (show annotations)
Mon Jul 20 16:01:49 2015 UTC (8 years, 10 months ago) by guez
Original Path: trunk/Sources/dyn3d/vlxqs.f
File size: 7860 byte(s)
Just encapsulated SUBROUTINE vlsplt in a module and cleaned it.

In procedure vlx, local variables dxqu and adxqu only need indices
iip2:ip1jm. Otherwise, just cleaned vlx.

Procedures dynredem0 and dynredem1 no longer have argument fichnom,
they just operate on a file named "restart.nc". The programming
guideline here is that gcm should not be more complex than it needs by
itself, other programs (ce0l etc.) just have to adapt to gcm. So ce0l
now creates files "restart.nc" and "restartphy.nc".

In order to facilitate decentralizing the writing of "restartphy.nc",
created a procedure phyredem0 out of phyredem. phyredem0 creates the
NetCDF header of "restartphy.nc" while phyredem writes the NetCDF
variables. As the global attribute itau_phy needs to be filled in
phyredem0, at the beginnig of the run, we must compute its value
instead of just using itap. So we have a dummy argument lmt_pas of
phyredem0. Also, the ncid of "startphy.nc" is upgraded from local
variable of phyetat0 to dummy argument. phyetat0 no longer closes
"startphy.nc".

Following the same decentralizing objective, the ncid of "restart.nc"
is upgraded from local variable of dynredem0 to module variable of
dynredem0_m. "restart.nc" is not closed at the end of dynredem0 nor
opened at the beginning of dynredem1.

In procedure etat0, instead of creating many vectors of size klon
which will be filled with zeroes, just create one array null_array.

In procedure phytrac, instead of writing trs(: 1) to a text file,
write it to "restartphy.nc" (following LMDZ). This is better because
now trs(: 1) is next to its coordinates. We can write to
"restartphy.nc" from phytrac directly, and not add trs(: 1) to the
long list of variables in physiq, thanks to the decentralizing of
"restartphy.nc".

In procedure phyetat0, we no longer write to standard output the
minimum and maximum values of read arrays. It is ok to check input and
abort on invalid values but just printing statistics on input seems too
much useless computation and out of place clutter.

1 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 use disvert_m
14 use conf_gcm_m
15 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 Logical first
39 SAVE first
40
41 REAL temps1,temps2,temps3,temps4,temps5
42 SAVE temps1,temps2,temps3,temps4,temps5
43
44
45 DATA first/.true./
46
47 IF(first) THEN
48 temps1=0.
49 temps2=0.
50 temps3=0.
51 temps4=0.
52 temps5=0.
53 first=.false.
54 ENDIF
55
56 ! calcul de la pente a droite et a gauche de la maille
57
58
59 IF (pente_max.gt.-1.e-5) THEN
60 ! IF (pente_max.gt.10) THEN
61
62 ! calcul des pentes avec limitation, Van Leer scheme I:
63 ! -----------------------------------------------------
64
65 ! calcul de la pente aux points u
66 DO l = 1, llm
67 DO ij=iip2,ip1jm-1
68 dxqu(ij)=q(ij+1,l)-q(ij,l)
69 ! IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
70 ! sigu(ij)=u_m(ij,l)/masse(ij,l)
71 ENDDO
72 DO ij=iip1+iip1,ip1jm,iip1
73 dxqu(ij)=dxqu(ij-iim)
74 ! sigu(ij)=sigu(ij-iim)
75 ENDDO
76
77 DO ij=iip2,ip1jm
78 adxqu(ij)=abs(dxqu(ij))
79 ENDDO
80
81 ! calcul de la pente maximum dans la maille en valeur absolue
82
83 DO ij=iip2+1,ip1jm
84 dxqmax(ij,l)=pente_max* &
85 min(adxqu(ij-1),adxqu(ij))
86 ! limitation subtile
87 ! , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
88
89
90 ENDDO
91
92 DO ij=iip1+iip1,ip1jm,iip1
93 dxqmax(ij-iim,l)=dxqmax(ij,l)
94 ENDDO
95
96 DO ij=iip2+1,ip1jm
97 IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
98 dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
99 ELSE
100 ! extremum local
101 dxq(ij,l)=0.
102 ENDIF
103 dxq(ij,l)=0.5*dxq(ij,l)
104 dxq(ij,l)= &
105 sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
106 ENDDO
107
108 ENDDO ! l=1,llm
109
110 ELSE ! (pente_max.lt.-1.e-5)
111
112 ! Pentes produits:
113 ! ----------------
114
115 DO l = 1, llm
116 DO ij=iip2,ip1jm-1
117 dxqu(ij)=q(ij+1,l)-q(ij,l)
118 ENDDO
119 DO ij=iip1+iip1,ip1jm,iip1
120 dxqu(ij)=dxqu(ij-iim)
121 ENDDO
122
123 DO ij=iip2+1,ip1jm
124 zz(ij)=dxqu(ij-1)*dxqu(ij)
125 zz(ij)=zz(ij)+zz(ij)
126 IF(zz(ij).gt.0) THEN
127 dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
128 ELSE
129 ! extremum local
130 dxq(ij,l)=0.
131 ENDIF
132 ENDDO
133
134 ENDDO
135
136 ENDIF ! (pente_max.lt.-1.e-5)
137
138 ! bouclage de la pente en iip1:
139 ! -----------------------------
140
141 DO l=1,llm
142 DO ij=iip1+iip1,ip1jm,iip1
143 dxq(ij-iim,l)=dxq(ij,l)
144 ENDDO
145
146 DO ij=1,ip1jmp1
147 iadvplus(ij,l)=0
148 ENDDO
149
150 ENDDO
151
152
153 ! calcul des flux a gauche et a droite
154
155 ! on cumule le flux correspondant a toutes les mailles dont la masse
156 ! au travers de la paroi pENDant le pas de temps.
157 ! le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
158 DO l=1,llm
159 DO ij=iip2,ip1jm-1
160 IF (u_m(ij,l).gt.0.) THEN
161 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
162 u_mq(ij,l)=u_m(ij,l)* &
163 min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
164 ELSE
165 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
166 u_mq(ij,l)=u_m(ij,l)* &
167 min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
168 ENDIF
169 ENDDO
170 ENDDO
171
172
173 ! detection des points ou on advecte plus que la masse de la
174 ! maille
175 DO l=1,llm
176 DO ij=iip2,ip1jm-1
177 IF(zdum(ij,l).lt.0) THEN
178 iadvplus(ij,l)=1
179 u_mq(ij,l)=0.
180 ENDIF
181 ENDDO
182 ENDDO
183 DO l=1,llm
184 DO ij=iip1+iip1,ip1jm,iip1
185 iadvplus(ij,l)=iadvplus(ij-iim,l)
186 ENDDO
187 ENDDO
188
189
190
191 ! traitement special pour le cas ou on advecte en longitude plus que le
192 ! contenu de la maille.
193 ! cette partie est mal vectorisee.
194
195 ! pas d'influence de la pression saturante (pour l'instant)
196
197 ! calcul du nombre de maille sur lequel on advecte plus que la maille.
198
199 n0=0
200 DO l=1,llm
201 nl(l)=0
202 DO ij=iip2,ip1jm
203 nl(l)=nl(l)+iadvplus(ij,l)
204 ENDDO
205 n0=n0+nl(l)
206 ENDDO
207
208 IF(n0.gt.0) THEN
209 DO l=1,llm
210 IF(nl(l).gt.0) THEN
211 iju=0
212 ! indicage des mailles concernees par le traitement special
213 DO ij=iip2,ip1jm
214 IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
215 iju=iju+1
216 indu(iju)=ij
217 ENDIF
218 ENDDO
219 niju=iju
220
221 ! traitement des mailles
222 DO iju=1,niju
223 ij=indu(iju)
224 j=(ij-1)/iip1+1
225 zu_m=u_m(ij,l)
226 u_mq(ij,l)=0.
227 IF(zu_m.gt.0.) THEN
228 ijq=ij
229 i=ijq-(j-1)*iip1
230 ! accumulation pour les mailles completements advectees
231 do while(zu_m.gt.masse(ijq,l))
232 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
233 zu_m=zu_m-masse(ijq,l)
234 i=mod(i-2+iim,iim)+1
235 ijq=(j-1)*iip1+i
236 ENDDO
237 ! ajout de la maille non completement advectee
238 u_mq(ij,l)=u_mq(ij,l)+zu_m* &
239 (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
240 ELSE
241 ijq=ij+1
242 i=ijq-(j-1)*iip1
243 ! accumulation pour les mailles completements advectees
244 do while(-zu_m.gt.masse(ijq,l))
245 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
246 zu_m=zu_m+masse(ijq,l)
247 i=mod(i,iim)+1
248 ijq=(j-1)*iip1+i
249 ENDDO
250 ! ajout de la maille non completement advectee
251 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)- &
252 0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
253 ENDIF
254 ENDDO
255 ENDIF
256 ENDDO
257 ENDIF ! n0.gt.0
258
259
260
261 ! bouclage en latitude
262
263 DO l=1,llm
264 DO ij=iip1+iip1,ip1jm,iip1
265 u_mq(ij,l)=u_mq(ij-iim,l)
266 ENDDO
267 ENDDO
268
269
270 ! calcul des tendances
271
272 DO l=1,llm
273 DO ij=iip2+1,ip1jm
274 new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
275 q(ij,l)=(q(ij,l)*masse(ij,l)+ &
276 u_mq(ij-1,l)-u_mq(ij,l)) &
277 /new_m
278 masse(ij,l)=new_m
279 ENDDO
280 ! Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
281 DO ij=iip1+iip1,ip1jm,iip1
282 q(ij-iim,l)=q(ij,l)
283 masse(ij-iim,l)=masse(ij,l)
284 ENDDO
285 ENDDO
286
287 ! CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
288 ! CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
289
290
291 RETURN
292 END

  ViewVC Help
Powered by ViewVC 1.1.21