1 | ! |
---|
2 | ! $Header$ |
---|
3 | ! |
---|
4 | SUBROUTINE pentes_ini (q,w,masse,pbaru,pbarv,mode) |
---|
5 | |
---|
6 | USE comconst_mod, ONLY: pi, dtvr |
---|
7 | |
---|
8 | IMPLICIT NONE |
---|
9 | |
---|
10 | c======================================================================= |
---|
11 | c Adaptation LMDZ: A.Armengaud (LGGE) |
---|
12 | c ---------------- |
---|
13 | c |
---|
14 | c ******************************************************************** |
---|
15 | c Transport des traceurs par la methode des pentes |
---|
16 | c ******************************************************************** |
---|
17 | c Reference possible : Russel. G.L., Lerner J.A.: |
---|
18 | c A new Finite-Differencing Scheme for Traceur Transport |
---|
19 | c Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81 |
---|
20 | c ******************************************************************** |
---|
21 | c q,w,masse,pbaru et pbarv |
---|
22 | c sont des arguments d'entree pour le s-pg .... |
---|
23 | c |
---|
24 | c======================================================================= |
---|
25 | |
---|
26 | |
---|
27 | include "dimensions.h" |
---|
28 | include "paramet.h" |
---|
29 | include "comgeom2.h" |
---|
30 | |
---|
31 | c Arguments: |
---|
32 | c ---------- |
---|
33 | integer mode |
---|
34 | REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm ) |
---|
35 | REAL q( iip1,jjp1,llm,0:3) |
---|
36 | REAL w( ip1jmp1,llm ) |
---|
37 | REAL masse( iip1,jjp1,llm) |
---|
38 | c Local: |
---|
39 | c ------ |
---|
40 | LOGICAL limit |
---|
41 | REAL sm ( iip1,jjp1, llm ) |
---|
42 | REAL s0( iip1,jjp1,llm ), sx( iip1,jjp1,llm ) |
---|
43 | REAL sy( iip1,jjp1,llm ), sz( iip1,jjp1,llm ) |
---|
44 | real masn,mass,zz |
---|
45 | INTEGER i,j,l,iq |
---|
46 | |
---|
47 | c modif Fred 24 03 96 |
---|
48 | |
---|
49 | real sinlon(iip1),sinlondlon(iip1) |
---|
50 | real coslon(iip1),coslondlon(iip1) |
---|
51 | save sinlon,coslon,sinlondlon,coslondlon |
---|
52 | real dyn1,dyn2,dys1,dys2 |
---|
53 | real qpn,qps,dqzpn,dqzps |
---|
54 | real smn,sms,s0n,s0s,sxn(iip1),sxs(iip1) |
---|
55 | real qmin,zq,pente_max |
---|
56 | c |
---|
57 | REAL SSUM |
---|
58 | integer ismax,ismin,lati,latf |
---|
59 | EXTERNAL SSUM, ismin,ismax |
---|
60 | logical first |
---|
61 | save first |
---|
62 | c fin modif |
---|
63 | |
---|
64 | c EXTERNAL masskg |
---|
65 | EXTERNAL advx |
---|
66 | EXTERNAL advy |
---|
67 | EXTERNAL advz |
---|
68 | |
---|
69 | c modif Fred 24 03 96 |
---|
70 | data first/.true./ |
---|
71 | |
---|
72 | limit = .TRUE. |
---|
73 | pente_max=2 |
---|
74 | c if (mode.eq.1.or.mode.eq.3) then |
---|
75 | c if (mode.eq.1) then |
---|
76 | if (mode.ge.1) then |
---|
77 | lati=2 |
---|
78 | latf=jjm |
---|
79 | else |
---|
80 | lati=1 |
---|
81 | latf=jjp1 |
---|
82 | endif |
---|
83 | |
---|
84 | qmin=0.4995 |
---|
85 | qmin=0. |
---|
86 | if(first) then |
---|
87 | print*,'SCHEMA AMONT NOUVEAU' |
---|
88 | first=.false. |
---|
89 | do i=2,iip1 |
---|
90 | coslon(i)=cos(rlonv(i)) |
---|
91 | sinlon(i)=sin(rlonv(i)) |
---|
92 | coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi |
---|
93 | sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi |
---|
94 | print*,coslondlon(i),sinlondlon(i) |
---|
95 | enddo |
---|
96 | coslon(1)=coslon(iip1) |
---|
97 | coslondlon(1)=coslondlon(iip1) |
---|
98 | sinlon(1)=sinlon(iip1) |
---|
99 | sinlondlon(1)=sinlondlon(iip1) |
---|
100 | print*,'sum sinlondlon ',ssum(iim,sinlondlon,1)/sinlondlon(1) |
---|
101 | print*,'sum coslondlon ',ssum(iim,coslondlon,1)/coslondlon(1) |
---|
102 | DO l = 1,llm |
---|
103 | DO j = 1,jjp1 |
---|
104 | DO i = 1,iip1 |
---|
105 | q ( i,j,l,1 )=0. |
---|
106 | q ( i,j,l,2 )=0. |
---|
107 | q ( i,j,l,3 )=0. |
---|
108 | ENDDO |
---|
109 | ENDDO |
---|
110 | ENDDO |
---|
111 | |
---|
112 | endif |
---|
113 | c Fin modif Fred |
---|
114 | |
---|
115 | c *** q contient les qqtes de traceur avant l'advection |
---|
116 | |
---|
117 | c *** Affectation des tableaux S a partir de Q |
---|
118 | c *** Rem : utilisation de SCOPY ulterieurement |
---|
119 | |
---|
120 | DO l = 1,llm |
---|
121 | DO j = 1,jjp1 |
---|
122 | DO i = 1,iip1 |
---|
123 | s0( i,j,llm+1-l ) = q ( i,j,l,0 ) |
---|
124 | sx( i,j,llm+1-l ) = q ( i,j,l,1 ) |
---|
125 | sy( i,j,llm+1-l ) = q ( i,j,l,2 ) |
---|
126 | sz( i,j,llm+1-l ) = q ( i,j,l,3 ) |
---|
127 | ENDDO |
---|
128 | ENDDO |
---|
129 | ENDDO |
---|
130 | |
---|
131 | c PRINT*,'----- S0 just before conversion -------' |
---|
132 | c PRINT*,'S0(16,12,1)=',s0(16,12,1) |
---|
133 | c PRINT*,'Q(16,12,1,4)=',q(16,12,1,4) |
---|
134 | |
---|
135 | c *** On calcule la masse d'air en kg |
---|
136 | |
---|
137 | DO l = 1,llm |
---|
138 | DO j = 1,jjp1 |
---|
139 | DO i = 1,iip1 |
---|
140 | sm ( i,j,llm+1-l)=masse( i,j,l ) |
---|
141 | ENDDO |
---|
142 | ENDDO |
---|
143 | ENDDO |
---|
144 | |
---|
145 | c *** On converti les champs S en atome (resp. kg) |
---|
146 | c *** Les routines d'advection traitent les champs |
---|
147 | c *** a advecter si ces derniers sont en atome (resp. kg) |
---|
148 | c *** A optimiser !!! |
---|
149 | |
---|
150 | DO l = 1,llm |
---|
151 | DO j = 1,jjp1 |
---|
152 | DO i = 1,iip1 |
---|
153 | s0(i,j,l) = s0(i,j,l) * sm ( i,j,l ) |
---|
154 | sx(i,j,l) = sx(i,j,l) * sm ( i,j,l ) |
---|
155 | sy(i,j,l) = sy(i,j,l) * sm ( i,j,l ) |
---|
156 | sz(i,j,l) = sz(i,j,l) * sm ( i,j,l ) |
---|
157 | ENDDO |
---|
158 | ENDDO |
---|
159 | ENDDO |
---|
160 | |
---|
161 | c ss0 = 0. |
---|
162 | c DO l = 1,llm |
---|
163 | c DO j = 1,jjp1 |
---|
164 | c DO i = 1,iim |
---|
165 | c ss0 = ss0 + s0 ( i,j,l ) |
---|
166 | c ENDDO |
---|
167 | c ENDDO |
---|
168 | c ENDDO |
---|
169 | c PRINT*, 'valeur tot s0 avant advection=',ss0 |
---|
170 | |
---|
171 | c *** Appel des subroutines d'advection en X, en Y et en Z |
---|
172 | c *** Advection avec "time-splitting" |
---|
173 | |
---|
174 | c----------------------------------------------------------- |
---|
175 | c PRINT*,'----- S0 just before ADVX -------' |
---|
176 | c PRINT*,'S0(16,12,1)=',s0(16,12,1) |
---|
177 | |
---|
178 | c----------------------------------------------------------- |
---|
179 | c do l=1,llm |
---|
180 | c do j=1,jjp1 |
---|
181 | c do i=1,iip1 |
---|
182 | c zq=s0(i,j,l)/sm(i,j,l) |
---|
183 | c if(zq.lt.qmin) |
---|
184 | c , print*,'avant advx1, s0(',i,',',j,',',l,')=',zq |
---|
185 | c enddo |
---|
186 | c enddo |
---|
187 | c enddo |
---|
188 | CCC |
---|
189 | if(mode.eq.2) then |
---|
190 | do l=1,llm |
---|
191 | s0s=0. |
---|
192 | s0n=0. |
---|
193 | dyn1=0. |
---|
194 | dys1=0. |
---|
195 | dyn2=0. |
---|
196 | dys2=0. |
---|
197 | smn=0. |
---|
198 | sms=0. |
---|
199 | do i=1,iim |
---|
200 | smn=smn+sm(i,1,l) |
---|
201 | sms=sms+sm(i,jjp1,l) |
---|
202 | s0n=s0n+s0(i,1,l) |
---|
203 | s0s=s0s+s0(i,jjp1,l) |
---|
204 | zz=sy(i,1,l)/sm(i,1,l) |
---|
205 | dyn1=dyn1+sinlondlon(i)*zz |
---|
206 | dyn2=dyn2+coslondlon(i)*zz |
---|
207 | zz=sy(i,jjp1,l)/sm(i,jjp1,l) |
---|
208 | dys1=dys1+sinlondlon(i)*zz |
---|
209 | dys2=dys2+coslondlon(i)*zz |
---|
210 | enddo |
---|
211 | do i=1,iim |
---|
212 | sy(i,1,l)=dyn1*sinlon(i)+dyn2*coslon(i) |
---|
213 | sy(i,jjp1,l)=dys1*sinlon(i)+dys2*coslon(i) |
---|
214 | enddo |
---|
215 | do i=1,iim |
---|
216 | s0(i,1,l)=s0n/smn+sy(i,1,l) |
---|
217 | s0(i,jjp1,l)=s0s/sms-sy(i,jjp1,l) |
---|
218 | enddo |
---|
219 | |
---|
220 | s0(iip1,1,l)=s0(1,1,l) |
---|
221 | s0(iip1,jjp1,l)=s0(1,jjp1,l) |
---|
222 | |
---|
223 | do i=1,iim |
---|
224 | sxn(i)=s0(i+1,1,l)-s0(i,1,l) |
---|
225 | sxs(i)=s0(i+1,jjp1,l)-s0(i,jjp1,l) |
---|
226 | c on rerentre les masses |
---|
227 | enddo |
---|
228 | do i=1,iim |
---|
229 | sy(i,1,l)=sy(i,1,l)*sm(i,1,l) |
---|
230 | sy(i,jjp1,l)=sy(i,jjp1,l)*sm(i,jjp1,l) |
---|
231 | s0(i,1,l)=s0(i,1,l)*sm(i,1,l) |
---|
232 | s0(i,jjp1,l)=s0(i,jjp1,l)*sm(i,jjp1,l) |
---|
233 | enddo |
---|
234 | sxn(iip1)=sxn(1) |
---|
235 | sxs(iip1)=sxs(1) |
---|
236 | do i=1,iim |
---|
237 | sx(i+1,1,l)=0.25*(sxn(i)+sxn(i+1))*sm(i+1,1,l) |
---|
238 | sx(i+1,jjp1,l)=0.25*(sxs(i)+sxs(i+1))*sm(i+1,jjp1,l) |
---|
239 | enddo |
---|
240 | s0(iip1,1,l)=s0(1,1,l) |
---|
241 | s0(iip1,jjp1,l)=s0(1,jjp1,l) |
---|
242 | sy(iip1,1,l)=sy(1,1,l) |
---|
243 | sy(iip1,jjp1,l)=sy(1,jjp1,l) |
---|
244 | sx(1,1,l)=sx(iip1,1,l) |
---|
245 | sx(1,jjp1,l)=sx(iip1,jjp1,l) |
---|
246 | enddo |
---|
247 | endif |
---|
248 | |
---|
249 | if (mode.eq.4) then |
---|
250 | do l=1,llm |
---|
251 | do i=1,iip1 |
---|
252 | sx(i,1,l)=0. |
---|
253 | sx(i,jjp1,l)=0. |
---|
254 | sy(i,1,l)=0. |
---|
255 | sy(i,jjp1,l)=0. |
---|
256 | enddo |
---|
257 | enddo |
---|
258 | endif |
---|
259 | call limx(s0,sx,sm,pente_max) |
---|
260 | c call minmaxq(zq,1.e33,-1.e33,'avant advx ') |
---|
261 | call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf) |
---|
262 | c call minmaxq(zq,1.e33,-1.e33,'avant advy ') |
---|
263 | if (mode.eq.4) then |
---|
264 | do l=1,llm |
---|
265 | do i=1,iip1 |
---|
266 | sx(i,1,l)=0. |
---|
267 | sx(i,jjp1,l)=0. |
---|
268 | sy(i,1,l)=0. |
---|
269 | sy(i,jjp1,l)=0. |
---|
270 | enddo |
---|
271 | enddo |
---|
272 | endif |
---|
273 | call limy(s0,sy,sm,pente_max) |
---|
274 | call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz ) |
---|
275 | c call minmaxq(zq,1.e33,-1.e33,'avant advz ') |
---|
276 | do j=1,jjp1 |
---|
277 | do i=1,iip1 |
---|
278 | sz(i,j,1)=0. |
---|
279 | sz(i,j,llm)=0. |
---|
280 | enddo |
---|
281 | enddo |
---|
282 | call limz(s0,sz,sm,pente_max) |
---|
283 | call advz( limit,dtvr,w,sm,s0,sx,sy,sz ) |
---|
284 | if (mode.eq.4) then |
---|
285 | do l=1,llm |
---|
286 | do i=1,iip1 |
---|
287 | sx(i,1,l)=0. |
---|
288 | sx(i,jjp1,l)=0. |
---|
289 | sy(i,1,l)=0. |
---|
290 | sy(i,jjp1,l)=0. |
---|
291 | enddo |
---|
292 | enddo |
---|
293 | endif |
---|
294 | call limy(s0,sy,sm,pente_max) |
---|
295 | call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz ) |
---|
296 | do l=1,llm |
---|
297 | do j=1,jjp1 |
---|
298 | sm(iip1,j,l)=sm(1,j,l) |
---|
299 | s0(iip1,j,l)=s0(1,j,l) |
---|
300 | sx(iip1,j,l)=sx(1,j,l) |
---|
301 | sy(iip1,j,l)=sy(1,j,l) |
---|
302 | sz(iip1,j,l)=sz(1,j,l) |
---|
303 | enddo |
---|
304 | enddo |
---|
305 | |
---|
306 | |
---|
307 | c call minmaxq(zq,1.e33,-1.e33,'avant advx ') |
---|
308 | if (mode.eq.4) then |
---|
309 | do l=1,llm |
---|
310 | do i=1,iip1 |
---|
311 | sx(i,1,l)=0. |
---|
312 | sx(i,jjp1,l)=0. |
---|
313 | sy(i,1,l)=0. |
---|
314 | sy(i,jjp1,l)=0. |
---|
315 | enddo |
---|
316 | enddo |
---|
317 | endif |
---|
318 | call limx(s0,sx,sm,pente_max) |
---|
319 | call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf) |
---|
320 | c call minmaxq(zq,1.e33,-1.e33,'apres advx ') |
---|
321 | c do l=1,llm |
---|
322 | c do j=1,jjp1 |
---|
323 | c do i=1,iip1 |
---|
324 | c zq=s0(i,j,l)/sm(i,j,l) |
---|
325 | c if(zq.lt.qmin) |
---|
326 | c , print*,'apres advx2, s0(',i,',',j,',',l,')=',zq |
---|
327 | c enddo |
---|
328 | c enddo |
---|
329 | c enddo |
---|
330 | c *** On repasse les S dans la variable q directement 14/10/94 |
---|
331 | c On revient a des rapports de melange en divisant par la masse |
---|
332 | |
---|
333 | c En dehors des poles: |
---|
334 | |
---|
335 | DO l = 1,llm |
---|
336 | DO j = 1,jjp1 |
---|
337 | DO i = 1,iim |
---|
338 | q(i,j,llm+1-l,0)=s0(i,j,l)/sm(i,j,l) |
---|
339 | q(i,j,llm+1-l,1)=sx(i,j,l)/sm(i,j,l) |
---|
340 | q(i,j,llm+1-l,2)=sy(i,j,l)/sm(i,j,l) |
---|
341 | q(i,j,llm+1-l,3)=sz(i,j,l)/sm(i,j,l) |
---|
342 | ENDDO |
---|
343 | ENDDO |
---|
344 | ENDDO |
---|
345 | |
---|
346 | c Traitements specifiques au pole |
---|
347 | |
---|
348 | if(mode.ge.1) then |
---|
349 | DO l=1,llm |
---|
350 | c filtrages aux poles |
---|
351 | masn=ssum(iim,sm(1,1,l),1) |
---|
352 | mass=ssum(iim,sm(1,jjp1,l),1) |
---|
353 | qpn=ssum(iim,s0(1,1,l),1)/masn |
---|
354 | qps=ssum(iim,s0(1,jjp1,l),1)/mass |
---|
355 | dqzpn=ssum(iim,sz(1,1,l),1)/masn |
---|
356 | dqzps=ssum(iim,sz(1,jjp1,l),1)/mass |
---|
357 | do i=1,iip1 |
---|
358 | q( i,1,llm+1-l,3)=dqzpn |
---|
359 | q( i,jjp1,llm+1-l,3)=dqzps |
---|
360 | q( i,1,llm+1-l,0)=qpn |
---|
361 | q( i,jjp1,llm+1-l,0)=qps |
---|
362 | enddo |
---|
363 | if(mode.eq.3) then |
---|
364 | dyn1=0. |
---|
365 | dys1=0. |
---|
366 | dyn2=0. |
---|
367 | dys2=0. |
---|
368 | do i=1,iim |
---|
369 | dyn1=dyn1+sinlondlon(i)*sy(i,1,l)/sm(i,1,l) |
---|
370 | dyn2=dyn2+coslondlon(i)*sy(i,1,l)/sm(i,1,l) |
---|
371 | dys1=dys1+sinlondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l) |
---|
372 | dys2=dys2+coslondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l) |
---|
373 | enddo |
---|
374 | do i=1,iim |
---|
375 | q(i,1,llm+1-l,2)= |
---|
376 | s (sinlon(i)*dyn1+coslon(i)*dyn2) |
---|
377 | q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2) |
---|
378 | q(i,jjp1,llm+1-l,2)= |
---|
379 | s (sinlon(i)*dys1+coslon(i)*dys2) |
---|
380 | q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0) |
---|
381 | s -q(i,jjp1,llm+1-l,2) |
---|
382 | enddo |
---|
383 | endif |
---|
384 | if(mode.eq.1) then |
---|
385 | c on filtre les valeurs au bord de la "grande maille pole" |
---|
386 | dyn1=0. |
---|
387 | dys1=0. |
---|
388 | dyn2=0. |
---|
389 | dys2=0. |
---|
390 | do i=1,iim |
---|
391 | zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0) |
---|
392 | dyn1=dyn1+sinlondlon(i)*zz |
---|
393 | dyn2=dyn2+coslondlon(i)*zz |
---|
394 | zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l) |
---|
395 | dys1=dys1+sinlondlon(i)*zz |
---|
396 | dys2=dys2+coslondlon(i)*zz |
---|
397 | enddo |
---|
398 | do i=1,iim |
---|
399 | q(i,1,llm+1-l,2)= |
---|
400 | s (sinlon(i)*dyn1+coslon(i)*dyn2)/2. |
---|
401 | q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2) |
---|
402 | q(i,jjp1,llm+1-l,2)= |
---|
403 | s (sinlon(i)*dys1+coslon(i)*dys2)/2. |
---|
404 | q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0) |
---|
405 | s -q(i,jjp1,llm+1-l,2) |
---|
406 | enddo |
---|
407 | q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0) |
---|
408 | q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0) |
---|
409 | |
---|
410 | do i=1,iim |
---|
411 | sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0) |
---|
412 | sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0) |
---|
413 | enddo |
---|
414 | sxn(iip1)=sxn(1) |
---|
415 | sxs(iip1)=sxs(1) |
---|
416 | do i=1,iim |
---|
417 | q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1)) |
---|
418 | q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1)) |
---|
419 | enddo |
---|
420 | q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1) |
---|
421 | q(1,jjp1,llm+1-l,1)=q(iip1,jjp1,llm+1-l,1) |
---|
422 | |
---|
423 | endif |
---|
424 | |
---|
425 | ENDDO |
---|
426 | endif |
---|
427 | |
---|
428 | c bouclage en longitude |
---|
429 | do iq=0,3 |
---|
430 | do l=1,llm |
---|
431 | do j=1,jjp1 |
---|
432 | q(iip1,j,l,iq)=q(1,j,l,iq) |
---|
433 | enddo |
---|
434 | enddo |
---|
435 | enddo |
---|
436 | |
---|
437 | c PRINT*, ' SORTIE DE PENTES --- ca peut glisser ....' |
---|
438 | |
---|
439 | DO l = 1,llm |
---|
440 | DO j = 1,jjp1 |
---|
441 | DO i = 1,iip1 |
---|
442 | IF (q(i,j,l,0).lt.0.) THEN |
---|
443 | c PRINT*,'------------ BIP-----------' |
---|
444 | c PRINT*,'Q0(',i,j,l,')=',q(i,j,l,0) |
---|
445 | c PRINT*,'QX(',i,j,l,')=',q(i,j,l,1) |
---|
446 | c PRINT*,'QY(',i,j,l,')=',q(i,j,l,2) |
---|
447 | c PRINT*,'QZ(',i,j,l,')=',q(i,j,l,3) |
---|
448 | c PRINT*,' PBL EN SORTIE DE PENTES' |
---|
449 | q(i,j,l,0)=0. |
---|
450 | c STOP |
---|
451 | ENDIF |
---|
452 | ENDDO |
---|
453 | ENDDO |
---|
454 | ENDDO |
---|
455 | |
---|
456 | c PRINT*, '-------------------------------------------' |
---|
457 | |
---|
458 | do l=1,llm |
---|
459 | do j=1,jjp1 |
---|
460 | do i=1,iip1 |
---|
461 | if(q(i,j,l,0).lt.qmin) |
---|
462 | , print*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0) |
---|
463 | enddo |
---|
464 | enddo |
---|
465 | enddo |
---|
466 | RETURN |
---|
467 | END |
---|
468 | |
---|
469 | |
---|
470 | |
---|
471 | |
---|
472 | |
---|
473 | |
---|
474 | |
---|
475 | |
---|
476 | |
---|
477 | |
---|
478 | |
---|
479 | |
---|