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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 157 - (hide annotations)
Mon Jul 20 16:01:49 2015 UTC (8 years, 11 months ago) by guez
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 guez 157 module vlx_m
2 guez 31
3 guez 157 IMPLICIT NONE
4 guez 31
5 guez 157 contains
6 guez 31
7 guez 157 SUBROUTINE vlx(q, pente_max, masse, u_m)
8 guez 31
9 guez 157 ! Authors: P. Le Van, F. Hourdin, F. Forget
10 guez 31
11 guez 157 ! Sch\'ema d'advection "pseudo-amont".
12 guez 31
13 guez 157 use dimens_m, only: iim, llm
14     use paramet_m, only: ip1jmp1, iip1, iip2, ip1jm
15 guez 31
16 guez 157 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 guez 31
21 guez 157 ! Local:
22 guez 31
23 guez 157 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 guez 31
31 guez 157 !-----------------------------------------------------
32 guez 31
33 guez 157 ! calcul de la pente a droite et a gauche de la maille
34 guez 31
35 guez 157 IF (pente_max > - 1e-5) THEN
36     ! calcul des pentes avec limitation, Van Leer scheme I:
37 guez 31
38 guez 157 ! 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 guez 31
47 guez 157 DO ij = iip2, ip1jm
48     adxqu(ij) = abs(dxqu(ij))
49     ENDDO
50 guez 31
51 guez 157 ! calcul de la pente maximum dans la maille en valeur absolue
52 guez 31
53 guez 157 DO ij = iip2 + 1, ip1jm
54     dxqmax(ij, l) = pente_max * min(adxqu(ij - 1), adxqu(ij))
55     ENDDO
56 guez 31
57 guez 157 DO ij = iip1 + iip1, ip1jm, iip1
58     dxqmax(ij - iim, l) = dxqmax(ij, l)
59     ENDDO
60 guez 31
61 guez 157 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 guez 31
75 guez 157 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 guez 31
83 guez 157 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 guez 31
96 guez 157 ! bouclage de la pente en iip1:
97 guez 31
98 guez 157 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 guez 31
107 guez 157 ! calcul des flux a gauche et a droite
108 guez 31
109 guez 157 ! on cumule le flux correspondant a toutes les mailles dont la masse
110     ! au travers de la paroi pENDant le pas de temps.
111 guez 31
112 guez 157 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 guez 31 ELSE
119 guez 157 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 guez 31 ENDIF
123     ENDDO
124 guez 157 ENDDO
125 guez 31
126 guez 157 ! 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 guez 31 ENDDO
135 guez 157 ENDDO
136 guez 31
137 guez 157 DO l = 1, llm
138     DO ij = iip1 + iip1, ip1jm, iip1
139     iadvplus(ij, l) = iadvplus(ij - iim, l)
140     ENDDO
141     ENDDO
142 guez 31
143 guez 157 ! traitement special pour le cas ou on advecte en longitude plus
144     ! que le contenu de la maille. cette partie est mal vectorisee.
145 guez 31
146 guez 157 ! calcul du nombre de maille sur lequel on advecte plus que la maille.
147 guez 31
148 guez 157 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 guez 31
157 guez 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 guez 31
170 guez 157 ! 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 guez 31
208 guez 157 ! 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 guez 31
215 guez 157 ! calcul des tENDances
216 guez 31
217 guez 157 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 guez 31
225 guez 157 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 guez 31
231 guez 157 END SUBROUTINE vlx
232 guez 31
233 guez 157 end module vlx_m

  ViewVC Help
Powered by ViewVC 1.1.21