GCC Code Coverage Report


Directory: ./
File: phys/thermcell_dq.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 35 90 38.9%
Branches: 42 94 44.7%

Line Branch Exec Source
1 2880 subroutine thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr, &
2 & masse,q,dq,qa,lev_out)
3 USE print_control_mod, ONLY: prt_level
4 implicit none
5
6 !=======================================================================
7 !
8 ! Calcul du transport verticale dans la couche limite en presence
9 ! de "thermiques" explicitement representes
10 ! calcul du dq/dt une fois qu'on connait les ascendances
11 !
12 ! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
13 ! Introduction of an implicit computation of vertical advection in
14 ! the environment of thermal plumes in thermcell_dq
15 ! impl = 0 : explicit, 1 : implicit, -1 : old version
16 !
17 !=======================================================================
18
19 integer ngrid,nlay,impl
20
21 real ptimestep
22 real masse(ngrid,nlay),fm(ngrid,nlay+1)
23 real entr(ngrid,nlay)
24 real q(ngrid,nlay)
25 real dq(ngrid,nlay)
26 integer lev_out ! niveau pour les print
27
28 5760 real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
29
30 real zzm
31
32 integer ig,k
33 real cfl
34
35 5760 real qold(ngrid,nlay),fqa(ngrid,nlay+1)
36 integer niter,iter
37 CHARACTER (LEN=20) :: modname='thermcell_dq'
38 CHARACTER (LEN=80) :: abort_message
39
40
41 ! Old explicite scheme
42
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2880 times.
2880 if (impl<=-1) then
43 call thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr, &
44 & masse,q,dq,qa,lev_out)
45 return
46 endif
47
48 ! Calcul du critere CFL pour l'advection dans la subsidence
49 cfl = 0.
50
2/2
✓ Branch 0 taken 112320 times.
✓ Branch 1 taken 2880 times.
115200 do k=1,nlay
51
2/2
✓ Branch 0 taken 111646080 times.
✓ Branch 1 taken 112320 times.
111761280 do ig=1,ngrid
52 111646080 zzm=masse(ig,k)/ptimestep
53 cfl=max(cfl,fm(ig,k)/zzm)
54
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 111646080 times.
111758400 if (entr(ig,k).gt.zzm) then
55 print*,'entr*dt>m,1',k,entr(ig,k)*ptimestep,masse(ig,k)
56 abort_message = 'entr dt > m, 1st'
57 CALL abort_physic (modname,abort_message,1)
58 endif
59 enddo
60 enddo
61
62
4/4
✓ Branch 0 taken 112320 times.
✓ Branch 1 taken 2880 times.
✓ Branch 2 taken 111646080 times.
✓ Branch 3 taken 112320 times.
111761280 qold=q
63
64
65
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2880 times.
2880 if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
66
67 ! calcul du detrainement
68
2/2
✓ Branch 0 taken 112320 times.
✓ Branch 1 taken 2880 times.
115200 do k=1,nlay
69
2/2
✓ Branch 0 taken 111646080 times.
✓ Branch 1 taken 112320 times.
111761280 do ig=1,ngrid
70 111646080 detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
71 ! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
72 !test
73
2/2
✓ Branch 0 taken 48966 times.
✓ Branch 1 taken 111597114 times.
111646080 if (detr(ig,k).lt.0.) then
74 48966 entr(ig,k)=entr(ig,k)-detr(ig,k)
75 48966 detr(ig,k)=0.
76 ! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
77 ! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
78 endif
79 if (fm(ig,k+1).lt.0.) then
80 ! print*,'fm2<0!!!'
81 endif
82 112320 if (entr(ig,k).lt.0.) then
83 ! print*,'entr2<0!!!'
84 endif
85 enddo
86 enddo
87
88 ! Computation of tracer concentrations in the ascending plume
89
2/2
✓ Branch 0 taken 2862720 times.
✓ Branch 1 taken 2880 times.
2865600 do ig=1,ngrid
90 2865600 qa(ig,1)=q(ig,1)
91 enddo
92
93
2/2
✓ Branch 0 taken 109440 times.
✓ Branch 1 taken 2880 times.
112320 do k=2,nlay
94
2/2
✓ Branch 0 taken 108783360 times.
✓ Branch 1 taken 109440 times.
108895680 do ig=1,ngrid
95
2/2
✓ Branch 0 taken 7820046 times.
✓ Branch 1 taken 100963314 times.
108783360 if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
96 & 1.e-5*masse(ig,k)) then
97 qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k)) &
98 7820046 & /(fm(ig,k+1)+detr(ig,k))
99 else
100 100963314 qa(ig,k)=q(ig,k)
101 endif
102 if (qa(ig,k).lt.0.) then
103 ! print*,'qa<0!!!'
104 endif
105 109440 if (q(ig,k).lt.0.) then
106 ! print*,'q<0!!!'
107 endif
108 enddo
109 enddo
110
111 ! Plume vertical flux
112
2/2
✓ Branch 0 taken 106560 times.
✓ Branch 1 taken 2880 times.
109440 do k=2,nlay-1
113
2/2
✓ Branch 0 taken 105920640 times.
✓ Branch 1 taken 106560 times.
106030080 fqa(:,k)=fm(:,k)*qa(:,k-1)
114 enddo
115
4/4
✓ Branch 0 taken 2862720 times.
✓ Branch 1 taken 2880 times.
✓ Branch 2 taken 2862720 times.
✓ Branch 3 taken 2880 times.
5728320 fqa(:,1)=0. ; fqa(:,nlay)=0.
116
117
118 ! Trace species evolution
119
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2880 times.
2880 if (impl==0) then
120 do k=1,nlay-1
121 q(:,k)=q(:,k)+(fqa(:,k)-fqa(:,k+1)-fm(:,k)*q(:,k)+fm(:,k+1)*q(:,k+1)) &
122 & *ptimestep/masse(:,k)
123 enddo
124 else
125
2/2
✓ Branch 0 taken 2880 times.
✓ Branch 1 taken 109440 times.
112320 do k=nlay-1,1,-1
126 ! FH debut de modif : le calcul ci dessous modifiait numériquement
127 ! la concentration quand le flux de masse etait nul car on divisait
128 ! puis multipliait par masse/ptimestep.
129 ! q(:,k)=(masse(:,k)*q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1)) &
130 ! & /(fm(:,k)+masse(:,k)/ptimestep)
131 q(:,k)=(q(:,k)+ptimestep/masse(:,k)*(fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1))) &
132
2/2
✓ Branch 0 taken 108783360 times.
✓ Branch 1 taken 109440 times.
108895680 & /(1.+fm(:,k)*ptimestep/masse(:,k))
133 ! FH fin de modif.
134 enddo
135 endif
136
137 ! Tendencies
138
2/2
✓ Branch 0 taken 2880 times.
✓ Branch 1 taken 112320 times.
115200 do k=1,nlay
139
2/2
✓ Branch 0 taken 111646080 times.
✓ Branch 1 taken 112320 times.
111761280 do ig=1,ngrid
140 111646080 dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
141 111758400 q(ig,k)=qold(ig,k)
142 enddo
143 enddo
144
145 return
146 end
147
148 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
149 ! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
150 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
151
152 subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr, &
153 & masse,q,dq,qa,lev_out)
154 USE print_control_mod, ONLY: prt_level
155 implicit none
156
157 !=======================================================================
158 !
159 ! Calcul du transport verticale dans la couche limite en presence
160 ! de "thermiques" explicitement representes
161 ! calcul du dq/dt une fois qu'on connait les ascendances
162 !
163 !=======================================================================
164
165 integer ngrid,nlay,impl
166
167 real ptimestep
168 real masse(ngrid,nlay),fm(ngrid,nlay+1)
169 real entr(ngrid,nlay)
170 real q(ngrid,nlay)
171 real dq(ngrid,nlay)
172 integer lev_out ! niveau pour les print
173
174 real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
175
176 real zzm
177
178 integer ig,k
179 real cfl
180
181 real qold(ngrid,nlay)
182 real ztimestep
183 integer niter,iter
184 CHARACTER (LEN=20) :: modname='thermcell_dq'
185 CHARACTER (LEN=80) :: abort_message
186
187
188
189 ! Calcul du critere CFL pour l'advection dans la subsidence
190 cfl = 0.
191 do k=1,nlay
192 do ig=1,ngrid
193 zzm=masse(ig,k)/ptimestep
194 cfl=max(cfl,fm(ig,k)/zzm)
195 if (entr(ig,k).gt.zzm) then
196 print*,'entr*dt>m,2',k,entr(ig,k)*ptimestep,masse(ig,k)
197 abort_message = 'entr dt > m, 2nd'
198 CALL abort_physic (modname,abort_message,1)
199 endif
200 enddo
201 enddo
202
203 !IM 090508 print*,'CFL CFL CFL CFL ',cfl
204
205 niter=1
206
207 ztimestep=ptimestep/niter
208 qold=q
209
210
211 do iter=1,niter
212 if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
213
214 ! calcul du detrainement
215 do k=1,nlay
216 do ig=1,ngrid
217 detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
218 ! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
219 !test
220 if (detr(ig,k).lt.0.) then
221 entr(ig,k)=entr(ig,k)-detr(ig,k)
222 detr(ig,k)=0.
223 ! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
224 ! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
225 endif
226 if (fm(ig,k+1).lt.0.) then
227 ! print*,'fm2<0!!!'
228 endif
229 if (entr(ig,k).lt.0.) then
230 ! print*,'entr2<0!!!'
231 endif
232 enddo
233 enddo
234
235 ! calcul de la valeur dans les ascendances
236 do ig=1,ngrid
237 qa(ig,1)=q(ig,1)
238 enddo
239
240 do k=2,nlay
241 do ig=1,ngrid
242 if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt. &
243 & 1.e-5*masse(ig,k)) then
244 qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k)) &
245 & /(fm(ig,k+1)+detr(ig,k))
246 else
247 qa(ig,k)=q(ig,k)
248 endif
249 if (qa(ig,k).lt.0.) then
250 ! print*,'qa<0!!!'
251 endif
252 if (q(ig,k).lt.0.) then
253 ! print*,'q<0!!!'
254 endif
255 enddo
256 enddo
257
258 ! Calcul du flux subsident
259
260 do k=2,nlay
261 do ig=1,ngrid
262
263 ! Schema avec advection sur plus qu'une maille.
264 zzm=masse(ig,k)/ztimestep
265 if (fm(ig,k)>zzm) then
266 wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1)
267 else
268 wqd(ig,k)=fm(ig,k)*q(ig,k)
269 endif
270
271 if (wqd(ig,k).lt.0.) then
272 ! print*,'wqd<0!!!'
273 endif
274 enddo
275 enddo
276 do ig=1,ngrid
277 wqd(ig,1)=0.
278 wqd(ig,nlay+1)=0.
279 enddo
280
281
282 ! Calcul des tendances
283 do k=1,nlay
284 do ig=1,ngrid
285 q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k) &
286 & -wqd(ig,k)+wqd(ig,k+1)) &
287 & *ztimestep/masse(ig,k)
288 ! if (dq(ig,k).lt.0.) then
289 ! print*,'dq<0!!!'
290 ! endif
291 enddo
292 enddo
293
294
295 enddo
296
297
298 ! Calcul des tendances
299 do k=1,nlay
300 do ig=1,ngrid
301 dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
302 q(ig,k)=qold(ig,k)
303 enddo
304 enddo
305
306 return
307 end
308