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

Contents of /trunk/dyn3d/Vlsplt/vlx.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/Vlsplt/vlx.f
File size: 6925 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 module vlx_m
2
3 IMPLICIT NONE
4
5 contains
6
7 SUBROUTINE vlx(q, pente_max, masse, u_m)
8
9 ! Authors: P. Le Van, F. Hourdin, F. Forget
10
11 ! Sch\'ema d'advection "pseudo-amont".
12
13 use dimens_m, only: iim, llm
14 use paramet_m, only: ip1jmp1, iip1, iip2, ip1jm
15
16 REAL, intent(inout):: q(ip1jmp1, llm)
17 REAL, intent(in):: pente_max
18 real, intent(inout):: masse(ip1jmp1, llm)
19 REAL, intent(in):: u_m(ip1jmp1, llm)
20
21 ! Local:
22
23 INTEGER ij, l, j, i, iju, ijq, indu(ip1jmp1), niju
24 INTEGER n0, iadvplus(ip1jmp1, llm), nl(llm)
25 REAL new_m, zu_m, zdum(ip1jmp1, llm)
26 REAL dxq(ip1jmp1, llm), dxqu(iip2:ip1jm)
27 REAL zz(ip1jmp1)
28 REAL adxqu(iip2:ip1jm), dxqmax(ip1jmp1, llm)
29 REAL u_mq(ip1jmp1, llm)
30
31 !-----------------------------------------------------
32
33 ! calcul de la pente a droite et a gauche de la maille
34
35 IF (pente_max > - 1e-5) THEN
36 ! calcul des pentes avec limitation, Van Leer scheme I:
37
38 ! calcul de la pente aux points u
39 DO l = 1, llm
40 DO ij = iip2, ip1jm - 1
41 dxqu(ij) = q(ij + 1, l) - q(ij, l)
42 ENDDO
43 DO ij = iip1 + iip1, ip1jm, iip1
44 dxqu(ij) = dxqu(ij - iim)
45 ENDDO
46
47 DO ij = iip2, ip1jm
48 adxqu(ij) = abs(dxqu(ij))
49 ENDDO
50
51 ! calcul de la pente maximum dans la maille en valeur absolue
52
53 DO ij = iip2 + 1, ip1jm
54 dxqmax(ij, l) = pente_max * min(adxqu(ij - 1), adxqu(ij))
55 ENDDO
56
57 DO ij = iip1 + iip1, ip1jm, iip1
58 dxqmax(ij - iim, l) = dxqmax(ij, l)
59 ENDDO
60
61 DO ij = iip2 + 1, ip1jm
62 IF (dxqu(ij - 1) * dxqu(ij) > 0.) THEN
63 dxq(ij, l) = dxqu(ij - 1) + dxqu(ij)
64 ELSE
65 ! extremum local
66 dxq(ij, l) = 0.
67 ENDIF
68 dxq(ij, l) = 0.5 * dxq(ij, l)
69 dxq(ij, l) = sign(min(abs(dxq(ij, l)), dxqmax(ij, l)), dxq(ij, l))
70 ENDDO
71 ENDDO
72 ELSE
73 ! Pentes produits:
74
75 DO l = 1, llm
76 DO ij = iip2, ip1jm - 1
77 dxqu(ij) = q(ij + 1, l) - q(ij, l)
78 ENDDO
79 DO ij = iip1 + iip1, ip1jm, iip1
80 dxqu(ij) = dxqu(ij - iim)
81 ENDDO
82
83 DO ij = iip2 + 1, ip1jm
84 zz(ij) = dxqu(ij - 1) * dxqu(ij)
85 zz(ij) = zz(ij) + zz(ij)
86 IF (zz(ij) > 0) THEN
87 dxq(ij, l) = zz(ij) / (dxqu(ij - 1) + dxqu(ij))
88 ELSE
89 ! extremum local
90 dxq(ij, l) = 0.
91 ENDIF
92 ENDDO
93 ENDDO
94 ENDIF
95
96 ! bouclage de la pente en iip1:
97
98 DO l = 1, llm
99 DO ij = iip1 + iip1, ip1jm, iip1
100 dxq(ij - iim, l) = dxq(ij, l)
101 ENDDO
102 DO ij = 1, ip1jmp1
103 iadvplus(ij, l) = 0
104 ENDDO
105 ENDDO
106
107 ! calcul des flux a gauche et a droite
108
109 ! on cumule le flux correspondant a toutes les mailles dont la masse
110 ! au travers de la paroi pENDant le pas de temps.
111
112 DO l = 1, llm
113 DO ij = iip2, ip1jm - 1
114 IF (u_m(ij, l) > 0.) THEN
115 zdum(ij, l) = 1. - u_m(ij, l) / masse(ij, l)
116 u_mq(ij, l) = u_m(ij, l) &
117 * (q(ij, l) + 0.5 * zdum(ij, l) * dxq(ij, l))
118 ELSE
119 zdum(ij, l) = 1. + u_m(ij, l) / masse(ij + 1, l)
120 u_mq(ij, l) = u_m(ij, l) &
121 * (q(ij + 1, l) - 0.5 * zdum(ij, l) * dxq(ij + 1, l))
122 ENDIF
123 ENDDO
124 ENDDO
125
126 ! detection des points ou on advecte plus que la masse de la
127 ! maille
128 DO l = 1, llm
129 DO ij = iip2, ip1jm - 1
130 IF (zdum(ij, l) < 0) THEN
131 iadvplus(ij, l) = 1
132 u_mq(ij, l) = 0.
133 ENDIF
134 ENDDO
135 ENDDO
136
137 DO l = 1, llm
138 DO ij = iip1 + iip1, ip1jm, iip1
139 iadvplus(ij, l) = iadvplus(ij - iim, l)
140 ENDDO
141 ENDDO
142
143 ! traitement special pour le cas ou on advecte en longitude plus
144 ! que le contenu de la maille. cette partie est mal vectorisee.
145
146 ! calcul du nombre de maille sur lequel on advecte plus que la maille.
147
148 n0 = 0
149 DO l = 1, llm
150 nl(l) = 0
151 DO ij = iip2, ip1jm
152 nl(l) = nl(l) + iadvplus(ij, l)
153 ENDDO
154 n0 = n0 + nl(l)
155 ENDDO
156
157 IF (n0 > 0) THEN
158 DO l = 1, llm
159 IF (nl(l) > 0) THEN
160 iju = 0
161 ! indicage des mailles concernees par le traitement special
162 DO ij = iip2, ip1jm
163 IF (iadvplus(ij, l) == 1.and.mod(ij, iip1) /= 0) THEN
164 iju = iju + 1
165 indu(iju) = ij
166 ENDIF
167 ENDDO
168 niju = iju
169
170 ! traitement des mailles
171 DO iju = 1, niju
172 ij = indu(iju)
173 j = (ij - 1) / iip1 + 1
174 zu_m = u_m(ij, l)
175 u_mq(ij, l) = 0.
176 IF (zu_m > 0.) THEN
177 ijq = ij
178 i = ijq - (j - 1) * iip1
179 ! accumulation pour les mailles completements advectees
180 do while(zu_m > masse(ijq, l))
181 u_mq(ij, l) = u_mq(ij, l) + q(ijq, l) * masse(ijq, l)
182 zu_m = zu_m - masse(ijq, l)
183 i = mod(i - 2 + iim, iim) + 1
184 ijq = (j - 1) * iip1 + i
185 ENDDO
186 ! ajout de la maille non completement advectee
187 u_mq(ij, l) = u_mq(ij, l) + zu_m * (q(ijq, l) + 0.5 * (1. &
188 - zu_m / masse(ijq, l)) * dxq(ijq, l))
189 ELSE
190 ijq = ij + 1
191 i = ijq - (j - 1) * iip1
192 ! accumulation pour les mailles completements advectees
193 do while(- zu_m > masse(ijq, l))
194 u_mq(ij, l) = u_mq(ij, l) - q(ijq, l) * masse(ijq, l)
195 zu_m = zu_m + masse(ijq, l)
196 i = mod(i, iim) + 1
197 ijq = (j - 1) * iip1 + i
198 ENDDO
199 ! ajout de la maille non completement advectee
200 u_mq(ij, l) = u_mq(ij, l) + zu_m * (q(ijq, l) - 0.5 * (1. &
201 + zu_m / masse(ijq, l)) * dxq(ijq, l))
202 ENDIF
203 ENDDO
204 ENDIF
205 ENDDO
206 ENDIF
207
208 ! bouclage en latitude
209 DO l = 1, llm
210 DO ij = iip1 + iip1, ip1jm, iip1
211 u_mq(ij, l) = u_mq(ij - iim, l)
212 ENDDO
213 ENDDO
214
215 ! calcul des tENDances
216
217 DO l = 1, llm
218 DO ij = iip2 + 1, ip1jm
219 new_m = masse(ij, l) + u_m(ij - 1, l) - u_m(ij, l)
220 q(ij, l) = (q(ij, l) * masse(ij, l) + u_mq(ij - 1, l) &
221 - u_mq(ij, l)) / new_m
222 masse(ij, l) = new_m
223 ENDDO
224
225 DO ij = iip1 + iip1, ip1jm, iip1
226 q(ij - iim, l) = q(ij, l)
227 masse(ij - iim, l) = masse(ij, l)
228 ENDDO
229 ENDDO
230
231 END SUBROUTINE vlx
232
233 end module vlx_m

  ViewVC Help
Powered by ViewVC 1.1.21