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 |