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