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

Annotation of /trunk/dyn3d/vlxqs.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 157 - (hide 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 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 guez 157 Logical first
39     SAVE first
40 guez 43
41 guez 105 REAL temps1,temps2,temps3,temps4,temps5
42     SAVE temps1,temps2,temps3,temps4,temps5
43 guez 43
44    
45 guez 157 DATA first/.true./
46 guez 43
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