1 |
SUBROUTINE advnx(q, qg, qd, masse, u_m, mode) |
2 |
|
3 |
! Auteur : F. Hourdin |
4 |
|
5 |
! ******************************************************************** |
6 |
! Shema d'advection " pseudo amont " . |
7 |
! ******************************************************************** |
8 |
! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... |
9 |
|
10 |
|
11 |
! -------------------------------------------------------------------- |
12 |
USE dimensions |
13 |
USE paramet_m |
14 |
USE comconst |
15 |
USE disvert_m |
16 |
USE conf_gcm_m |
17 |
IMPLICIT NONE |
18 |
|
19 |
|
20 |
|
21 |
! Arguments: |
22 |
! ---------- |
23 |
INTEGER mode |
24 |
REAL masse(ip1jmp1, llm) |
25 |
REAL u_m(ip1jmp1, llm) |
26 |
REAL q(ip1jmp1, llm), qd(ip1jmp1, llm), qg(ip1jmp1, llm) |
27 |
|
28 |
! Local |
29 |
! --------- |
30 |
|
31 |
INTEGER i, j, ij, l, indu(ip1jmp1), niju, iju, ijq |
32 |
INTEGER n0, nl(llm) |
33 |
|
34 |
REAL new_m, zu_m, zdq, zz |
35 |
REAL zsigg(ip1jmp1, llm), zsigd(ip1jmp1, llm), zsig |
36 |
REAL u_mq(ip1jmp1, llm) |
37 |
|
38 |
REAL zm, zq, zsigm, zsigp, zqm, zqp, zu |
39 |
|
40 |
LOGICAL ladvplus(ip1jmp1, llm) |
41 |
|
42 |
REAL prec |
43 |
SAVE prec |
44 |
|
45 |
DATA prec/1.E-15/ |
46 |
|
47 |
DO l = 1, llm |
48 |
DO ij = iip2, ip1jm |
49 |
zdq = qd(ij, l) - qg(ij, l) |
50 |
IF (abs(zdq)>prec) THEN |
51 |
zsigd(ij, l) = (q(ij,l)-qg(ij,l))/zdq |
52 |
zsigg(ij, l) = 1. - zsigd(ij, l) |
53 |
ELSE |
54 |
zsigd(ij, l) = 0.5 |
55 |
zsigg(ij, l) = 0.5 |
56 |
qd(ij, l) = q(ij, l) |
57 |
qg(ij, l) = q(ij, l) |
58 |
END IF |
59 |
END DO |
60 |
END DO |
61 |
|
62 |
! calcul de la pente maximum dans la maille en valeur absolue |
63 |
|
64 |
DO l = 1, llm |
65 |
DO ij = iip2, ip1jm - 1 |
66 |
IF (u_m(ij,l)>=0.) THEN |
67 |
zsigp = zsigd(ij, l) |
68 |
zsigm = zsigg(ij, l) |
69 |
zqp = qd(ij, l) |
70 |
zqm = qg(ij, l) |
71 |
zm = masse(ij, l) |
72 |
zq = q(ij, l) |
73 |
ELSE |
74 |
zsigm = zsigd(ij+1, l) |
75 |
zsigp = zsigg(ij+1, l) |
76 |
zqm = qd(ij+1, l) |
77 |
zqp = qg(ij+1, l) |
78 |
zm = masse(ij+1, l) |
79 |
zq = q(ij+1, l) |
80 |
END IF |
81 |
zu = abs(u_m(ij,l)) |
82 |
ladvplus(ij, l) = zu > zm |
83 |
zsig = zu/zm |
84 |
IF (zsig==0.) zsigp = 0.1 |
85 |
IF (mode==1) THEN |
86 |
IF (zsig<=zsigp) THEN |
87 |
u_mq(ij, l) = u_m(ij, l)*zqp |
88 |
ELSE IF (mode==1) THEN |
89 |
u_mq(ij, l) = sign(zm, u_m(ij,l))*(zsigp*zqp+(zsig-zsigp)*zqm) |
90 |
END IF |
91 |
ELSE |
92 |
IF (zsig<=zsigp) THEN |
93 |
u_mq(ij, l) = u_m(ij, l)*(zqp-0.5*zsig/zsigp*(zqp-zq)) |
94 |
ELSE |
95 |
zz = 0.5*(zsig-zsigp)/zsigm |
96 |
u_mq(ij, l) = sign(zm, u_m(ij,l))*(0.5*(zq+zqp)*zsigp+(zsig-zsigp)* & |
97 |
(zq+zz*(zqm-zq))) |
98 |
END IF |
99 |
END IF |
100 |
END DO |
101 |
END DO |
102 |
|
103 |
DO l = 1, llm |
104 |
DO ij = iip1 + iip1, ip1jm, iip1 |
105 |
u_mq(ij, l) = u_mq(ij-iim, l) |
106 |
ladvplus(ij, l) = ladvplus(ij-iim, l) |
107 |
END DO |
108 |
END DO |
109 |
|
110 |
! ================================================================= |
111 |
! SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES |
112 |
! ================================================================= |
113 |
! tris des regions a traiter |
114 |
n0 = 0 |
115 |
DO l = 1, llm |
116 |
nl(l) = 0 |
117 |
DO ij = iip2, ip1jm |
118 |
IF (ladvplus(ij,l)) THEN |
119 |
nl(l) = nl(l) + 1 |
120 |
u_mq(ij, l) = 0. |
121 |
END IF |
122 |
END DO |
123 |
n0 = n0 + nl(l) |
124 |
END DO |
125 |
|
126 |
IF (n0>1) THEN |
127 |
DO l = 1, llm |
128 |
IF (nl(l)>0) THEN |
129 |
iju = 0 |
130 |
! indicage des mailles concernees par le traitement special |
131 |
DO ij = iip2, ip1jm |
132 |
IF (ladvplus(ij,l) .AND. mod(ij,iip1)/=0) THEN |
133 |
iju = iju + 1 |
134 |
indu(iju) = ij |
135 |
END IF |
136 |
END DO |
137 |
niju = iju |
138 |
|
139 |
! traitement des mailles |
140 |
DO iju = 1, niju |
141 |
ij = indu(iju) |
142 |
j = (ij-1)/iip1 + 1 |
143 |
zu_m = u_m(ij, l) |
144 |
u_mq(ij, l) = 0. |
145 |
IF (zu_m>0.) THEN |
146 |
ijq = ij |
147 |
i = ijq - (j-1)*iip1 |
148 |
! accumulation pour les mailles completements advectees |
149 |
DO WHILE (zu_m>masse(ijq,l)) |
150 |
u_mq(ij, l) = u_mq(ij, l) + q(ijq, l)*masse(ijq, l) |
151 |
zu_m = zu_m - masse(ijq, l) |
152 |
i = mod(i-2+iim, iim) + 1 |
153 |
ijq = (j-1)*iip1 + i |
154 |
END DO |
155 |
! MODIFS SPECIFIQUES DU SCHEMA |
156 |
! ajout de la maille non completement advectee |
157 |
zsig = zu_m/masse(ijq, l) |
158 |
IF (zsig<=zsigd(ijq,l)) THEN |
159 |
u_mq(ij, l) = u_mq(ij, l) + zu_m*(qd(ijq,l)-0.5*zsig/zsigd(ijq, & |
160 |
l)*(qd(ijq,l)-q(ijq,l))) |
161 |
ELSE |
162 |
! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l) |
163 |
! goto 8888 |
164 |
zz = 0.5*(zsig-zsigd(ijq,l))/zsigg(ijq, l) |
165 |
IF (.NOT. (zz>0. .AND. zz<=0.5)) THEN |
166 |
PRINT *, 'probleme2 au point ij=', ij, ' l=', l |
167 |
PRINT *, 'zz=', zz |
168 |
STOP |
169 |
END IF |
170 |
u_mq(ij, l) = u_mq(ij, l) + masse(ijq, l)*(0.5*(q(ijq, & |
171 |
l)+qd(ijq,l))*zsigd(ijq,l)+(zsig-zsigd(ijq,l))*(q(ijq, & |
172 |
l)+zz*(qg(ijq,l)-q(ijq,l)))) |
173 |
END IF |
174 |
ELSE |
175 |
ijq = ij + 1 |
176 |
i = ijq - (j-1)*iip1 |
177 |
! accumulation pour les mailles completements advectees |
178 |
DO WHILE (-zu_m>masse(ijq,l)) |
179 |
u_mq(ij, l) = u_mq(ij, l) - q(ijq, l)*masse(ijq, l) |
180 |
zu_m = zu_m + masse(ijq, l) |
181 |
i = mod(i, iim) + 1 |
182 |
ijq = (j-1)*iip1 + i |
183 |
END DO |
184 |
! ajout de la maille non completement advectee |
185 |
! 2eme MODIF SPECIFIQUE |
186 |
zsig = -zu_m/masse(ij+1, l) |
187 |
IF (zsig<=zsigg(ijq,l)) THEN |
188 |
u_mq(ij, l) = u_mq(ij, l) + zu_m*(qg(ijq,l)-0.5*zsig/zsigg(ijq, & |
189 |
l)*(qg(ijq,l)-q(ijq,l))) |
190 |
ELSE |
191 |
! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l) |
192 |
! goto 9999 |
193 |
zz = 0.5*(zsig-zsigg(ijq,l))/zsigd(ijq, l) |
194 |
IF (.NOT. (zz>0. .AND. zz<=0.5)) THEN |
195 |
PRINT *, 'probleme22 au point ij=', ij, ' l=', l |
196 |
PRINT *, 'zz=', zz |
197 |
STOP |
198 |
END IF |
199 |
u_mq(ij, l) = u_mq(ij, l) - masse(ijq, l)*(0.5*(q(ijq, & |
200 |
l)+qg(ijq,l))*zsigg(ijq,l)+(zsig-zsigg(ijq,l))*(q(ijq, & |
201 |
l)+zz*(qd(ijq,l)-q(ijq,l)))) |
202 |
END IF |
203 |
! fin de la modif |
204 |
END IF |
205 |
END DO |
206 |
END IF |
207 |
END DO |
208 |
END IF ! n0.gt.0 |
209 |
|
210 |
! bouclage en latitude |
211 |
DO l = 1, llm |
212 |
DO ij = iip1 + iip1, ip1jm, iip1 |
213 |
u_mq(ij, l) = u_mq(ij-iim, l) |
214 |
END DO |
215 |
END DO |
216 |
|
217 |
! ================================================================= |
218 |
! CALCUL DE LA CONVERGENCE DES FLUX |
219 |
! ================================================================= |
220 |
|
221 |
DO l = 1, llm |
222 |
DO ij = iip2 + 1, ip1jm |
223 |
new_m = masse(ij, l) + u_m(ij-1, l) - u_m(ij, l) |
224 |
q(ij, l) = (q(ij,l)*masse(ij,l)+u_mq(ij-1,l)-u_mq(ij,l))/new_m |
225 |
masse(ij, l) = new_m |
226 |
END DO |
227 |
! Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous) |
228 |
DO ij = iip1 + iip1, ip1jm, iip1 |
229 |
q(ij-iim, l) = q(ij, l) |
230 |
masse(ij-iim, l) = masse(ij, l) |
231 |
END DO |
232 |
END DO |
233 |
|
234 |
RETURN |
235 |
END SUBROUTINE advnx |