GCC Code Coverage Report


Directory: ./
File: dyn3d_common/prather.f
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 178 0.0%
Branches: 0 82 0.0%

Line Branch Exec Source
1 !
2 ! $Header$
3 !
4 SUBROUTINE prather (q,w,masse,pbaru,pbarv,nt,dt)
5
6 USE comconst_mod, ONLY: pi
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 de prather
16 c Ref :
17 c
18 c ************************************************
19 c q,w,pext,pbaru et pbarv : arguments d'entree pour le s-pg
20 c
21 c=======================================================================
22
23
24 include "dimensions.h"
25 include "paramet.h"
26 include "comgeom2.h"
27
28 c Arguments:
29 c ----------
30 INTEGER iq,nt
31 REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
32 REAL masse(iip1,jjp1,llm)
33 REAL q( iip1,jjp1,llm,0:9)
34 REAL w( ip1jmp1,llm )
35 integer ordre,ilim
36
37 c Local:
38 c ------
39 LOGICAL limit
40 real zq(iip1,jjp1,llm)
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 sxx( iip1,jjp1,llm)
45 REAL sxy( iip1,jjp1,llm)
46 REAL sxz( iip1,jjp1,llm)
47 REAL syy( iip1,jjp1,llm )
48 REAL syz( iip1,jjp1,llm )
49 REAL szz( iip1,jjp1,llm ),zz
50 INTEGER i,j,l,indice
51 real sxn(iip1),sxs(iip1)
52
53 real sinlon(iip1),sinlondlon(iip1)
54 real coslon(iip1),coslondlon(iip1)
55 real qmin,qmax
56 save qmin,qmax
57 save sinlon,coslon,sinlondlon,coslondlon
58 real dyn1,dyn2,dys1,dys2,qpn,qps,dqzpn,dqzps
59 real masn,mass
60 c
61 REAL SSUM
62 integer ismax,ismin
63 EXTERNAL SSUM, ismin,ismax
64 logical first
65 save first
66 EXTERNAL advxp,advyp,advzp
67
68
69 data first/.true./
70 data qmin,qmax/-1.e33,1.e33/
71
72
73 c==========================================================================
74 c==========================================================================
75 c MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
76 c==========================================================================
77 c==========================================================================
78 REAL dt
79 c==========================================================================
80 limit = .TRUE.
81
82 if(first) then
83 print*,'SCHEMA PRATHER'
84 first=.false.
85 do i=2,iip1
86 coslon(i)=cos(rlonv(i))
87 sinlon(i)=sin(rlonv(i))
88 coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
89 sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
90 enddo
91 coslon(1)=coslon(iip1)
92 coslondlon(1)=coslondlon(iip1)
93 sinlon(1)=sinlon(iip1)
94 sinlondlon(1)=sinlondlon(iip1)
95
96 DO l = 1,llm
97 DO j = 1,jjp1
98 DO i = 1,iip1
99 q( i,j,l,1 )=0.
100 q( i,j,l,2)=0.
101 q( i,j,l,3)=0.
102 q( i,j,l,4)=0.
103 q( i,j,l,5)=0.
104 q( i,j,l,6)=0.
105 q( i,j,l,7)=0.
106 q( i,j,l,8)=0.
107 q( i,j,l,9)=0.
108 ENDDO
109 ENDDO
110 ENDDO
111 endif
112 c Fin modif Fred
113
114 c *** On calcule la masse d'air en kg
115
116 DO l = 1,llm
117 DO j = 1,jjp1
118 DO i = 1,iip1
119 sm( i,j,llm+1-l ) =masse(i,j,l)
120 ENDDO
121 ENDDO
122 ENDDO
123
124 c *** q contient les qqtes de traceur avant l'advection
125
126 c *** Affectation des tableaux S a partir de Q
127
128 DO l = 1,llm
129 DO j = 1,jjp1
130 DO i = 1,iip1
131 s0( i,j,l) = q ( i,j,llm+1-l,0 )*sm(i,j,l)
132 sx( i,j,l) = q( i,j,llm+1-l,1 )*sm(i,j,l)
133 sy( i,j,l) = q( i,j,llm+1-l,2)*sm(i,j,l)
134 sz( i,j,l) = q( i,j,llm+1-l,3)*sm(i,j,l)
135 sxx( i,j,l) = q( i,j,llm+1-l,4)*sm(i,j,l)
136 sxy( i,j,l) = q( i,j,llm+1-l,5)*sm(i,j,l)
137 sxz( i,j,l) = q( i,j,llm+1-l,6)*sm(i,j,l)
138 syy( i,j,l) = q( i,j,llm+1-l,7)*sm(i,j,l)
139 syz( i,j,l) = q( i,j,llm+1-l,8)*sm(i,j,l)
140 szz( i,j,l) = q( i,j,llm+1-l,9)*sm(i,j,l)
141 ENDDO
142 ENDDO
143 ENDDO
144 c *** Appel des subroutines d'advection en X, en Y et en Z
145 c *** Advection avec "time-splitting"
146
147 c-----------------------------------------------------------
148 do indice =1,nt
149 call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz
150 . ,sxx,sxy,sxz,syy,syz,szz,1 )
151 end do
152 do l=1,llm
153 do i=1,iip1
154 sy(i,1,l)=0.
155 sy(i,jjp1,l)=0.
156 enddo
157 enddo
158 c---------------------------------------------------------
159 call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz
160 . ,sxx,sxy,sxz,syy,syz,szz,1 )
161 c---------------------------------------------------------
162
163 c---------------------------------------------------------
164 do j=1,jjp1
165 do i=1,iip1
166 sz(i,j,1)=0.
167 sz(i,j,llm)=0.
168 sxz(i,j,1)=0.
169 sxz(i,j,llm)=0.
170 syz(i,j,1)=0.
171 syz(i,j,llm)=0.
172 szz(i,j,1)=0.
173 szz(i,j,llm)=0.
174 enddo
175 enddo
176 call advzp( limit,dt*nt,w,sm,s0,sx,sy,sz
177 . ,sxx,sxy,sxz,syy,syz,szz,1 )
178 do l=1,llm
179 do i=1,iip1
180 sy(i,1,l)=0.
181 sy(i,jjp1,l)=0.
182 enddo
183 enddo
184
185 c---------------------------------------------------------
186
187 c---------------------------------------------------------
188 call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz
189 . ,sxx,sxy,sxz,syy,syz,szz,1 )
190 c---------------------------------------------------------
191 DO l = 1,llm
192 DO j = 1,jjp1
193 s0( iip1,j,l)=s0( 1,j,l )
194 sx( iip1,j,l)=sx( 1,j,l )
195 sy( iip1,j,l)=sy( 1,j,l )
196 sz( iip1,j,l)=sz( 1,j,l )
197 sxx( iip1,j,l)=sxx( 1,j,l )
198 sxy( iip1,j,l)=sxy( 1,j,l)
199 sxz( iip1,j,l)=sxz( 1,j,l )
200 syy( iip1,j,l)=syy( 1,j,l )
201 syz( iip1,j,l)=syz( 1,j,l)
202 szz( iip1,j,l)=szz( 1,j,l )
203 ENDDO
204 ENDDO
205 do indice=1,nt
206 call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz
207 . ,sxx,sxy,sxz,syy,syz,szz,1 )
208 end do
209 c---------------------------------------------------------
210 c---------------------------------------------------------
211 c *** On repasse les S dans la variable qpr
212 c *** On repasse les S dans la variable q directement 14/10/94
213
214 DO l = 1,llm
215 DO j = 1,jjp1
216 DO i = 1,iip1
217 q( i,j,llm+1-l,0 )=s0( i,j,l )/sm(i,j,l)
218 q( i,j,llm+1-l,1 ) = sx( i,j,l )/sm(i,j,l)
219 q( i,j,llm+1-l,2 ) = sy( i,j,l )/sm(i,j,l)
220 q( i,j,llm+1-l,3 ) = sz( i,j,l )/sm(i,j,l)
221 q( i,j,llm+1-l,4 ) = sxx( i,j,l )/sm(i,j,l)
222 q( i,j,llm+1-l,5 ) = sxy( i,j,l )/sm(i,j,l)
223 q( i,j,llm+1-l,6 ) = sxz( i,j,l )/sm(i,j,l)
224 q( i,j,llm+1-l,7 ) = syy( i,j,l )/sm(i,j,l)
225 q( i,j,llm+1-l,8 ) = syz( i,j,l )/sm(i,j,l)
226 q( i,j,llm+1-l,9 ) = szz( i,j,l )/sm(i,j,l)
227 ENDDO
228 ENDDO
229 ENDDO
230
231 c---------------------------------------------------------
232 c go to 777
233 c filtrages aux poles
234
235 c Traitements specifiques au pole
236
237 c filtrages aux poles
238 DO l=1,llm
239 c filtrages aux poles
240 masn=ssum(iim,sm(1,1,l),1)
241 mass=ssum(iim,sm(1,jjp1,l),1)
242 qpn=ssum(iim,s0(1,1,l),1)/masn
243 qps=ssum(iim,s0(1,jjp1,l),1)/mass
244 dqzpn=ssum(iim,sz(1,1,l),1)/masn
245 dqzps=ssum(iim,sz(1,jjp1,l),1)/mass
246 do i=1,iip1
247 q( i,1,llm+1-l,3)=dqzpn
248 q( i,jjp1,llm+1-l,3)=dqzps
249 q( i,1,llm+1-l,0)=qpn
250 q( i,jjp1,llm+1-l,0)=qps
251 enddo
252 c enddo
253 c print*,'qpn',qpn,'qps',qps
254 c print*,'dqzpn',dqzpn,'dqzps',dqzps
255 c enddo
256 dyn1=0.
257 dys1=0.
258 dyn2=0.
259 dys2=0.
260 do i=1,iim
261 zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)
262 dyn1=dyn1+sinlondlon(i)*zz
263 dyn2=dyn2+coslondlon(i)*zz
264 zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)
265 dys1=dys1+sinlondlon(i)*zz
266 dys2=dys2+coslondlon(i)*zz
267 enddo
268 do i=1,iim
269 q(i,1,llm+1-l,2)=
270 $ (sinlon(i)*dyn1+coslon(i)*dyn2)/2.
271 q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)
272 $ +q(i,1,llm+1-l,2)
273 q(i,jjp1,llm+1-l,2)=
274 $ (sinlon(i)*dys1+coslon(i)*dys2)/2.
275 q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
276 $ -q(i,jjp1,llm+1-l,2)
277 enddo
278 q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)
279 q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)
280 do i=1,iim
281 sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)
282 sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)
283 enddo
284 sxn(iip1)=sxn(1)
285 sxs(iip1)=sxs(1)
286 do i=1,iim
287 q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))
288 q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))
289 END DO
290 q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)
291 q(1,jjp1,llm+1-l,1)=
292 $ q(iip1,jjp1,llm+1-l,1)
293 enddo
294 do l=1,llm
295 do i=1,iim
296 q( i,1,llm+1-l,4)=0.
297 q( i,jjp1,llm+1-l,4)=0.
298 q( i,1,llm+1-l,5)=0.
299 q( i,jjp1,llm+1-l,5)=0.
300 q( i,1,llm+1-l,6)=0.
301 q( i,jjp1,llm+1-l,6)=0.
302 q( i,1,llm+1-l,7)=0.
303 q( i,jjp1,llm+1-l,7)=0.
304 q( i,1,llm+1-l,8)=0.
305 q( i,jjp1,llm+1-l,8)=0.
306 q( i,1,llm+1-l,9)=0.
307 q( i,jjp1,llm+1-l,9)=0.
308 enddo
309 ENDDO
310
311 777 continue
312 c
313 c bouclage en longitude
314 do l=1,llm
315 do j=1,jjp1
316 q(iip1,j,l,0)=q(1,j,l,0)
317 q(iip1,j,llm+1-l,0)=q(1,j,llm+1-l,0)
318 q(iip1,j,llm+1-l,1)=q(1,j,llm+1-l,1)
319 q(iip1,j,llm+1-l,2)=q(1,j,llm+1-l,2)
320 q(iip1,j,llm+1-l,3)=q(1,j,llm+1-l,3)
321 q(iip1,j,llm+1-l,4)=q(1,j,llm+1-l,4)
322 q(iip1,j,llm+1-l,5)=q(1,j,llm+1-l,5)
323 q(iip1,j,llm+1-l,6)=q(1,j,llm+1-l,6)
324 q(iip1,j,llm+1-l,7)=q(1,j,llm+1-l,7)
325 q(iip1,j,llm+1-l,8)=q(1,j,llm+1-l,8)
326 q(iip1,j,llm+1-l,9)=q(1,j,llm+1-l,9)
327 enddo
328 enddo
329 DO l = 1,llm
330 DO j = 2,jjm
331 DO i = 1,iip1
332 IF (q(i,j,l,0).lt.0.) THEN
333 PRINT*,'------------ BIP-----------'
334 PRINT*,'S0(',i,j,l,')=',q(i,j,l,0),
335 $ q(i,j-1,l,0)
336 PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)
337 PRINT*,'SY(',i,j,l,')=',q(i,j,l,2),
338 $ q(i,j-1,l,2)
339 PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)
340 c PRINT*,' PBL EN SORTIE D'' ADVZP'
341 q(i,j,l,0)=0.
342 c STOP
343 ENDIF
344 ENDDO
345 ENDDO
346 do j=1,jjp1,jjm
347 do i=1,iip1
348 IF (q(i,j,l,0).lt.0.) THEN
349 PRINT*,'------------ BIP 2-----------'
350 PRINT*,'S0(',i,j,l,')=',q(i,j,l,0)
351 PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)
352 PRINT*,'SY(',i,j,l,')=',q(i,j,l,2)
353 PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)
354
355 q(i,j,l,0)=0.
356 c STOP
357 ENDIF
358 enddo
359 enddo
360 ENDDO
361 RETURN
362 END
363