GCC Code Coverage Report


Directory: ./
File: dyn3d_common/pentes_ini.f
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 207 0.0%
Branches: 0 128 0.0%

Line Branch Exec Source
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
480