GCC Code Coverage Report

Directory: ./
File: phys/cv3_routines.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 1070 1721 62.2%
Branches: 921 1656 55.6%
Line Branch Exec Source
1
2 ! $Id: cv3_routines.F90 3670 2020-04-27 08:49:09Z jyg $
3
4
5
6
7 480 SUBROUTINE cv3_param(nd, k_upper, delt)
8
9 USE ioipsl_getin_p_mod, ONLY : getin_p
10 use mod_phys_lmdz_para
11 IMPLICIT NONE
12
13 !------------------------------------------------------------
14 !Set parameters for convectL for iflag_con = 3
15 !------------------------------------------------------------
16
17
18 !*** PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE ***
19 !*** PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO ***
20 !*** PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. ***
21 !*** EFFICIENCY IS ASSUMED TO BE UNITY ***
22 !*** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT ***
23 !*** SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE ***
24 !*** OF CLOUD ***
25
26 ![TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA]
27 !*** ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF ***
28 !*** APPROACH TO QUASI-EQUILIBRIUM ***
29 !*** (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) ***
30 !*** (BETA MUST BE LESS THAN OR EQUAL TO 1) ***
31
32 !*** DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE ***
33 !*** APPROACH TO QUASI-EQUILIBRIUM ***
34 !*** IT MUST BE LESS THAN 0 ***
35
36 include "cv3param.h"
37 include "cvflag.h"
38 include "conema3.h"
39
40 INTEGER, INTENT(IN) :: nd
41 INTEGER, INTENT(IN) :: k_upper
42 REAL, INTENT(IN) :: delt ! timestep (seconds)
43
44 ! Local variables
45 CHARACTER (LEN=20) :: modname = 'cv3_param'
46 CHARACTER (LEN=80) :: abort_message
47
48 LOGICAL, SAVE :: first = .TRUE.
49 !$OMP THREADPRIVATE(first)
50
51 !glb noff: integer limit for convection (nd-noff)
52 ! minorig: First level of convection
53
54 ! -- limit levels for convection:
55
56 !jyg<
57 ! noff is chosen such that nl = k_upper so that upmost loops end at about 22 km
58 !
59 240 noff = min(max(nd-k_upper, 1), (nd+1)/2)
60 !! noff = 1
61 !>jyg
62 240 minorig = 1
63 240 nl = nd - noff
64 240 nlp = nl + 1
65 240 nlm = nl - 1
66
67
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 239 times.
 240 IF (first) THEN
68 ! -- "microphysical" parameters:
69 ! IM beg: ajout fis. reglage ep
70 ! CR+JYG: shedding coefficient (used when iflag_mix_adiab=1)
71 ! IM lu dans physiq.def via conf_phys.F90 epmax = 0.993
72
73 1 omtrain = 45.0 ! used also for snow (no disctinction rain/snow)
74 ! -- misc:
75 1 dtovsh = -0.2 ! dT for overshoot
76 ! cc dttrig = 5. ! (loose) condition for triggering
77 1 dttrig = 10. ! (loose) condition for triggering
78 1 dtcrit = -2.0
79 ! -- end of convection
80 ! -- interface cloud parameterization:
81 1 delta = 0.01 ! cld
82 ! -- interface with boundary-layer (gust factor): (sb)
83 1 betad = 10.0 ! original value (from convect 4.3)
84
85 ! Var interm pour le getin
86 1 cv_flag_feed=1
87 1 CALL getin_p('cv_flag_feed',cv_flag_feed)
88 1 T_top_max = 1000.
89 1 CALL getin_p('t_top_max',T_top_max)
90 1 dpbase=-40.
91 1 CALL getin_p('dpbase',dpbase)
92 1 pbcrit=150.0
93 1 CALL getin_p('pbcrit',pbcrit)
94 1 ptcrit=500.0
95 1 CALL getin_p('ptcrit',ptcrit)
96 1 sigdz=0.01
97 1 CALL getin_p('sigdz',sigdz)
98 1 spfac=0.15
99 1 CALL getin_p('spfac',spfac)
100 1 tau=8000.
101 1 CALL getin_p('tau',tau)
102 1 flag_wb=1
103 1 CALL getin_p('flag_wb',flag_wb)
104 1 wbmax=6.
105 1 CALL getin_p('wbmax',wbmax)
106 1 ok_convstop=.False.
107 1 CALL getin_p('ok_convstop',ok_convstop)
108 1 tau_stop=15000.
109 1 CALL getin_p('tau_stop',tau_stop)
110 1 ok_intermittent=.False.
111 1 CALL getin_p('ok_intermittent',ok_intermittent)
112 1 ok_optim_yield=.False.
113 1 CALL getin_p('ok_optim_yield',ok_optim_yield)
114 1 ok_homo_tend=.TRUE.
115 1 CALL getin_p('ok_homo_tend',ok_homo_tend)
116 1 ok_entrain=.TRUE.
117 1 CALL getin_p('ok_entrain',ok_entrain)
118
119 1 coef_peel=0.25
120 1 CALL getin_p('coef_peel',coef_peel)
121
122 1 flag_epKEorig=1
123 1 CALL getin_p('flag_epKEorig',flag_epKEorig)
124 1 elcrit=0.0003
125 1 CALL getin_p('elcrit',elcrit)
126 1 tlcrit=-55.0
127 1 CALL getin_p('tlcrit',tlcrit)
128 1 ejectliq=0.
129 1 CALL getin_p('ejectliq',ejectliq)
130 1 ejectice=0.
131 1 CALL getin_p('ejectice',ejectice)
132 1 cvflag_prec_eject = .FALSE.
133 1 CALL getin_p('cvflag_prec_eject',cvflag_prec_eject)
134 1 qsat_depends_on_qt = .FALSE.
135 1 CALL getin_p('qsat_depends_on_qt',qsat_depends_on_qt)
136 1 adiab_ascent_mass_flux_depends_on_ejectliq = .FALSE.
137 1 CALL getin_p('adiab_ascent_mass_flux_depends_on_ejectliq',adiab_ascent_mass_flux_depends_on_ejectliq)
138 1 keepbug_ice_frac = .TRUE.
139 1 CALL getin_p('keepbug_ice_frac', keepbug_ice_frac)
140
141 1 WRITE (*, *) 't_top_max=', t_top_max
142 1 WRITE (*, *) 'dpbase=', dpbase
143 1 WRITE (*, *) 'pbcrit=', pbcrit
144 1 WRITE (*, *) 'ptcrit=', ptcrit
145 1 WRITE (*, *) 'sigdz=', sigdz
146 1 WRITE (*, *) 'spfac=', spfac
147 1 WRITE (*, *) 'tau=', tau
148 1 WRITE (*, *) 'flag_wb=', flag_wb
149 1 WRITE (*, *) 'wbmax=', wbmax
150 1 WRITE (*, *) 'ok_convstop=', ok_convstop
151 1 WRITE (*, *) 'tau_stop=', tau_stop
152 1 WRITE (*, *) 'ok_intermittent=', ok_intermittent
153 1 WRITE (*, *) 'ok_optim_yield =', ok_optim_yield
154 1 WRITE (*, *) 'coef_peel=', coef_peel
155
156 1 WRITE (*, *) 'flag_epKEorig=', flag_epKEorig
157 1 WRITE (*, *) 'elcrit=', elcrit
158 1 WRITE (*, *) 'tlcrit=', tlcrit
159 1 WRITE (*, *) 'ejectliq=', ejectliq
160 1 WRITE (*, *) 'ejectice=', ejectice
161 1 WRITE (*, *) 'cvflag_prec_eject =', cvflag_prec_eject
162 1 WRITE (*, *) 'qsat_depends_on_qt =', qsat_depends_on_qt
163 1 WRITE (*, *) 'adiab_ascent_mass_flux_depends_on_ejectliq =', adiab_ascent_mass_flux_depends_on_ejectliq
164 1 WRITE (*, *) 'keepbug_ice_frac =', keepbug_ice_frac
165
166 1 first = .FALSE.
167 END IF ! (first)
168
169 240 beta = 1.0 - delt/tau
170 alpha1 = 1.5E-3
171 !JYG Correction bug alpha
172 240 alpha1 = alpha1*1.5
173 240 alpha = alpha1*delt/tau
174 !JYG Bug
175 ! cc increase alpha to compensate W decrease:
176 ! c alpha = alpha*1.5
177
178 240 noconv_stop = max(2.,tau_stop/delt)
179
180 240 RETURN
181 END SUBROUTINE cv3_param
182
183 240 SUBROUTINE cv3_incrcount(len, nd, delt, sig)
184
185 IMPLICIT NONE
186
187 ! =====================================================================
188 ! Increment the counter sig(nd)
189 ! =====================================================================
190
191 include "cv3param.h"
192 include "cvflag.h"
193
194 !inputs:
195 INTEGER, INTENT(IN) :: len
196 INTEGER, INTENT(IN) :: nd
197 REAL, INTENT(IN) :: delt ! timestep (seconds)
198
199 !input/output
200 REAL, DIMENSION(len,nd), INTENT(INOUT) :: sig
201
202 !local variables
203 INTEGER il
204
205 ! print *,'cv3_incrcount : noconv_stop ',noconv_stop
206 ! print *,'cv3_incrcount in, sig(1,nd) ',sig(1,nd)
207
1/2
✓ Branch 0 taken 240 times.
✗ Branch 1 not taken.
 240 IF(ok_convstop) THEN
208 ✗ DO il = 1, len
209 ✗ sig(il, nd) = sig(il, nd) + 1.
210 ✗ sig(il, nd) = min(sig(il,nd), noconv_stop+0.1)
211 END DO
212 ELSE
213
2/2
✓ Branch 0 taken 240 times.
✓ Branch 1 taken 238560 times.
 238800 DO il = 1, len
214 238560 sig(il, nd) = sig(il, nd) + 1.
215 238800 sig(il, nd) = min(sig(il,nd), 12.1)
216 END DO
217 ENDIF ! (ok_convstop)
218 ! print *,'cv3_incrcount out, sig(1,nd) ',sig(1,nd)
219
220 240 RETURN
221 END SUBROUTINE cv3_incrcount
222
223 25764960 SUBROUTINE cv3_prelim(len, nd, ndp1, t, q, p, ph, &
224 480 lv, lf, cpn, tv, gz, h, hm, th)
225 IMPLICIT NONE
226
227 ! =====================================================================
228 ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
229 ! "ori": from convect4.3 (vectorized)
230 ! "convect3": to be exactly consistent with convect3
231 ! =====================================================================
232
233 ! inputs:
234 INTEGER len, nd, ndp1
235 REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1)
236
237 ! outputs:
238 REAL lv(len, nd), lf(len, nd), cpn(len, nd), tv(len, nd)
239 REAL gz(len, nd), h(len, nd), hm(len, nd)
240 REAL th(len, nd)
241
242 ! local variables:
243 INTEGER k, i
244 REAL rdcp
245 REAL tvx, tvy ! convect3
246 960 REAL cpx(len, nd)
247
248 include "cvthermo.h"
249 include "cv3param.h"
250
251
252 ! ori do 110 k=1,nlp
253 ! abderr do 110 k=1,nl ! convect3
254
2/2
✓ Branch 0 taken 13440 times.
✓ Branch 1 taken 480 times.
 13920 DO k = 1, nlp
255
256
2/2
✓ Branch 0 taken 13359360 times.
✓ Branch 1 taken 13440 times.
 13373280 DO i = 1, len
257 ! debug lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
258 13359360 lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15)
259 !! lf(i, k) = lf0 - clmci*(t(i,k)-273.15) ! erreur de signe !!
260 13359360 lf(i, k) = lf0 + clmci*(t(i,k)-273.15)
261 13359360 cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
262 13359360 cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k)
263 ! ori tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1)
264 13359360 tv(i, k) = t(i, k)*(1.0+q(i,k)/eps-q(i,k))
265 13359360 rdcp = (rrd*(1.-q(i,k))+q(i,k)*rrv)/cpn(i, k)
266 13372800 th(i, k) = t(i, k)*(1000.0/p(i,k))**rdcp
267 END DO
268 END DO
269
270 ! gz = phi at the full levels (same as p).
271
272 !! DO i = 1, len !jyg
273 !! gz(i, 1) = 0.0 !jyg
274 !! END DO !jyg
275
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 gz(:,:) = 0. !jyg: initialization of the whole array
276 ! ori do 140 k=2,nlp
277
2/2
✓ Branch 0 taken 12480 times.
✓ Branch 1 taken 480 times.
 12960 DO k = 2, nl ! convect3
278
2/2
✓ Branch 0 taken 12405120 times.
✓ Branch 1 taken 12480 times.
 12418080 DO i = 1, len
279 12405120 tvx = t(i, k)*(1.+q(i,k)/eps-q(i,k)) !convect3
280 12405120 tvy = t(i, k-1)*(1.+q(i,k-1)/eps-q(i,k-1)) !convect3
281 gz(i, k) = gz(i, k-1) + 0.5*rrd*(tvx+tvy)* & !convect3
282 12417600 (p(i,k-1)-p(i,k))/ph(i, k) !convect3
283
284 ! c print *,' gz(',k,')',gz(i,k),' tvx',tvx,' tvy ',tvy
285
286 ! ori gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
287 ! ori & *(p(i,k-1)-p(i,k))/ph(i,k)
288 END DO
289 END DO
290
291 ! h = phi + cpT (dry static energy).
292 ! hm = phi + cp(T-Tbase)+Lq
293
294 ! ori do 170 k=1,nlp
295
2/2
✓ Branch 0 taken 12960 times.
✓ Branch 1 taken 480 times.
 13440 DO k = 1, nl ! convect3
296
2/2
✓ Branch 0 taken 12882240 times.
✓ Branch 1 taken 12960 times.
 12895680 DO i = 1, len
297 12882240 h(i, k) = gz(i, k) + cpn(i, k)*t(i, k)
298 12895200 hm(i, k) = gz(i, k) + cpx(i, k)*(t(i,k)-t(i,1)) + lv(i, k)*q(i, k)
299 END DO
300 END DO
301
302 480 RETURN
303 END SUBROUTINE cv3_prelim
304
305 240 SUBROUTINE cv3_feed(len, nd, ok_conserv_q, &
306 240 t, q, u, v, p, ph, h, gz, &
307 p1feed, p2feed, wght, &
308 wghti, tnk, thnk, qnk, qsnk, unk, vnk, &
309 240 cpnk, hnk, nk, icb, icbmax, iflag, gznk, plcl)
310
311 USE mod_phys_lmdz_transfert_para, ONLY : bcast
312 USE add_phys_tend_mod, ONLY: fl_cor_ebil
313 USE print_control_mod, ONLY: prt_level
314 IMPLICIT NONE
315
316 ! ================================================================
317 ! Purpose: CONVECTIVE FEED
318
319 ! Main differences with cv_feed:
320 ! - ph added in input
321 ! - here, nk(i)=minorig
322 ! - icb defined differently (plcl compared with ph instead of p)
323 ! - dry static energy as argument instead of moist static energy
324
325 ! Main differences with convect3:
326 ! - we do not compute dplcldt and dplcldr of CLIFT anymore
327 ! - values iflag different (but tests identical)
328 ! - A,B explicitely defined (!...)
329 ! ================================================================
330
331 include "cv3param.h"
332 include "cvthermo.h"
333
334 !inputs:
335 INTEGER, INTENT (IN) :: len, nd
336 LOGICAL, INTENT (IN) :: ok_conserv_q
337 REAL, DIMENSION (len, nd), INTENT (IN) :: t, q, p
338 REAL, DIMENSION (len, nd), INTENT (IN) :: u, v
339 REAL, DIMENSION (len, nd), INTENT (IN) :: h, gz
340 REAL, DIMENSION (len, nd+1), INTENT (IN) :: ph
341 REAL, DIMENSION (len), INTENT (IN) :: p1feed
342 REAL, DIMENSION (nd), INTENT (IN) :: wght
343 !input-output
344 REAL, DIMENSION (len), INTENT (INOUT) :: p2feed
345 !outputs:
346 INTEGER, INTENT (OUT) :: icbmax
347 INTEGER, DIMENSION (len), INTENT (OUT) :: iflag, nk, icb
348 REAL, DIMENSION (len, nd), INTENT (OUT) :: wghti
349 REAL, DIMENSION (len), INTENT (OUT) :: tnk, thnk, qnk, qsnk
350 REAL, DIMENSION (len), INTENT (OUT) :: unk, vnk
351 REAL, DIMENSION (len), INTENT (OUT) :: cpnk, hnk, gznk
352 REAL, DIMENSION (len), INTENT (OUT) :: plcl
353
354 !local variables:
355 INTEGER i, k, iter, niter
356 INTEGER ihmin(len)
357 REAL work(len)
358 480 REAL pup(len), plo(len), pfeed(len)
359 480 REAL plclup(len), plcllo(len), plclfeed(len)
360 480 REAL pfeedmin(len)
361 480 REAL posit(len)
362 480 LOGICAL nocond(len)
363
364 !jyg20140217<
365 INTEGER iostat
366 LOGICAL, SAVE :: first
367 LOGICAL, SAVE :: ok_new_feed
368 REAL, SAVE :: dp_lcl_feed
369 !$OMP THREADPRIVATE (first,ok_new_feed,dp_lcl_feed)
370 DATA first/.TRUE./
371 DATA dp_lcl_feed/2./
372
373
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 239 times.
 240 IF (first) THEN
374 !$OMP MASTER
375 1 ok_new_feed = ok_conserv_q
376 1 OPEN (98, FILE='cv3feed_param.data', STATUS='old', FORM='formatted', IOSTAT=iostat)
377
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 IF (iostat==0) THEN
378 ✗ READ (98, *, END=998) ok_new_feed
379 998 CONTINUE
380 ✗ CLOSE (98)
381 END IF
382 1 PRINT *, ' ok_new_feed: ', ok_new_feed
383 !$OMP END MASTER
384 1 call bcast(ok_new_feed)
385 1 first = .FALSE.
386 END IF
387 !jyg>
388 ! -------------------------------------------------------------------
389 ! --- Origin level of ascending parcels for convect3:
390 ! -------------------------------------------------------------------
391
392
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
393 238560 nk(i) = minorig
394 238800 gznk(i) = gz(i, nk(i))
395 END DO
396
397 ! -------------------------------------------------------------------
398 ! --- Adjust feeding layer thickness so that lifting up to the top of
399 ! --- the feeding layer does not induce condensation (i.e. so that
400 ! --- plcl < p2feed).
401 ! --- Method : iterative secant method.
402 ! -------------------------------------------------------------------
403
404 ! 1- First bracketing of the solution : ph(nk+1), p2feed
405
406 ! 1.a- LCL associated with p2feed
407
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
408 238800 pup(i) = p2feed(i)
409 END DO
410
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 240 times.
 240 IF (fl_cor_ebil >=2 ) THEN
411 CALL cv3_estatmix(len, nd, iflag, p1feed, pup, p, ph, &
412 t, q, u, v, h, gz, wght, &
413 ✗ wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclup)
414 ELSE
415 CALL cv3_enthalpmix(len, nd, iflag, p1feed, pup, p, ph, &
416 t, q, u, v, wght, &
417 240 wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclup)
418 ENDIF ! (fl_cor_ebil >=2 )
419 ! 1.b- LCL associated with ph(nk+1)
420
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
421 238800 plo(i) = ph(i, nk(i)+1)
422 END DO
423
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 240 times.
 240 IF (fl_cor_ebil >=2 ) THEN
424 CALL cv3_estatmix(len, nd, iflag, p1feed, plo, p, ph, &
425 t, q, u, v, h, gz, wght, &
426 ✗ wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plcllo)
427 ELSE
428 CALL cv3_enthalpmix(len, nd, iflag, p1feed, plo, p, ph, &
429 t, q, u, v, wght, &
430 240 wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plcllo)
431 ENDIF ! (fl_cor_ebil >=2 )
432 ! 2- Iterations
433 niter = 5
434
2/2
✓ Branch 0 taken 1200 times.
✓ Branch 1 taken 240 times.
 1440 DO iter = 1, niter
435
2/2
✓ Branch 0 taken 1192800 times.
✓ Branch 1 taken 1200 times.
 1194000 DO i = 1, len
436 1192800 plcllo(i) = min(plo(i), plcllo(i))
437 1192800 plclup(i) = max(pup(i), plclup(i))
438 1194000 nocond(i) = plclup(i) <= pup(i)
439 END DO
440
2/2
✓ Branch 0 taken 1192800 times.
✓ Branch 1 taken 1200 times.
 1194000 DO i = 1, len
441
2/2
✓ Branch 0 taken 869425 times.
✓ Branch 1 taken 323375 times.
 1194000 IF (nocond(i)) THEN
442 869425 pfeed(i) = pup(i)
443 ELSE
444 !JYG20140217<
445
1/2
✓ Branch 0 taken 323375 times.
✗ Branch 1 not taken.
 323375 IF (ok_new_feed) THEN
446 pfeed(i) = (pup(i)*(plo(i)-plcllo(i)-dp_lcl_feed)+ &
447 plo(i)*(plclup(i)-pup(i)+dp_lcl_feed))/ &
448 323375 (plo(i)-plcllo(i)+plclup(i)-pup(i))
449 ELSE
450 pfeed(i) = (pup(i)*(plo(i)-plcllo(i))+ &
451 plo(i)*(plclup(i)-pup(i)))/ &
452 ✗ (plo(i)-plcllo(i)+plclup(i)-pup(i))
453 END IF
454 !JYG>
455 END IF
456 END DO
457 !jyg20140217<
458 ! For the last iteration, make sure that the top of the feeding layer
459 ! and LCL are not in the same layer:
460
1/2
✓ Branch 0 taken 1200 times.
✗ Branch 1 not taken.
 1200 IF (ok_new_feed) THEN
461
2/2
✓ Branch 0 taken 240 times.
✓ Branch 1 taken 960 times.
 1200 IF (iter==niter) THEN
462
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1,len !jyg
463 238800 pfeedmin(i) = ph(i,minorig+1) !jyg
464 ENDDO !jyg
465
2/2
✓ Branch 0 taken 6240 times.
✓ Branch 1 taken 240 times.
 6480 DO k = minorig+1, nl !jyg
466 !! DO k = minorig, nl !jyg
467
2/2
✓ Branch 0 taken 6202560 times.
✓ Branch 1 taken 6240 times.
 6209040 DO i = 1, len
468
2/2
✓ Branch 0 taken 711123 times.
✓ Branch 1 taken 5491437 times.
 6208800 IF (ph(i,k)>=plclfeed(i)) pfeedmin(i) = ph(i, k)
469 END DO
470 END DO
471
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
472 238800 pfeed(i) = max(pfeedmin(i), pfeed(i))
473 END DO
474 END IF
475 END IF
476 !jyg>
477
478
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1200 times.
 1200 IF (fl_cor_ebil >=2 ) THEN
479 CALL cv3_estatmix(len, nd, iflag, p1feed, pfeed, p, ph, &
480 t, q, u, v, h, gz, wght, &
481 ✗ wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclfeed)
482 ELSE
483 CALL cv3_enthalpmix(len, nd, iflag, p1feed, pfeed, p, ph, &
484 t, q, u, v, wght, &
485 1200 wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclfeed)
486 ENDIF ! (fl_cor_ebil >=2 )
487 !jyg20140217<
488
1/2
✓ Branch 0 taken 1200 times.
✗ Branch 1 not taken.
 1200 IF (ok_new_feed) THEN
489
2/2
✓ Branch 0 taken 1200 times.
✓ Branch 1 taken 1192800 times.
 1194000 DO i = 1, len
490 1192800 posit(i) = (sign(1.,plclfeed(i)-pfeed(i)+dp_lcl_feed)+1.)*0.5
491
2/2
✓ Branch 0 taken 258 times.
✓ Branch 1 taken 1192542 times.
 1194000 IF (plclfeed(i)-pfeed(i)+dp_lcl_feed==0.) posit(i) = 1.
492 END DO
493 ELSE
494 ✗ DO i = 1, len
495 ✗ posit(i) = (sign(1.,plclfeed(i)-pfeed(i))+1.)*0.5
496 ✗ IF (plclfeed(i)==pfeed(i)) posit(i) = 1.
497 END DO
498 END IF
499 !jyg>
500
2/2
✓ Branch 0 taken 1192800 times.
✓ Branch 1 taken 1200 times.
 1194240 DO i = 1, len
501 ! - posit = 1 when lcl is below top of feeding layer (plclfeed>pfeed)
502 ! - => pup=pfeed
503 ! - posit = 0 when lcl is above top of feeding layer (plclfeed<pfeed)
504 ! - => plo=pfeed
505 1192800 pup(i) = posit(i)*pfeed(i) + (1.-posit(i))*pup(i)
506 1192800 plo(i) = (1.-posit(i))*pfeed(i) + posit(i)*plo(i)
507 1192800 plclup(i) = posit(i)*plclfeed(i) + (1.-posit(i))*plclup(i)
508 1194000 plcllo(i) = (1.-posit(i))*plclfeed(i) + posit(i)*plcllo(i)
509 END DO
510 END DO ! iter
511
512
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
513 238560 p2feed(i) = pfeed(i)
514 238800 plcl(i) = plclfeed(i)
515 END DO
516
517
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
518 238560 cpnk(i) = cpd*(1.0-qnk(i)) + cpv*qnk(i)
519 238800 hnk(i) = gz(i, 1) + cpnk(i)*tnk(i)
520 END DO
521
522 ! -------------------------------------------------------------------
523 ! --- Check whether parcel level temperature and specific humidity
524 ! --- are reasonable
525 ! -------------------------------------------------------------------
526
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 240 times.
 240 IF (cv_flag_feed == 1) THEN
527 ✗ DO i = 1, len
528 IF (((tnk(i)<250.0) .OR. &
529 ✗ (qnk(i)<=0.0)) .AND. &
530 ✗ (iflag(i)==0)) iflag(i) = 7
531 END DO
532
1/2
✓ Branch 0 taken 240 times.
✗ Branch 1 not taken.
 240 ELSEIF (cv_flag_feed >= 2) THEN
533 ! --- and demand that LCL be high enough
534
2/2
✓ Branch 0 taken 240 times.
✓ Branch 1 taken 238560 times.
 238800 DO i = 1, len
535 IF (((tnk(i)<250.0) .OR. &
536 (qnk(i)<=0.0) .OR. &
537
6/8
✓ Branch 0 taken 224671 times.
✓ Branch 1 taken 13889 times.
✓ Branch 2 taken 224671 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 68244 times.
✓ Branch 5 taken 156427 times.
✓ Branch 6 taken 82133 times.
✗ Branch 7 not taken.
 238560 (plcl(i)>min(0.99*ph(i,1),ph(i,3)))) .AND. &
538 82373 (iflag(i)==0)) iflag(i) = 7
539 END DO
540 ENDIF
541
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 240 times.
 240 IF (prt_level .GE. 10) THEN
542 ✗ print *,'cv3_feed : iflag(1), pfeed(1), plcl(1), wghti(1,k) ', &
543 ✗ iflag(1), pfeed(1), plcl(1), (wghti(1,k),k=1,10)
544 ENDIF
545
546 ! -------------------------------------------------------------------
547 ! --- Calculate first level above lcl (=icb)
548 ! -------------------------------------------------------------------
549
550 !@ do 270 i=1,len
551 !@ icb(i)=nlm
552 !@ 270 continue
553 !@c
554 !@ do 290 k=minorig,nl
555 !@ do 280 i=1,len
556 !@ if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
557 !@ & icb(i)=min(icb(i),k)
558 !@ 280 continue
559 !@ 290 continue
560 !@c
561 !@ do 300 i=1,len
562 !@ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
563 !@ 300 continue
564
565
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
566 238800 icb(i) = nlm
567 END DO
568
569 ! la modification consiste a comparer plcl a ph et non a p:
570 ! icb est defini par : ph(icb)<plcl<ph(icb-1)
571 !@ do 290 k=minorig,nl
572
2/2
✓ Branch 0 taken 5760 times.
✓ Branch 1 taken 240 times.
 6000 DO k = 3, nl - 1 ! modif pour que icb soit sup/egal a 2
573
2/2
✓ Branch 0 taken 5725440 times.
✓ Branch 1 taken 5760 times.
 5731440 DO i = 1, len
574
2/2
✓ Branch 0 taken 5192126 times.
✓ Branch 1 taken 533314 times.
 5731200 IF (ph(i,k)<plcl(i)) icb(i) = min(icb(i), k)
575 END DO
576 END DO
577
578
579 ! print*,'icb dans cv3_feed '
580 ! write(*,'(64i2)') icb(2:len-1)
581 ! call dump2d(64,43,'plcl dans cv3_feed ',plcl(2:len-1))
582
583
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
584 !@ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
585
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 238560 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 238800 IF ((icb(i)==nlm) .AND. (iflag(i)==0)) iflag(i) = 9
586 END DO
587
588
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
589 238800 icb(i) = icb(i) - 1 ! icb sup ou egal a 2
590 END DO
591
592 ! Compute icbmax.
593
594 240 icbmax = 2
595
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
596 !! icbmax=max(icbmax,icb(i))
597
2/2
✓ Branch 0 taken 156427 times.
✓ Branch 1 taken 82133 times.
 238800 IF (iflag(i)<7) icbmax = max(icbmax, icb(i)) ! sb Jun7th02
598 END DO
599
600 240 RETURN
601 END SUBROUTINE cv3_feed
602
603 3369900 SUBROUTINE cv3_undilute1(len, nd, t, qs, gz, plcl, p, icb, tnk, qnk, gznk, &
604 tp, tvp, clw, icbs)
605 IMPLICIT NONE
606
607 ! ----------------------------------------------------------------
608 ! Equivalent de TLIFT entre NK et ICB+1 inclus
609
610 ! Differences with convect4:
611 ! - specify plcl in input
612 ! - icbs is the first level above LCL (may differ from icb)
613 ! - in the iterations, used x(icbs) instead x(icb)
614 ! - many minor differences in the iterations
615 ! - tvp is computed in only one time
616 ! - icbs: first level above Plcl (IMIN de TLIFT) in output
617 ! - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)
618 ! ----------------------------------------------------------------
619
620 include "cvthermo.h"
621 include "cv3param.h"
622
623 ! inputs:
624 INTEGER, INTENT (IN) :: len, nd
625 INTEGER, DIMENSION (len), INTENT (IN) :: icb
626 REAL, DIMENSION (len, nd), INTENT (IN) :: t, qs, gz
627 REAL, DIMENSION (len), INTENT (IN) :: tnk, qnk, gznk
628 REAL, DIMENSION (len, nd), INTENT (IN) :: p
629 REAL, DIMENSION (len), INTENT (IN) :: plcl ! convect3
630
631 ! outputs:
632 INTEGER, DIMENSION (len), INTENT (OUT) :: icbs
633 REAL, DIMENSION (len, nd), INTENT (OUT) :: tp, tvp, clw
634
635 ! local variables:
636 INTEGER i, k
637 480 INTEGER icb1(len), icbsmax2 ! convect3
638 REAL tg, qg, alv, s, ahg, tc, denom, es, rg
639 480 REAL ah0(len), cpp(len)
640 480 REAL ticb(len), gzicb(len)
641 480 REAL qsicb(len) ! convect3
642 240 REAL cpinv(len) ! convect3
643
644 ! -------------------------------------------------------------------
645 ! --- Calculates the lifted parcel virtual temperature at nk,
646 ! --- the actual temperature, and the adiabatic
647 ! --- liquid water content. The procedure is to solve the equation.
648 ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
649 ! -------------------------------------------------------------------
650
651
652 ! *** Calculate certain parcel quantities, including static energy ***
653
654
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
655 238560 ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
656 238560 cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv
657 238800 cpinv(i) = 1./cpp(i)
658 END DO
659
660 ! *** Calculate lifted parcel quantities below cloud base ***
661
662
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len !convect3
663 238560 icb1(i) = min(max(icb(i), 2), nl)
664 ! if icb is below LCL, start loop at ICB+1:
665 ! (icbs est le premier niveau au-dessus du LCL)
666 238560 icbs(i) = icb1(i) !convect3
667
2/2
✓ Branch 0 taken 92010 times.
✓ Branch 1 taken 146550 times.
 238800 IF (plcl(i)<p(i,icb1(i))) THEN
668 92010 icbs(i) = min(icbs(i)+1, nl) !convect3
669 END IF
670 END DO !convect3
671
672
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len !convect3
673 238560 ticb(i) = t(i, icbs(i)) !convect3
674 238560 gzicb(i) = gz(i, icbs(i)) !convect3
675 238800 qsicb(i) = qs(i, icbs(i)) !convect3
676 END DO !convect3
677
678
679 ! Re-compute icbsmax (icbsmax2): !convect3
680 ! !convect3
681 icbsmax2 = 2 !convect3
682
2/2
✓ Branch 0 taken 240 times.
✓ Branch 1 taken 238560 times.
 238800 DO i = 1, len !convect3
683 238800 icbsmax2 = max(icbsmax2, icbs(i)) !convect3
684 END DO !convect3
685
686 ! initialization outputs:
687
688
2/2
✓ Branch 0 taken 2670 times.
✓ Branch 1 taken 240 times.
 2910 DO k = 1, icbsmax2 ! convect3
689
2/2
✓ Branch 0 taken 2653980 times.
✓ Branch 1 taken 2670 times.
 2656890 DO i = 1, len ! convect3
690 2653980 tp(i, k) = 0.0 ! convect3
691 2653980 tvp(i, k) = 0.0 ! convect3
692 2656650 clw(i, k) = 0.0 ! convect3
693 END DO ! convect3
694 END DO ! convect3
695
696 ! tp and tvp below cloud base:
697
698
2/2
✓ Branch 0 taken 2430 times.
✓ Branch 1 taken 240 times.
 2670 DO k = minorig, icbsmax2 - 1
699
2/2
✓ Branch 0 taken 2415420 times.
✓ Branch 1 taken 2430 times.
 2418090 DO i = 1, len
700 2415420 tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))*cpinv(i)
701 2417850 tvp(i, k) = tp(i, k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)
702 END DO
703 END DO
704
705 ! *** Find lifted parcel quantities above cloud base ***
706
707
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
708 238560 tg = ticb(i)
709 ! ori qg=qs(i,icb(i))
710 238560 qg = qsicb(i) ! convect3
711 ! debug alv=lv0-clmcpv*(ticb(i)-t0)
712 238560 alv = lv0 - clmcpv*(ticb(i)-273.15)
713
714 ! First iteration.
715
716 ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
717 s = cpd*(1.-qnk(i)) + cl*qnk(i) + & ! convect3
718 238560 alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
719 238560 s = 1./s
720 ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
721 238560 ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
722 238560 tg = tg + s*(ah0(i)-ahg)
723 ! ori tg=max(tg,35.0)
724 ! debug tc=tg-t0
725 238560 tc = tg - 273.15
726 238560 denom = 243.5 + tc
727 238560 denom = max(denom, 1.0) ! convect3
728 ! ori if(tc.ge.0.0)then
729 238560 es = 6.112*exp(17.67*tc/denom)
730 ! ori else
731 ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
732 ! ori endif
733 ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
734 238560 qg = eps*es/(p(i,icbs(i))-es*(1.-eps))
735
736 ! Second iteration.
737
738
739 ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
740 ! ori s=1./s
741 ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
742 238560 ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
743 238560 tg = tg + s*(ah0(i)-ahg)
744 ! ori tg=max(tg,35.0)
745 ! debug tc=tg-t0
746 238560 tc = tg - 273.15
747 238560 denom = 243.5 + tc
748 238560 denom = max(denom, 1.0) ! convect3
749 ! ori if(tc.ge.0.0)then
750 238560 es = 6.112*exp(17.67*tc/denom)
751 ! ori else
752 ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
753 ! ori end if
754 ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
755 238560 qg = eps*es/(p(i,icbs(i))-es*(1.-eps))
756
757 alv = lv0 - clmcpv*(ticb(i)-273.15)
758
759 ! ori c approximation here:
760 ! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
761 ! ori & -gz(i,icb(i))-alv*qg)/cpd
762
763 ! convect3: no approximation:
764 238560 tp(i, icbs(i)) = (ah0(i)-gz(i,icbs(i))-alv*qg)/(cpd+(cl-cpd)*qnk(i))
765
766 ! ori clw(i,icb(i))=qnk(i)-qg
767 ! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
768 238560 clw(i, icbs(i)) = qnk(i) - qg
769 238560 clw(i, icbs(i)) = max(0.0, clw(i,icbs(i)))
770
771 rg = qg/(1.-qnk(i))
772 ! ori tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
773 ! convect3: (qg utilise au lieu du vrai mixing ratio rg)
774 238800 tvp(i, icbs(i)) = tp(i, icbs(i))*(1.+qg/eps-qnk(i)) !whole thing
775
776 END DO
777
778 ! ori do 380 k=minorig,icbsmax2
779 ! ori do 370 i=1,len
780 ! ori tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
781 ! ori 370 continue
782 ! ori 380 continue
783
784
785 ! -- The following is only for convect3:
786
787 ! * icbs is the first level above the LCL:
788 ! if plcl<p(icb), then icbs=icb+1
789 ! if plcl>p(icb), then icbs=icb
790
791 ! * the routine above computes tvp from minorig to icbs (included).
792
793 ! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
794 ! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
795
796 ! * therefore, in the case icbs=icb, we compute tvp at level icb+1
797 ! (tvp at other levels will be computed in cv3_undilute2.F)
798
799
800
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
801 238560 ticb(i) = t(i, icb(i)+1)
802 238560 gzicb(i) = gz(i, icb(i)+1)
803 238800 qsicb(i) = qs(i, icb(i)+1)
804 END DO
805
806
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
807 238560 tg = ticb(i)
808 238560 qg = qsicb(i) ! convect3
809 ! debug alv=lv0-clmcpv*(ticb(i)-t0)
810 238560 alv = lv0 - clmcpv*(ticb(i)-273.15)
811
812 ! First iteration.
813
814 ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
815 s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3
816 238560 +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
817 238560 s = 1./s
818 ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
819 238560 ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
820 238560 tg = tg + s*(ah0(i)-ahg)
821 ! ori tg=max(tg,35.0)
822 ! debug tc=tg-t0
823 238560 tc = tg - 273.15
824 238560 denom = 243.5 + tc
825 238560 denom = max(denom, 1.0) ! convect3
826 ! ori if(tc.ge.0.0)then
827 238560 es = 6.112*exp(17.67*tc/denom)
828 ! ori else
829 ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
830 ! ori endif
831 ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
832 238560 qg = eps*es/(p(i,icb(i)+1)-es*(1.-eps))
833
834 ! Second iteration.
835
836
837 ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
838 ! ori s=1./s
839 ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
840 238560 ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
841 238560 tg = tg + s*(ah0(i)-ahg)
842 ! ori tg=max(tg,35.0)
843 ! debug tc=tg-t0
844 238560 tc = tg - 273.15
845 238560 denom = 243.5 + tc
846 238560 denom = max(denom, 1.0) ! convect3
847 ! ori if(tc.ge.0.0)then
848 238560 es = 6.112*exp(17.67*tc/denom)
849 ! ori else
850 ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
851 ! ori end if
852 ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
853 238560 qg = eps*es/(p(i,icb(i)+1)-es*(1.-eps))
854
855 alv = lv0 - clmcpv*(ticb(i)-273.15)
856
857 ! ori c approximation here:
858 ! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
859 ! ori & -gz(i,icb(i))-alv*qg)/cpd
860
861 ! convect3: no approximation:
862 238560 tp(i, icb(i)+1) = (ah0(i)-gz(i,icb(i)+1)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
863
864 ! ori clw(i,icb(i))=qnk(i)-qg
865 ! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
866 238560 clw(i, icb(i)+1) = qnk(i) - qg
867 238560 clw(i, icb(i)+1) = max(0.0, clw(i,icb(i)+1))
868
869 rg = qg/(1.-qnk(i))
870 ! ori tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
871 ! convect3: (qg utilise au lieu du vrai mixing ratio rg)
872 238800 tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
873
874 END DO
875
876 240 RETURN
877 END SUBROUTINE cv3_undilute1
878
879 8814891 SUBROUTINE cv3_trigger(len, nd, icb, plcl, p, th, tv, tvp, thnk, &
880 240 pbase, buoybase, iflag, sig, w0)
881 IMPLICIT NONE
882
883 ! -------------------------------------------------------------------
884 ! --- TRIGGERING
885
886 ! - computes the cloud base
887 ! - triggering (crude in this version)
888 ! - relaxation of sig and w0 when no convection
889
890 ! Caution1: if no convection, we set iflag=14
891 ! (it used to be 0 in convect3)
892
893 ! Caution2: at this stage, tvp (and thus buoy) are know up
894 ! through icb only!
895 ! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
896 ! -------------------------------------------------------------------
897
898 include "cv3param.h"
899
900 ! input:
901 INTEGER len, nd
902 INTEGER icb(len)
903 REAL plcl(len), p(len, nd)
904 REAL th(len, nd), tv(len, nd), tvp(len, nd)
905 REAL thnk(len)
906
907 ! output:
908 REAL pbase(len), buoybase(len)
909
910 ! input AND output:
911 INTEGER iflag(len)
912 REAL sig(len, nd), w0(len, nd)
913
914 ! local variables:
915 INTEGER i, k
916 REAL tvpbase, tvbase, tdif, ath, ath1
917
918
919 ! *** set cloud base buoyancy at (plcl+dpbase) level buoyancy
920
921
2/2
✓ Branch 0 taken 238560 times.
✓ Branch 1 taken 240 times.
 238800 DO i = 1, len
922 238560 pbase(i) = plcl(i) + dpbase
923 tvpbase = tvp(i, icb(i)) *(pbase(i)-p(i,icb(i)+1))/(p(i,icb(i))-p(i,icb(i)+1)) + &
924 238560 tvp(i, icb(i)+1)*(p(i,icb(i))-pbase(i)) /(p(i,icb(i))-p(i,icb(i)+1))
925 tvbase = tv(i, icb(i)) *(pbase(i)-p(i,icb(i)+1))/(p(i,icb(i))-p(i,icb(i)+1)) + &
926 238560 tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i)) /(p(i,icb(i))-p(i,icb(i)+1))
927 238800 buoybase(i) = tvpbase - tvbase
928 END DO
929
930
931 ! *** make sure that column is dry adiabatic between the surface ***
932 ! *** and cloud base, and that lifted air is positively buoyant ***
933 ! *** at cloud base ***
934 ! *** if not, return to calling program after resetting ***
935 ! *** sig(i) and w0(i) ***
936
937
938 ! oct3 do 200 i=1,len
939 ! oct3
940 ! oct3 tdif = buoybase(i)
941 ! oct3 ath1 = th(i,1)
942 ! oct3 ath = th(i,icb(i)-1) - dttrig
943 ! oct3
944 ! oct3 if (tdif.lt.dtcrit .or. ath.gt.ath1) then
945 ! oct3 do 60 k=1,nl
946 ! oct3 sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif
947 ! oct3 sig(i,k) = AMAX1(sig(i,k),0.0)
948 ! oct3 w0(i,k) = beta*w0(i,k)
949 ! oct3 60 continue
950 ! oct3 iflag(i)=4 ! pour version vectorisee
951 ! oct3c convect3 iflag(i)=0
952 ! oct3cccc return
953 ! oct3 endif
954 ! oct3
955 ! oct3200 continue
956
957 ! -- oct3: on reecrit la boucle 200 (pour la vectorisation)
958
959
2/2
✓ Branch 0 taken 6480 times.
✓ Branch 1 taken 240 times.
 6720 DO k = 1, nl
960
2/2
✓ Branch 0 taken 6441120 times.
✓ Branch 1 taken 6480 times.
 6447840 DO i = 1, len
961
962 6441120 tdif = buoybase(i)
963 6441120 ath1 = thnk(i)
964 6441120 ath = th(i, icb(i)-1) - dttrig
965
966
3/4
✓ Branch 0 taken 4306149 times.
✓ Branch 1 taken 2134971 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4306149 times.
 6447600 IF (tdif<dtcrit .OR. ath>ath1) THEN
967 2134971 sig(i, k) = beta*sig(i, k) - 2.*alpha*tdif*tdif
968 2134971 sig(i, k) = amax1(sig(i,k), 0.0)
969 2134971 w0(i, k) = beta*w0(i, k)
970 2134971 iflag(i) = 14 ! pour version vectorisee
971 ! convect3 iflag(i)=0
972 END IF
973
974 END DO
975 END DO
976
977 ! fin oct3 --
978
979 240 RETURN
980 END SUBROUTINE cv3_trigger
981
982 ✗ SUBROUTINE cv3_compress(len, nloc, ncum, nd, ntra, &
983 ✗ iflag1, nk1, icb1, icbs1, &
984 plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, &
985 t1, q1, qs1, u1, v1, gz1, th1, &
986 tra1, &
987 h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, &
988 ✗ sig1, w01, &
989 iflag, nk, icb, icbs, &
990 plcl, tnk, qnk, gznk, pbase, buoybase, &
991 t, q, qs, u, v, gz, th, &
992 tra, &
993 h, lv, cpn, p, ph, tv, tp, tvp, clw, &
994 ✗ sig, w0)
995 USE print_control_mod, ONLY: lunout
996 IMPLICIT NONE
997
998 include "cv3param.h"
999
1000 !inputs:
1001 INTEGER len, ncum, nd, ntra, nloc
1002 INTEGER iflag1(len), nk1(len), icb1(len), icbs1(len)
1003 REAL plcl1(len), tnk1(len), qnk1(len), gznk1(len)
1004 REAL pbase1(len), buoybase1(len)
1005 REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
1006 REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
1007 REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
1008 REAL tvp1(len, nd), clw1(len, nd)
1009 REAL th1(len, nd)
1010 REAL sig1(len, nd), w01(len, nd)
1011 REAL tra1(len, nd, ntra)
1012
1013 !outputs:
1014 ! en fait, on a nloc=len pour l'instant (cf cv_driver)
1015 INTEGER iflag(nloc), nk(nloc), icb(nloc), icbs(nloc)
1016 REAL plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
1017 REAL pbase(nloc), buoybase(nloc)
1018 REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
1019 REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
1020 REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
1021 REAL tvp(nloc, nd), clw(nloc, nd)
1022 REAL th(nloc, nd)
1023 REAL sig(nloc, nd), w0(nloc, nd)
1024 REAL tra(nloc, nd, ntra)
1025
1026 !local variables:
1027 INTEGER i, k, nn, j
1028
1029 CHARACTER (LEN=20) :: modname = 'cv3_compress'
1030 CHARACTER (LEN=80) :: abort_message
1031
1032 ✗ DO k = 1, nl + 1
1033 ✗ nn = 0
1034 ✗ DO i = 1, len
1035 ✗ IF (iflag1(i)==0) THEN
1036 ✗ nn = nn + 1
1037 ✗ sig(nn, k) = sig1(i, k)
1038 ✗ w0(nn, k) = w01(i, k)
1039 ✗ t(nn, k) = t1(i, k)
1040 ✗ q(nn, k) = q1(i, k)
1041 ✗ qs(nn, k) = qs1(i, k)
1042 ✗ u(nn, k) = u1(i, k)
1043 ✗ v(nn, k) = v1(i, k)
1044 ✗ gz(nn, k) = gz1(i, k)
1045 ✗ h(nn, k) = h1(i, k)
1046 ✗ lv(nn, k) = lv1(i, k)
1047 ✗ cpn(nn, k) = cpn1(i, k)
1048 ✗ p(nn, k) = p1(i, k)
1049 ✗ ph(nn, k) = ph1(i, k)
1050 ✗ tv(nn, k) = tv1(i, k)
1051 ✗ tp(nn, k) = tp1(i, k)
1052 ✗ tvp(nn, k) = tvp1(i, k)
1053 ✗ clw(nn, k) = clw1(i, k)
1054 ✗ th(nn, k) = th1(i, k)
1055 END IF
1056 END DO
1057 END DO
1058
1059 !AC! do 121 j=1,ntra
1060 !AC!ccccc do 111 k=1,nl+1
1061 !AC! do 111 k=1,nd
1062 !AC! nn=0
1063 !AC! do 101 i=1,len
1064 !AC! if(iflag1(i).eq.0)then
1065 !AC! nn=nn+1
1066 !AC! tra(nn,k,j)=tra1(i,k,j)
1067 !AC! endif
1068 !AC! 101 continue
1069 !AC! 111 continue
1070 !AC! 121 continue
1071
1072 ✗ IF (nn/=ncum) THEN
1073 ✗ WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
1074 ✗ abort_message = ''
1075 ✗ CALL abort_physic(modname, abort_message, 1)
1076 END IF
1077
1078 ✗ nn = 0
1079 ✗ DO i = 1, len
1080 ✗ IF (iflag1(i)==0) THEN
1081 ✗ nn = nn + 1
1082 ✗ pbase(nn) = pbase1(i)
1083 ✗ buoybase(nn) = buoybase1(i)
1084 ✗ plcl(nn) = plcl1(i)
1085 ✗ tnk(nn) = tnk1(i)
1086 ✗ qnk(nn) = qnk1(i)
1087 ✗ gznk(nn) = gznk1(i)
1088 ✗ nk(nn) = nk1(i)
1089 ✗ icb(nn) = icb1(i)
1090 ✗ icbs(nn) = icbs1(i)
1091 ✗ iflag(nn) = iflag1(i)
1092 END IF
1093 END DO
1094
1095 ✗ RETURN
1096 END SUBROUTINE cv3_compress
1097
1098 ✗ SUBROUTINE icefrac(t, clw, qi, nl, len)
1099 IMPLICIT NONE
1100
1101
1102 !JAM--------------------------------------------------------------------
1103 ! Calcul de la quantit� d'eau sous forme de glace
1104 ! --------------------------------------------------------------------
1105 INTEGER nl, len
1106 REAL qi(len, nl)
1107 REAL t(len, nl), clw(len, nl)
1108 REAL fracg
1109 INTEGER k, i
1110
1111 ✗ DO k = 3, nl
1112 ✗ DO i = 1, len
1113 ✗ IF (t(i,k)>263.15) THEN
1114 ✗ qi(i, k) = 0.
1115 ELSE
1116 ✗ IF (t(i,k)<243.15) THEN
1117 ✗ qi(i, k) = clw(i, k)
1118 ELSE
1119 ✗ fracg = (263.15-t(i,k))/20
1120 ✗ qi(i, k) = clw(i, k)*fracg
1121 END IF
1122 END IF
1123 ! print*,t(i,k),qi(i,k),'temp,testglace'
1124 END DO
1125 END DO
1126
1127 ✗ RETURN
1128
1129 END SUBROUTINE icefrac
1130
1131 20693345 SUBROUTINE cv3_undilute2(nloc, ncum, nd, iflag, icb, icbs, nk, &
1132 240 tnk, qnk, gznk, hnk, t, q, qs, gz, &
1133 p, ph, h, tv, lv, lf, pbase, buoybase, plcl, &
1134 inb, tp, tvp, clw, hp, ep, sigp, buoy, &
1135 frac_a, frac_s, qpreca, qta)
1136 USE print_control_mod, ONLY: prt_level
1137 IMPLICIT NONE
1138
1139 ! ---------------------------------------------------------------------
1140 ! Purpose:
1141 ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
1142 ! &
1143 ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
1144 ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
1145 ! &
1146 ! FIND THE LEVEL OF NEUTRAL BUOYANCY
1147
1148 ! Main differences convect3/convect4:
1149 ! - icbs (input) is the first level above LCL (may differ from icb)
1150 ! - many minor differences in the iterations
1151 ! - condensed water not removed from tvp in convect3
1152 ! - vertical profile of buoyancy computed here (use of buoybase)
1153 ! - the determination of inb is different
1154 ! - no inb1, only inb in output
1155 ! ---------------------------------------------------------------------
1156
1157 include "cvthermo.h"
1158 include "cv3param.h"
1159 include "conema3.h"
1160 include "cvflag.h"
1161 include "YOMCST2.h"
1162
1163 !inputs:
1164 INTEGER, INTENT (IN) :: ncum, nd, nloc
1165 INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, icbs, nk
1166 REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, q, qs, gz
1167 REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
1168 REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
1169 REAL, DIMENSION (nloc), INTENT (IN) :: tnk, qnk, gznk
1170 REAL, DIMENSION (nloc), INTENT (IN) :: hnk
1171 REAL, DIMENSION (nloc, nd), INTENT (IN) :: lv, lf, tv, h
1172 REAL, DIMENSION (nloc), INTENT (IN) :: pbase, buoybase, plcl
1173
1174 !input/outputs:
1175 REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: tp, tvp, clw ! Input for k = 1, icb+1 (computed in cv3_undilute1)
1176 ! Output above
1177 INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag
1178
1179 !outputs:
1180 INTEGER, DIMENSION (nloc), INTENT (OUT) :: inb
1181 REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ep, sigp, hp
1182 REAL, DIMENSION (nloc, nd), INTENT (OUT) :: buoy
1183 REAL, DIMENSION (nloc, nd), INTENT (OUT) :: frac_a, frac_s
1184 REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qpreca
1185 REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qta
1186
1187 !local variables:
1188 INTEGER i, j, k
1189 REAL smallestreal
1190 REAL tg, qg, dqgdT, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
1191 REAL :: phinu2p
1192 REAL :: qhthreshold
1193 REAL :: als
1194 REAL :: qsat_new, snew
1195 480 REAL, DIMENSION (nloc,nd) :: qi
1196 480 REAL, DIMENSION (nloc,nd) :: ha ! moist static energy of adiabatic ascents
1197 ! taking into account precip ejection
1198 480 REAL, DIMENSION (nloc,nd) :: hla ! liquid water static energy of adiabatic ascents
1199 ! taking into account precip ejection
1200 480 REAL, DIMENSION (nloc,nd) :: qcld ! specific cloud water
1201 480 REAL, DIMENSION (nloc,nd) :: qhsat ! specific humidity at saturation
1202 REAL, DIMENSION (nloc,nd) :: dqhsatdT ! dqhsat/dT
1203 480 REAL, DIMENSION (nloc,nd) :: frac ! ice fraction function of envt temperature
1204 480 REAL, DIMENSION (nloc,nd) :: qps ! specific solid precipitation
1205 480 REAL, DIMENSION (nloc,nd) :: qpl ! specific liquid precipitation
1206 480 REAL, DIMENSION (nloc) :: ah0, cape, capem, byp
1207 LOGICAL, DIMENSION (nloc) :: lcape
1208 480 INTEGER, DIMENSION (nloc) :: iposit
1209 REAL :: denomm1
1210 REAL :: by, defrac, pden, tbis
1211 REAL :: fracg
1212 REAL :: deltap
1213 REAL, SAVE :: Tx, Tm
1214 DATA Tx/263.15/, Tm/243.15/
1215 !$OMP THREADPRIVATE(Tx, Tm)
1216 REAL :: aa, bb, dd, ddelta, discr
1217 REAL :: ff, fp
1218 REAL :: coefx, coefm, Zx, Zm, Ux, U, Um
1219
1220
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 240 times.
 240 IF (prt_level >= 10) THEN
1221 ✗ print *,'cv3_undilute2.0. icvflag_Tpa, t(1,k), q(1,k), qs(1,k) ', &
1222 ✗ icvflag_Tpa, (k, t(1,k), q(1,k), qs(1,k), k = 1,nl)
1223 ENDIF
1224 smallestreal=tiny(smallestreal)
1225
1226 ! =====================================================================
1227 ! --- SOME INITIALIZATIONS
1228 ! =====================================================================
1229
1230
2/2
✓ Branch 0 taken 6480 times.
✓ Branch 1 taken 240 times.
 6720 DO k = 1, nl
1231
2/2
✓ Branch 0 taken 3431619 times.
✓ Branch 1 taken 6480 times.
 3438339 DO i = 1, ncum
1232 3438099 qi(i, k) = 0.
1233 END DO
1234 END DO
1235
1236
1237 ! =====================================================================
1238 ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
1239 ! =====================================================================
1240
1241 ! --- The procedure is to solve the equation.
1242 ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
1243
1244 ! *** Calculate certain parcel quantities, including static energy ***
1245
1246
1247
2/2
✓ Branch 0 taken 127097 times.
✓ Branch 1 taken 240 times.
 127337 DO i = 1, ncum
1248 ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i)+ &
1249 ! debug qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i)
1250 127337 qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
1251 END DO
1252 !
1253 ! Ice fraction
1254 !
1255
1/2
✓ Branch 0 taken 240 times.
✗ Branch 1 not taken.
 240 IF (cvflag_ice) THEN
1256
2/2
✓ Branch 0 taken 6480 times.
✓ Branch 1 taken 240 times.
 6720 DO k = minorig, nl
1257
2/2
✓ Branch 0 taken 3431619 times.
✓ Branch 1 taken 6480 times.
 3438339 DO i = 1, ncum
1258 3431619 frac(i, k) = (Tx - t(i,k))/(Tx - Tm)
1259 3438099 frac(i, k) = min(max(frac(i,k),0.0), 1.0)
1260 END DO
1261 END DO
1262 ! Below cloud base, set ice fraction to cloud base value
1263
2/2
✓ Branch 0 taken 6480 times.
✓ Branch 1 taken 240 times.
 6720 DO k = 1, nl
1264
2/2
✓ Branch 0 taken 3431619 times.
✓ Branch 1 taken 6480 times.
 3438339 DO i = 1, ncum
1265
2/2
✓ Branch 0 taken 520451 times.
✓ Branch 1 taken 2911168 times.
 3438099 IF (k<icb(i)) THEN
1266 520451 frac(i,k) = frac(i,icb(i))
1267 END IF
1268 END DO
1269 END DO
1270 ELSE
1271 ✗ DO k = 1, nl
1272 ✗ DO i = 1, ncum
1273 ✗ frac(i,k) = 0.
1274 END DO
1275 END DO
1276 ENDIF ! (cvflag_ice)
1277
1278
1279
2/2
✓ Branch 0 taken 6480 times.
✓ Branch 1 taken 240 times.
 6720 DO k = minorig, nl
1280
2/2
✓ Branch 0 taken 3431619 times.
✓ Branch 1 taken 6480 times.
 3438339 DO i = 1,ncum
1281 3431619 ha(i,k) = ah0(i)
1282 3431619 hla(i,k) = hnk(i)
1283 3431619 qta(i,k) = qnk(i)
1284 3431619 qpreca(i,k) = 0.
1285 3431619 frac_a(i,k) = 0.
1286 3431619 frac_s(i,k) = frac(i,k)
1287 3431619 qpl(i,k) = 0.
1288 3431619 qps(i,k) = 0.
1289 3431619 qhsat(i,k) = qs(i,k)
1290 3431619 qcld(i,k) = max(qta(i,k)-qhsat(i,k),0.)
1291
2/2
✓ Branch 0 taken 774645 times.
✓ Branch 1 taken 2656974 times.
 3438099 IF (k <= icb(i)+1) THEN
1292 774645 qhsat(i,k) = qnk(i)-clw(i,k)
1293 774645 qcld(i,k) = clw(i,k)
1294 ENDIF
1295 ENDDO
1296 ENDDO
1297
1298 !jyg<
1299 ! =====================================================================
1300 ! --- SET THE THE FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
1301 ! =====================================================================
1302
2/2
✓ Branch 0 taken 6480 times.
✓ Branch 1 taken 240 times.
 6720 DO k = 1, nl
1303
2/2
✓ Branch 0 taken 3431619 times.
✓ Branch 1 taken 6480 times.
 3438339 DO i = 1, ncum
1304 3431619 ep(i, k) = 0.0
1305 3438099 sigp(i, k) = spfac
1306 END DO
1307 END DO
1308 !>jyg
1309 !
1310
1311 ! *** Find lifted parcel quantities above cloud base ***
1312
1313 !----------------------------------------------------------------------------
1314 !
1315
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 240 times.
 240 IF (icvflag_Tpa == 2) THEN
1316 !
1317 !----------------------------------------------------------------------------
1318 !
1319 ✗ DO k = minorig + 1, nl
1320 ✗ DO i = 1,ncum
1321 ✗ tp(i,k) = t(i,k)
1322 ENDDO
1323 !! alv = lv0 - clmcpv*(t(i,k)-273.15)
1324 !! alf = lf0 + clmci*(t(i,k)-273.15)
1325 !! als = alf + alv
1326 ✗ DO j = 1,4
1327 ✗ DO i = 1, ncum
1328 ! ori if(k.ge.(icb(i)+1))then
1329 ✗ IF (k>=(icbs(i)+1)) THEN ! convect3
1330 ✗ tg = tp(i, k)
1331 ✗ IF (tg .gt. Tx) THEN
1332 ✗ es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
1333 ✗ qg = eps*es/(p(i,k)-es*(1.-eps))
1334 ELSE
1335 ✗ esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
1336 ✗ qg = eps*esi/(p(i,k)-esi*(1.-eps))
1337 ENDIF
1338 ! Ice fraction
1339 ff = 0.
1340 ✗ fp = 1./(Tx - Tm)
1341 ✗ IF (tg < Tx) THEN
1342 ✗ IF (tg > Tm) THEN
1343 ✗ ff = (Tx - tg)*fp
1344 ELSE
1345 ff = 1.
1346 ENDIF ! (tg > Tm)
1347 ENDIF ! (tg < Tx)
1348 ! Intermediate variables
1349 ✗ aa = cpd + (cl-cpd)*qnk(i) + lv(i,k)*lv(i,k)*qg/(rrv*tg*tg)
1350 ahg = (cpd + (cl-cpd)*qnk(i))*tg + lv(i,k)*qg - &
1351 ✗ lf(i,k)*ff*(qnk(i) - qg) + gz(i,k)
1352 ✗ dd = lf(i,k)*lv(i,k)*qg/(rrv*tg*tg)
1353 ✗ ddelta = lf(i,k)*(qnk(i) - qg)
1354 ✗ bb = aa + ddelta*fp + dd*fp*(Tx-tg)
1355 ! Compute Zx and Zm
1356 coefx = aa
1357 ✗ coefm = aa + dd
1358 ✗ IF (tg .gt. Tx) THEN
1359 ✗ Zx = ahg + coefx*(Tx - tg)
1360 ✗ Zm = ahg - ddelta + coefm*(Tm - tg)
1361 ELSE
1362 ✗ IF (tg .gt. Tm) THEN
1363 ✗ Zx = ahg + (coefx +fp*ddelta)*(Tx - Tg)
1364 ✗ Zm = ahg + (coefm +fp*ddelta)*(Tm - Tg)
1365 ELSE
1366 ✗ Zx = ahg + ddelta + coefx*(Tx - tg)
1367 ✗ Zm = ahg + coefm*(Tm - tg)
1368 ENDIF ! (tg .gt. Tm)
1369 ENDIF ! (tg .gt. Tx)
1370 ! Compute the masks Um, U, Ux
1371 ✗ Um = (sign(1., Zm-ah0(i))+1.)/2.
1372 ✗ Ux = (sign(1., ah0(i)-Zx)+1.)/2.
1373 ✗ U = (1. - Um)*(1. - Ux)
1374 ! Compute the updated parcell temperature Tp : 3 cases depending on tg value
1375 ✗ IF (tg .gt. Tx) THEN
1376 ✗ discr = bb*bb - 4*dd*fp*(ah0(i) - ahg + ddelta*fp*(Tx-tg))
1377 Tp(i,k) = tg + &
1378 Um* (ah0(i) - ahg + ddelta) /(aa + dd) + &
1379 U *2*(ah0(i) - ahg + ddelta*fp*(Tx-tg))/(bb + sqrt(discr)) + &
1380 ✗ Ux* (ah0(i) - ahg) /aa
1381 ✗ ELSEIF (tg .gt. Tm) THEN
1382 ✗ discr = bb*bb - 4*dd*fp*(ah0(i) - ahg)
1383 Tp(i,k) = tg + &
1384 Um* (ah0(i) - ahg + ddelta*fp*(tg-Tm))/(aa + dd) + &
1385 U *2*(ah0(i) - ahg) /(bb + sqrt(discr)) + &
1386 ✗ Ux* (ah0(i) - ahg + ddelta*fp*(tg-Tx))/aa
1387 ELSE
1388 ✗ discr = bb*bb - 4*dd*fp*(ah0(i) - ahg + ddelta*fp*(Tm-tg))
1389 Tp(i,k) = tg + &
1390 Um* (ah0(i) - ahg) /(aa + dd) + &
1391 U *2*(ah0(i) - ahg + ddelta*fp*(Tm-tg))/(bb + sqrt(discr)) + &
1392 ✗ Ux* (ah0(i) - ahg - ddelta) /aa
1393 ENDIF ! (tg .gt. Tx)
1394 !
1395 !! print *,' j, k, Um, U, Ux, aa, bb, discr, dd, ddelta ', j, k, Um, U, Ux, aa, bb, discr, dd, ddelta
1396 !! print *,' j, k, ah0(i), ahg, tg, qg, tp(i,k), ff ', j, k, ah0(i), ahg, tg, qg, tp(i,k), ff
1397 END IF ! (k>=(icbs(i)+1))
1398 END DO ! i = 1, ncum
1399 END DO ! j = 1,4
1400 ✗ DO i = 1, ncum
1401 ✗ IF (k>=(icbs(i)+1)) THEN ! convect3
1402 ✗ tg = tp(i, k)
1403 ✗ IF (tg .gt. Tx) THEN
1404 ✗ es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
1405 ✗ qg = eps*es/(p(i,k)-es*(1.-eps))
1406 ELSE
1407 ✗ esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
1408 ✗ qg = eps*esi/(p(i,k)-esi*(1.-eps))
1409 ENDIF
1410 ✗ clw(i, k) = qnk(i) - qg
1411 ✗ clw(i, k) = max(0.0, clw(i,k))
1412 ✗ tvp(i, k) = max(0., tp(i,k)*(1.+qg/eps-qnk(i)))
1413 ! print*,tvp(i,k),'tvp'
1414 ✗ IF (clw(i,k)<1.E-11) THEN
1415 ✗ tp(i, k) = tv(i, k)
1416 ✗ tvp(i, k) = tv(i, k)
1417 END IF ! (clw(i,k)<1.E-11)
1418 END IF ! (k>=(icbs(i)+1))
1419 END DO ! i = 1, ncum
1420 END DO ! k = minorig + 1, nl
1421 !----------------------------------------------------------------------------
1422 !
1423
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 240 times.
 240 ELSE IF (icvflag_Tpa == 1) THEN ! (icvflag_Tpa == 2)
1424 !
1425 !----------------------------------------------------------------------------
1426 !
1427 ✗ DO k = minorig + 1, nl
1428 ✗ DO i = 1,ncum
1429 ✗ tp(i,k) = t(i,k)
1430 ENDDO
1431 !! alv = lv0 - clmcpv*(t(i,k)-273.15)
1432 !! alf = lf0 + clmci*(t(i,k)-273.15)
1433 !! als = alf + alv
1434 ✗ DO j = 1,4
1435 ✗ DO i = 1, ncum
1436 ! ori if(k.ge.(icb(i)+1))then
1437 ✗ IF (k>=(icbs(i)+1)) THEN ! convect3
1438 ✗ tg = tp(i, k)
1439 ✗ IF (tg .gt. Tx .OR. .NOT.cvflag_ice) THEN
1440 ✗ es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
1441 ✗ qg = eps*es/(p(i,k)-es*(1.-eps))
1442 ✗ dqgdT = lv(i,k)*qg/(rrv*tg*tg)
1443 ELSE
1444 ✗ esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
1445 ✗ qg = eps*esi/(p(i,k)-esi*(1.-eps))
1446 ✗ dqgdT = (lv(i,k)+lf(i,k))*qg/(rrv*tg*tg)
1447 ENDIF
1448 ✗ IF (qsat_depends_on_qt) THEN
1449 ✗ dqgdT = dqgdT*(1.-qta(i,k-1))/(1.-qg)**2
1450 ✗ qg = qg*(1.-qta(i,k-1))/(1.-qg)
1451 ENDIF
1452 ahg = (cpd + (cl-cpd)*qta(i,k-1))*tg + lv(i,k)*qg - &
1453 ✗ lf(i,k)*frac(i,k)*(qta(i,k-1) - qg) + gz(i,k)
1454 Tp(i,k) = tg + (ah0(i) - ahg)/ &
1455 ✗ (cpd + (cl-cpd)*qta(i,k-1) + (lv(i,k)+frac(i,k)*lf(i,k))*dqgdT)
1456 !! print *,'undilute2 iterations k, Tp(i,k), ah0(i), ahg ', &
1457 !! k, Tp(i,k), ah0(i), ahg
1458 END IF ! (k>=(icbs(i)+1))
1459 END DO ! i = 1, ncum
1460 END DO ! j = 1,4
1461 ✗ DO i = 1, ncum
1462 ✗ IF (k>=(icbs(i)+1)) THEN ! convect3
1463 ✗ tg = tp(i, k)
1464 ✗ IF (tg .gt. Tx .OR. .NOT.cvflag_ice) THEN
1465 ✗ es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
1466 ✗ qg = eps*es/(p(i,k)-es*(1.-eps))
1467 ELSE
1468 ✗ esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
1469 ✗ qg = eps*esi/(p(i,k)-esi*(1.-eps))
1470 ENDIF
1471 ✗ IF (qsat_depends_on_qt) THEN
1472 ✗ qg = qg*(1.-qta(i,k-1))/(1.-qg)
1473 ENDIF
1474 ✗ qhsat(i,k) = qg
1475 END IF ! (k>=(icbs(i)+1))
1476 END DO ! i = 1, ncum
1477 ✗ DO i = 1, ncum
1478 ✗ IF (k>=(icbs(i)+1)) THEN ! convect3
1479 ✗ clw(i, k) = qta(i,k-1) - qhsat(i,k)
1480 ✗ clw(i, k) = max(0.0, clw(i,k))
1481 ✗ tvp(i, k) = max(0., tp(i,k)*(1.+qhsat(i,k)/eps-qta(i,k-1)))
1482 ! print*,tvp(i,k),'tvp'
1483 ✗ IF (clw(i,k)<1.E-11) THEN
1484 ✗ tp(i, k) = tv(i, k)
1485 ✗ tvp(i, k) = tv(i, k)
1486 END IF ! (clw(i,k)<1.E-11)
1487 END IF ! (k>=(icbs(i)+1))
1488 END DO ! i = 1, ncum
1489 !
1490 ✗ IF (cvflag_prec_eject) THEN
1491 ✗ DO i = 1, ncum
1492 ✗ IF (k>=(icbs(i)+1)) THEN ! convect3
1493 ! Specific precipitation (liquid and solid) and ice content
1494 ! before ejection of precipitation !!jygprl
1495 ✗ elacrit = elcrit*min(max(1.-(tp(i,k)-T0)/Tlcrit, 0.), 1.) !!jygprl
1496 !!!! qcld(i,k) = min(clw(i,k), elacrit) !!jygprl
1497 ✗ qhthreshold = elacrit*(1.-qta(i,k-1))/(1.-elacrit)
1498 ✗ qcld(i,k) = min(clw(i,k), qhthreshold) !!jygprl
1499 !!!! phinu2p = max(qhsat(i,k-1) + qcld(i,k-1) - (qhsat(i,k) + qcld(i,k)),0.) !!jygprl
1500 ✗ phinu2p = max(clw(i,k) - max(qta(i,k-1) - qhsat(i,k-1), qhthreshold), 0.)
1501 ✗ qpl(i,k) = qpl(i,k-1) + (1.-frac(i,k))*phinu2p !!jygprl
1502 ✗ qps(i,k) = qps(i,k-1) + frac(i,k) *phinu2p !!jygprl
1503 qi(i,k) = (1.-ejectliq)*clw(i,k)*frac(i,k) + & !!jygprl
1504 ✗ ejectliq*(qps(i,k-1) + frac(i,k)*(phinu2p+qcld(i,k))) !!jygprl
1505 !!
1506 ! =====================================================================================
1507 ! Ejection of precipitation from adiabatic ascents if requested (cvflag_prec_eject=True):
1508 ! Compute the steps of total water (qta), of moist static energy (ha), of specific
1509 ! precipitation (qpl and qps) and of specific cloud water (qcld) associated with precipitation
1510 ! ejection.
1511 ! =====================================================================================
1512 !
1513 ! Verif
1514 ✗ qpreca(i,k) = ejectliq*qpl(i,k) + ejectice*qps(i,k) !!jygprl
1515 ✗ frac_a(i,k) = ejectice*qps(i,k)/max(qpreca(i,k),smallestreal) !!jygprl
1516 frac_s(i,k) = (1.-ejectliq)*frac(i,k) + & !!jygprl
1517 ✗ ejectliq*(1. - (qpl(i,k)+(1.-frac(i,k))*qcld(i,k))/max(clw(i,k),smallestreal)) !!jygprl
1518 !
1519 ✗ denomm1 = 1./(1. - qpreca(i,k))
1520 !
1521 qta(i,k) = qta(i,k-1) - &
1522 ✗ qpreca(i,k)*(1.-qta(i,k-1))*denomm1
1523 ha(i,k) = ha(i,k-1) + &
1524 ( qpreca(i,k)*(-(1.-qta(i,k-1))*(cl-cpd)*tp(i,k) + &
1525 lv(i,k)*qhsat(i,k) - lf(i,k)*(frac_s(i,k)*qcld(i,k)+qps(i,k))) + &
1526 ✗ lf(i,k)*ejectice*qps(i,k))*denomm1
1527 hla(i,k) = hla(i,k-1) + &
1528 ( qpreca(i,k)*(-(1.-qta(i,k-1))*(cpv-cpd)*tp(i,k) - &
1529 lv(i,k)*((1.-frac_s(i,k))*qcld(i,k)+qpl(i,k)) - &
1530 (lv(i,k)+lf(i,k))*(frac_s(i,k)*qcld(i,k)+qps(i,k))) + &
1531 ✗ lv(i,k)*ejectliq*qpl(i,k) + (lv(i,k)+lf(i,k))*ejectice*qps(i,k))*denomm1
1532 ✗ qpl(i,k) = qpl(i,k)*(1.-ejectliq)*denomm1
1533 ✗ qps(i,k) = qps(i,k)*(1.-ejectice)*denomm1
1534 ✗ qcld(i,k) = qcld(i,k)*denomm1
1535 ✗ qhsat(i,k) = qhsat(i,k)*(1.-qta(i,k))/(1.-qta(i,k-1))
1536 END IF ! (k>=(icbs(i)+1))
1537 END DO ! i = 1, ncum
1538 ENDIF ! (cvflag_prec_eject)
1539 !
1540 END DO ! k = minorig + 1, nl
1541 !
1542 !----------------------------------------------------------------------------
1543 !
1544
1/2
✓ Branch 0 taken 240 times.
✗ Branch 1 not taken.
 240 ELSE IF (icvflag_Tpa == 0) THEN! (icvflag_Tpa == 2) ELSE IF(icvflag_Tpa == 1)
1545 !
1546 !----------------------------------------------------------------------------
1547 !
1548
2/2
✓ Branch 0 taken 240 times.
✓ Branch 1 taken 6240 times.
 6480 DO k = minorig + 1, nl
1549
2/2
✓ Branch 0 taken 3304522 times.
✓ Branch 1 taken 6240 times.
 3311002 DO i = 1, ncum
1550 ! ori if(k.ge.(icb(i)+1))then
1551
2/2
✓ Branch 0 taken 2715198 times.
✓ Branch 1 taken 589324 times.
 3310762 IF (k>=(icbs(i)+1)) THEN ! convect3
1552 2715198 tg = t(i, k)
1553 2715198 qg = qs(i, k)
1554 ! debug alv=lv0-clmcpv*(t(i,k)-t0)
1555 2715198 alv = lv0 - clmcpv*(t(i,k)-273.15)
1556
1557 ! First iteration.
1558
1559 ! ori s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
1560 s = cpd*(1.-qnk(i)) + cl*qnk(i) + & ! convect3
1561 2715198 alv*alv*qg/(rrv*t(i,k)*t(i,k)) ! convect3
1562 2715198 s = 1./s
1563 ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
1564 2715198 ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3
1565 2715198 tg = tg + s*(ah0(i)-ahg)
1566 ! ori tg=max(tg,35.0)
1567 ! debug tc=tg-t0
1568 2715198 tc = tg - 273.15
1569 2715198 denom = 243.5 + tc
1570 2715198 denom = max(denom, 1.0) ! convect3
1571 ! ori if(tc.ge.0.0)then
1572 2715198 es = 6.112*exp(17.67*tc/denom)
1573 ! ori else
1574 ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
1575 ! ori endif
1576 2715198 qg = eps*es/(p(i,k)-es*(1.-eps))
1577
1578 ! Second iteration.
1579
1580 ! ori s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
1581 ! ori s=1./s
1582 ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
1583 2715198 ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3
1584 2715198 tg = tg + s*(ah0(i)-ahg)
1585 ! ori tg=max(tg,35.0)
1586 ! debug tc=tg-t0
1587 2715198 tc = tg - 273.15
1588 2715198 denom = 243.5 + tc
1589 2715198 denom = max(denom, 1.0) ! convect3
1590 ! ori if(tc.ge.0.0)then
1591 2715198 es = 6.112*exp(17.67*tc/denom)
1592 ! ori else
1593 ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
1594 ! ori endif
1595 2715198 qg = eps*es/(p(i,k)-es*(1.-eps))
1596
1597 ! debug alv=lv0-clmcpv*(t(i,k)-t0)
1598 alv = lv0 - clmcpv*(t(i,k)-273.15)
1599 ! print*,'cpd dans convect2 ',cpd
1600 ! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
1601 ! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
1602
1603 ! ori c approximation here:
1604 ! ori tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
1605
1606 ! convect3: no approximation:
1607
1/2
✓ Branch 0 taken 2715198 times.
✗ Branch 1 not taken.
 2715198 IF (cvflag_ice) THEN
1608 2715198 tp(i, k) = max(0., (ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i)))
1609 ELSE
1610 ✗ tp(i, k) = (ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
1611 END IF
1612
1613 2715198 clw(i, k) = qnk(i) - qg
1614 2715198 clw(i, k) = max(0.0, clw(i,k))
1615 rg = qg/(1.-qnk(i))
1616 ! ori tvp(i,k)=tp(i,k)*(1.+rg*epsi)
1617 ! convect3: (qg utilise au lieu du vrai mixing ratio rg):
1618 2715198 tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
1619
1/2
✓ Branch 0 taken 2715198 times.
✗ Branch 1 not taken.
 2715198 IF (cvflag_ice) THEN
1620
2/2
✓ Branch 0 taken 748 times.
✓ Branch 1 taken 2714450 times.
 2715198 IF (clw(i,k)<1.E-11) THEN
1621 748 tp(i, k) = tv(i, k)
1622 748 tvp(i, k) = tv(i, k)
1623 END IF
1624 END IF
1625 !jyg<
1626 !! END IF ! Endif moved to the end of the loop
1627 !>jyg
1628
1629
1/2
✓ Branch 0 taken 2715198 times.
✗ Branch 1 not taken.
 2715198 IF (cvflag_ice) THEN
1630 !CR:attention boucle en klon dans Icefrac
1631 ! Call Icefrac(t,clw,qi,nl,nloc)
1632
2/2
✓ Branch 0 taken 599190 times.
✓ Branch 1 taken 2116008 times.
 2715198 IF (t(i,k)>263.15) THEN
1633 599190 qi(i, k) = 0.
1634 ELSE
1635
2/2
✓ Branch 0 taken 1778803 times.
✓ Branch 1 taken 337205 times.
 2116008 IF (t(i,k)<243.15) THEN
1636 1778803 qi(i, k) = clw(i, k)
1637 ELSE
1638 337205 fracg = (263.15-t(i,k))/20
1639 337205 qi(i, k) = clw(i, k)*fracg
1640 END IF
1641 END IF
1642 !CR: fin test
1643
2/2
✓ Branch 0 taken 2116008 times.
✓ Branch 1 taken 599190 times.
 2715198 IF (t(i,k)<263.15) THEN
1644 !CR: on commente les calculs d'Arnaud car division par zero
1645 ! nouveau calcul propose par JYG
1646 ! alv=lv0-clmcpv*(t(i,k)-273.15)
1647 ! alf=lf0-clmci*(t(i,k)-273.15)
1648 ! tg=tp(i,k)
1649 ! tc=tp(i,k)-273.15
1650 ! denom=243.5+tc
1651 ! do j=1,3
1652 ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1653 ! il faudra que esi vienne en argument de la convection
1654 ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1655 ! tbis=t(i,k)+(tp(i,k)-tg)
1656 ! esi=exp(23.33086-(6111.72784/tbis) + &
1657 ! 0.15215*log(tbis))
1658 ! qsat_new=eps*esi/(p(i,k)-esi*(1.-eps))
1659 ! snew=cpd*(1.-qnk(i))+cl*qnk(i)+alv*alv*qsat_new/ &
1660 ! (rrv*tbis*tbis)
1661 ! snew=1./snew
1662 ! print*,esi,qsat_new,snew,'esi,qsat,snew'
1663 ! tp(i,k)=tg+(alf*qi(i,k)+alv*qg*(1.-(esi/es)))*snew
1664 ! print*,k,tp(i,k),qnk(i),'avec glace'
1665 ! print*,'tpNAN',tg,alf,qi(i,k),alv,qg,esi,es,snew
1666 ! enddo
1667
1668 alv = lv0 - clmcpv*(t(i,k)-273.15)
1669 2116008 alf = lf0 + clmci*(t(i,k)-273.15)
1670 2116008 als = alf + alv
1671 2116008 tg = tp(i, k)
1672 2116008 tp(i, k) = t(i, k)
1673
2/2
✓ Branch 0 taken 6348024 times.
✓ Branch 1 taken 2116008 times.
 8464032 DO j = 1, 3
1674 6348024 esi = exp(23.33086-(6111.72784/tp(i,k))+0.15215*log(tp(i,k)))
1675 6348024 qsat_new = eps*esi/(p(i,k)-esi*(1.-eps))
1676 snew = cpd*(1.-qnk(i)) + cl*qnk(i) + alv*als*qsat_new/ &
1677 6348024 (rrv*tp(i,k)*tp(i,k))
1678 6348024 snew = 1./snew
1679 ! c print*,esi,qsat_new,snew,'esi,qsat,snew'
1680 tp(i, k) = tp(i, k) + &
1681 ((cpd*(1.-qnk(i))+cl*qnk(i))*(tg-tp(i,k)) + &
1682 8464032 alv*(qg-qsat_new)+alf*qi(i,k))*snew
1683 ! print*,k,tp(i,k),qsat_new,qnk(i),qi(i,k), &
1684 ! 'k,tp,q,qt,qi avec glace'
1685 END DO
1686
1687 !CR:reprise du code AJ
1688 2116008 clw(i, k) = qnk(i) - qsat_new
1689 2116008 clw(i, k) = max(0.0, clw(i,k))
1690 2116008 tvp(i, k) = max(0., tp(i,k)*(1.+qsat_new/eps-qnk(i)))
1691 ! print*,tvp(i,k),'tvp'
1692 END IF
1693
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 2715189 times.
 2715198 IF (clw(i,k)<1.E-11) THEN
1694 9 tp(i, k) = tv(i, k)
1695 9 tvp(i, k) = tv(i, k)
1696 END IF
1697 END IF ! (cvflag_ice)
1698 !jyg<
1699 END IF ! (k>=(icbs(i)+1))
1700 !>jyg
1701 END DO
1702 END DO
1703
1704 !----------------------------------------------------------------------------
1705 !
1706 ENDIF ! (icvflag_Tpa == 2) ELSEIF (icvflag_Tpa == 1) ELSE (icvflag_Tpa == 0)
1707 !
1708 !----------------------------------------------------------------------------
1709 !
1710 ! =====================================================================
1711 ! --- SET THE PRECIPITATION EFFICIENCIES
1712 ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
1713 ! =====================================================================
1714 !
1715
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 240 times.
 240 IF (flag_epkeorig/=1) THEN
1716 ✗ DO k = 1, nl ! convect3
1717 ✗ DO i = 1, ncum
1718 !jyg<
1719 ✗ IF(k>=icb(i)) THEN
1720 !>jyg
1721 ✗ pden = ptcrit - pbcrit
1722 ✗ ep(i, k) = (plcl(i)-p(i,k)-pbcrit)/pden*epmax
1723 ✗ ep(i, k) = max(ep(i,k), 0.0)
1724 ✗ ep(i, k) = min(ep(i,k), epmax)
1725 !! sigp(i, k) = spfac ! jyg
1726 ENDIF ! (k>=icb(i))
1727 END DO
1728 END DO
1729 ELSE
1730
2/2
✓ Branch 0 taken 240 times.
✓ Branch 1 taken 6480 times.
 6720 DO k = 1, nl
1731
2/2
✓ Branch 0 taken 3431619 times.
✓ Branch 1 taken 6480 times.
 3438339 DO i = 1, ncum
1732
2/2
✓ Branch 0 taken 2911168 times.
✓ Branch 1 taken 520451 times.
 3438099 IF(k>=icb(i)) THEN
1733 !! IF (k>=(nk(i)+1)) THEN
1734 !>jyg
1735 2911168 tca = tp(i, k) - t0
1736
2/2
✓ Branch 0 taken 563593 times.
✓ Branch 1 taken 2347575 times.
 2911168 IF (tca>=0.0) THEN
1737 563593 elacrit = elcrit
1738 ELSE
1739 2347575 elacrit = elcrit*(1.0-tca/tlcrit)
1740 END IF
1741 2911168 elacrit = max(elacrit, 0.0)
1742 2911168 ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8)
1743 2911168 ep(i, k) = max(ep(i,k), 0.0)
1744 2911168 ep(i, k) = min(ep(i,k), epmax)
1745 !! sigp(i, k) = spfac ! jyg
1746 END IF ! (k>=icb(i))
1747 END DO
1748 END DO
1749 END IF
1750 !
1751 ! =========================================================================
1752
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 240 times.
 240 IF (prt_level >= 10) THEN
1753 ✗ print *,'cv3_undilute2.1. tp(1,k), tvp(1,k) ', &
1754 ✗ (k, tp(1,k), tvp(1,k), k = 1,nl)
1755 ENDIF
1756 !
1757 ! =====================================================================
1758 ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
1759 ! --- VIRTUAL TEMPERATURE
1760 ! =====================================================================
1761
1762 ! dans convect3, tvp est calcule en une seule fois, et sans retirer
1763 ! l'eau condensee (~> reversible CAPE)
1764
1765 ! ori do 340 k=minorig+1,nl
1766 ! ori do 330 i=1,ncum
1767 ! ori if(k.ge.(icb(i)+1))then
1768 ! ori tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
1769 ! oric print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
1770 ! oric print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
1771 ! ori endif
1772 ! ori 330 continue
1773 ! ori 340 continue
1774
1775 ! ori do 350 i=1,ncum
1776 ! ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
1777 ! ori 350 continue
1778
1779
2/2
✓ Branch 0 taken 127097 times.
✓ Branch 1 taken 240 times.
 127337 DO i = 1, ncum ! convect3
1780 127337 tp(i, nlp) = tp(i, nl) ! convect3
1781 END DO ! convect3
1782
1783 ! =====================================================================
1784 ! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
1785 ! =====================================================================
1786
1787 ! -- this is for convect3 only:
1788
1789 ! first estimate of buoyancy:
1790
1791 !jyg : k-loop outside i-loop (07042015)
1792
2/2
✓ Branch 0 taken 6480 times.
✓ Branch 1 taken 240 times.
 6720 DO k = 1, nl
1793
2/2
✓ Branch 0 taken 3431619 times.
✓ Branch 1 taken 6480 times.
 3438339 DO i = 1, ncum
1794 3438099 buoy(i, k) = tvp(i, k) - tv(i, k)
1795 END DO
1796 END DO
1797
1798 ! set buoyancy=buoybase for all levels below base
1799 ! for safety, set buoy(icb)=buoybase
1800
1801 !jyg : k-loop outside i-loop (07042015)
1802
2/2
✓ Branch 0 taken 6480 times.
✓ Branch 1 taken 240 times.
 6720 DO k = 1, nl
1803
2/2
✓ Branch 0 taken 3431619 times.
✓ Branch 1 taken 6480 times.
 3438339 DO i = 1, ncum
1804
4/4
✓ Branch 0 taken 2911168 times.
✓ Branch 1 taken 520451 times.
✓ Branch 2 taken 229384 times.
✓ Branch 3 taken 2681784 times.
 3438099 IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i,k)>=pbase(i))) THEN
1805 229384 buoy(i, k) = buoybase(i)
1806 END IF
1807 END DO
1808 END DO
1809
2/2
✓ Branch 0 taken 127097 times.
✓ Branch 1 taken 240 times.
 127337 DO i = 1, ncum
1810 ! buoy(icb(i),k)=buoybase(i)
1811 127337 buoy(i, icb(i)) = buoybase(i)
1812 END DO
1813
1814 ! -- end convect3
1815
1816 ! =====================================================================
1817 ! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
1818 ! --- LEVEL OF NEUTRAL BUOYANCY
1819 ! =====================================================================
1820
1821 ! -- this is for convect3 only:
1822
1823
2/2
✓ Branch 0 taken 127097 times.
✓ Branch 1 taken 240 times.
 127337 DO i = 1, ncum
1824 127097 inb(i) = nl - 1
1825 127337 iposit(i) = nl
1826 END DO
1827
1828
1829 ! -- iposit(i) = first level, above icb, with positive buoyancy
1830
2/2
✓ Branch 0 taken 6240 times.
✓ Branch 1 taken 240 times.
 6480 DO k = 1, nl - 1
1831
2/2
✓ Branch 0 taken 3304522 times.
✓ Branch 1 taken 6240 times.
 3311002 DO i = 1, ncum
1832
4/4
✓ Branch 0 taken 2784071 times.
✓ Branch 1 taken 520451 times.
✓ Branch 2 taken 992890 times.
✓ Branch 3 taken 1791181 times.
 3310762 IF (k>=icb(i) .AND. buoy(i,k)>0.) THEN
1833 992890 iposit(i) = min(iposit(i), k)
1834 END IF
1835 END DO
1836 END DO
1837
1838
2/2
✓ Branch 0 taken 127097 times.
✓ Branch 1 taken 240 times.
 127337 DO i = 1, ncum
1839
2/2
✓ Branch 0 taken 8529 times.
✓ Branch 1 taken 118568 times.
 127337 IF (iposit(i)==nl) THEN
1840 8529 iposit(i) = icb(i)
1841 END IF
1842 END DO
1843
1844
2/2
✓ Branch 0 taken 6240 times.
✓ Branch 1 taken 240 times.
 6480 DO k = 1, nl - 1
1845
2/2
✓ Branch 0 taken 3304522 times.
✓ Branch 1 taken 6240 times.
 3311002 DO i = 1, ncum
1846
4/4
✓ Branch 0 taken 2732454 times.
✓ Branch 1 taken 572068 times.
✓ Branch 2 taken 1715764 times.
✓ Branch 3 taken 1016690 times.
 3310762 IF ((k>=iposit(i)) .AND. (buoy(i,k)<dtovsh)) THEN
1847 1715764 inb(i) = min(inb(i), k)
1848 END IF
1849 END DO
1850 END DO
1851
1852 !CR fix computation of inb
1853 !keep flag or modify in all cases?
1854
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 240 times.
 240 IF (iflag_mix_adiab.eq.1) THEN
1855 ✗ DO i = 1, ncum
1856 ✗ cape(i)=0.
1857 ✗ inb(i)=icb(i)+1
1858 ENDDO
1859
1860 ✗ DO k = 2, nl
1861 ✗ DO i = 1, ncum
1862 ✗ IF ((k>=iposit(i))) THEN
1863 ✗ deltap = min(plcl(i), ph(i,k-1)) - min(plcl(i), ph(i,k))
1864 ✗ cape(i) = cape(i) + rrd*buoy(i, k-1)*deltap/p(i, k-1)
1865 ✗ IF (cape(i).gt.0.) THEN
1866 ✗ inb(i) = max(inb(i), k)
1867 END IF
1868 ENDIF
1869 ENDDO
1870 ENDDO
1871
1872 ! DO i = 1, ncum
1873 ! print*,"inb",inb(i)
1874 ! ENDDO
1875
1876 endif
1877
1878 ! -- end convect3
1879
1880 ! ori do 510 i=1,ncum
1881 ! ori cape(i)=0.0
1882 ! ori capem(i)=0.0
1883 ! ori inb(i)=icb(i)+1
1884 ! ori inb1(i)=inb(i)
1885 ! ori 510 continue
1886
1887 ! Originial Code
1888
1889 ! do 530 k=minorig+1,nl-1
1890 ! do 520 i=1,ncum
1891 ! if(k.ge.(icb(i)+1))then
1892 ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1893 ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1894 ! cape(i)=cape(i)+by
1895 ! if(by.ge.0.0)inb1(i)=k+1
1896 ! if(cape(i).gt.0.0)then
1897 ! inb(i)=k+1
1898 ! capem(i)=cape(i)
1899 ! endif
1900 ! endif
1901 !520 continue
1902 !530 continue
1903 ! do 540 i=1,ncum
1904 ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
1905 ! cape(i)=capem(i)+byp
1906 ! defrac=capem(i)-cape(i)
1907 ! defrac=max(defrac,0.001)
1908 ! frac(i)=-cape(i)/defrac
1909 ! frac(i)=min(frac(i),1.0)
1910 ! frac(i)=max(frac(i),0.0)
1911 !540 continue
1912
1913 ! K Emanuel fix
1914
1915 ! call zilch(byp,ncum)
1916 ! do 530 k=minorig+1,nl-1
1917 ! do 520 i=1,ncum
1918 ! if(k.ge.(icb(i)+1))then
1919 ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1920 ! cape(i)=cape(i)+by
1921 ! if(by.ge.0.0)inb1(i)=k+1
1922 ! if(cape(i).gt.0.0)then
1923 ! inb(i)=k+1
1924 ! capem(i)=cape(i)
1925 ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1926 ! endif
1927 ! endif
1928 !520 continue
1929 !530 continue
1930 ! do 540 i=1,ncum
1931 ! inb(i)=max(inb(i),inb1(i))
1932 ! cape(i)=capem(i)+byp(i)
1933 ! defrac=capem(i)-cape(i)
1934 ! defrac=max(defrac,0.001)
1935 ! frac(i)=-cape(i)/defrac
1936 ! frac(i)=min(frac(i),1.0)
1937 ! frac(i)=max(frac(i),0.0)
1938 !540 continue
1939
1940 ! J Teixeira fix
1941
1942 ! ori call zilch(byp,ncum)
1943 ! ori do 515 i=1,ncum
1944 ! ori lcape(i)=.true.
1945 ! ori 515 continue
1946 ! ori do 530 k=minorig+1,nl-1
1947 ! ori do 520 i=1,ncum
1948 ! ori if(cape(i).lt.0.0)lcape(i)=.false.
1949 ! ori if((k.ge.(icb(i)+1)).and.lcape(i))then
1950 ! ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1951 ! ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1952 ! ori cape(i)=cape(i)+by
1953 ! ori if(by.ge.0.0)inb1(i)=k+1
1954 ! ori if(cape(i).gt.0.0)then
1955 ! ori inb(i)=k+1
1956 ! ori capem(i)=cape(i)
1957 ! ori endif
1958 ! ori endif
1959 ! ori 520 continue
1960 ! ori 530 continue
1961 ! ori do 540 i=1,ncum
1962 ! ori cape(i)=capem(i)+byp(i)
1963 ! ori defrac=capem(i)-cape(i)
1964 ! ori defrac=max(defrac,0.001)
1965 ! ori frac(i)=-cape(i)/defrac
1966 ! ori frac(i)=min(frac(i),1.0)
1967 ! ori frac(i)=max(frac(i),0.0)
1968 ! ori 540 continue
1969
1970 ! --------------------------------------------------------------------
1971 ! Prevent convection when top is too hot
1972 ! --------------------------------------------------------------------
1973
2/2
✓ Branch 0 taken 127097 times.
✓ Branch 1 taken 240 times.
 127337 DO i = 1,ncum
1974
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 127097 times.
 127337 IF (t(i,inb(i)) > T_top_max) iflag(i) = 10
1975 ENDDO
1976
1977 ! =====================================================================
1978 ! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
1979 ! =====================================================================
1980
1981
2/2
✓ Branch 0 taken 6480 times.
✓ Branch 1 taken 240 times.
 6720 DO k = 1, nl
1982
2/2
✓ Branch 0 taken 3431619 times.
✓ Branch 1 taken 6480 times.
 3438339 DO i = 1, ncum
1983 3438099 hp(i, k) = h(i, k)
1984 END DO
1985 END DO
1986
1987 !jyg : cvflag_ice test outside the loops (07042015)
1988 !
1989
1/2
✓ Branch 0 taken 240 times.
✗ Branch 1 not taken.
 240 IF (cvflag_ice) THEN
1990 !
1991
1/2
✓ Branch 0 taken 240 times.
✗ Branch 1 not taken.
 240 IF (cvflag_prec_eject) THEN
1992 !! DO k = minorig + 1, nl
1993 !! DO i = 1, ncum
1994 !! IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
1995 !! frac_s(i,k) = qi(i,k)/max(clw(i,k),smallestreal)
1996 !! frac_s(i,k) = 1. - (qpl(i,k)+(1.-frac_s(i,k))*qcld(i,k))/max(clw(i,k),smallestreal)
1997 !! END IF
1998 !! END DO
1999 !! END DO
2000 ELSE ! (cvflag_prec_eject)
2001
2/2
✓ Branch 0 taken 6240 times.
✓ Branch 1 taken 240 times.
 6480 DO k = minorig + 1, nl
2002
2/2
✓ Branch 0 taken 3304522 times.
✓ Branch 1 taken 6240 times.
 3311002 DO i = 1, ncum
2003
4/4
✓ Branch 0 taken 2911168 times.
✓ Branch 1 taken 393354 times.
✓ Branch 2 taken 1144704 times.
✓ Branch 3 taken 1766464 times.
 3310762 IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
2004 !jyg< frac computation moved to beginning of cv3_undilute2.
2005 ! kept here for compatibility test with CMip6 version
2006 1144704 frac_s(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15)
2007 1144704 frac_s(i, k) = min(max(frac_s(i,k),0.0), 1.0)
2008 END IF
2009 END DO
2010 END DO
2011 ENDIF ! (cvflag_prec_eject) ELSE
2012
2/2
✓ Branch 0 taken 6240 times.
✓ Branch 1 taken 240 times.
 6480 DO k = minorig + 1, nl
2013
2/2
✓ Branch 0 taken 3304522 times.
✓ Branch 1 taken 6240 times.
 3311002 DO i = 1, ncum
2014
4/4
✓ Branch 0 taken 2911168 times.
✓ Branch 1 taken 393354 times.
✓ Branch 2 taken 1144704 times.
✓ Branch 3 taken 1766464 times.
 3310762 IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
2015 !! hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac_s(i,k)*lf(i,k))* & !!jygprl
2016 !! ep(i, k)*clw(i, k) !!jygprl
2017 hp(i, k) = hla(i,k-1) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac_s(i,k)*lf(i,k))* & !!jygprl
2018 1144704 ep(i, k)*clw(i, k) !!jygprl
2019 END IF
2020 END DO
2021 END DO
2022 !
2023 ELSE ! (cvflag_ice)
2024 !
2025 ✗ DO k = minorig + 1, nl
2026 ✗ DO i = 1, ncum
2027 ✗ IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
2028 !jyg< (energy conservation tests)
2029 !! hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*tp(i,k))*ep(i, k)*clw(i, k)
2030 !! hp(i, k) = ( hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k) ) / &
2031 !! (1. - ep(i,k)*clw(i,k))
2032 !! hp(i, k) = ( hnk(i) + (lv(i,k)+(cpd-cl)*t(i,k))*ep(i, k)*clw(i, k) ) / &
2033 !! (1. - ep(i,k)*clw(i,k))
2034 ✗ hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k)
2035 END IF
2036 END DO
2037 END DO
2038 !
2039 END IF ! (cvflag_ice)
2040
2041 240 RETURN
2042 END SUBROUTINE cv3_undilute2
2043
2044 ✗ SUBROUTINE cv3_closure(nloc, ncum, nd, icb, inb, &
2045 pbase, p, ph, tv, buoy, &
2046 ✗ sig, w0, cape, m, iflag)
2047 IMPLICIT NONE
2048
2049 ! ===================================================================
2050 ! --- CLOSURE OF CONVECT3
2051 !
2052 ! vectorization: S. Bony
2053 ! ===================================================================
2054
2055 include "cvthermo.h"
2056 include "cv3param.h"
2057
2058 !input:
2059 INTEGER ncum, nd, nloc
2060 INTEGER icb(nloc), inb(nloc)
2061 REAL pbase(nloc)
2062 REAL p(nloc, nd), ph(nloc, nd+1)
2063 REAL tv(nloc, nd), buoy(nloc, nd)
2064
2065 !input/output:
2066 REAL sig(nloc, nd), w0(nloc, nd)
2067 INTEGER iflag(nloc)
2068
2069 !output:
2070 REAL cape(nloc)
2071 REAL m(nloc, nd)
2072
2073 !local variables:
2074 INTEGER i, j, k, icbmax
2075 REAL deltap, fac, w, amu
2076 ✗ REAL dtmin(nloc, nd), sigold(nloc, nd)
2077 REAL cbmflast(nloc)
2078
2079
2080 ! -------------------------------------------------------
2081 ! -- Initialization
2082 ! -------------------------------------------------------
2083
2084 ✗ DO k = 1, nl
2085 ✗ DO i = 1, ncum
2086 ✗ m(i, k) = 0.0
2087 END DO
2088 END DO
2089
2090 ! -------------------------------------------------------
2091 ! -- Reset sig(i) and w0(i) for i>inb and i<icb
2092 ! -------------------------------------------------------
2093
2094 ! update sig and w0 above LNB:
2095
2096 ✗ DO k = 1, nl - 1
2097 ✗ DO i = 1, ncum
2098 ✗ IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
2099 sig(i, k) = beta*sig(i, k) + &
2100 ✗ 2.*alpha*buoy(i, inb(i))*abs(buoy(i,inb(i)))
2101 ✗ sig(i, k) = amax1(sig(i,k), 0.0)
2102 ✗ w0(i, k) = beta*w0(i, k)
2103 END IF
2104 END DO
2105 END DO
2106
2107 ! compute icbmax:
2108
2109 icbmax = 2
2110 ✗ DO i = 1, ncum
2111 ✗ icbmax = max(icbmax, icb(i))
2112 END DO
2113
2114 ! update sig and w0 below cloud base:
2115
2116 ✗ DO k = 1, icbmax
2117 ✗ DO i = 1, ncum
2118 ✗ IF (k<=icb(i)) THEN
2119 sig(i, k) = beta*sig(i, k) - &
2120 ✗ 2.*alpha*buoy(i, icb(i))*buoy(i, icb(i))
2121 ✗ sig(i, k) = max(sig(i,k), 0.0)
2122 ✗ w0(i, k) = beta*w0(i, k)
2123 END IF
2124 END DO
2125 END DO
2126
2127 !! if(inb.lt.(nl-1))then
2128 !! do 85 i=inb+1,nl-1
2129 !! sig(i)=beta*sig(i)+2.*alpha*buoy(inb)*
2130 !! 1 abs(buoy(inb))
2131 !! sig(i)=max(sig(i),0.0)
2132 !! w0(i)=beta*w0(i)
2133 !! 85 continue
2134 !! end if
2135
2136 !! do 87 i=1,icb
2137 !! sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb)
2138 !! sig(i)=max(sig(i),0.0)
2139 !! w0(i)=beta*w0(i)
2140 !! 87 continue
2141
2142 ! -------------------------------------------------------------
2143 ! -- Reset fractional areas of updrafts and w0 at initial time
2144 ! -- and after 10 time steps of no convection
2145 ! -------------------------------------------------------------
2146
2147 ✗ DO k = 1, nl - 1
2148 ✗ DO i = 1, ncum
2149 ✗ IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN
2150 ✗ sig(i, k) = 0.0
2151 ✗ w0(i, k) = 0.0
2152 END IF
2153 END DO
2154 END DO
2155
2156 ! -------------------------------------------------------------
2157 ! -- Calculate convective available potential energy (cape),
2158 ! -- vertical velocity (w), fractional area covered by
2159 ! -- undilute updraft (sig), and updraft mass flux (m)
2160 ! -------------------------------------------------------------
2161
2162 ✗ DO i = 1, ncum
2163 ✗ cape(i) = 0.0
2164 END DO
2165
2166 ! compute dtmin (minimum buoyancy between ICB and given level k):
2167
2168 ✗ DO i = 1, ncum
2169 ✗ DO k = 1, nl
2170 ✗ dtmin(i, k) = 100.0
2171 END DO
2172 END DO
2173
2174 ✗ DO i = 1, ncum
2175 ✗ DO k = 1, nl
2176 ✗ DO j = minorig, nl
2177 ✗ IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k-1))) THEN
2178 ✗ dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
2179 END IF
2180 END DO
2181 END DO
2182 END DO
2183
2184 ! the interval on which cape is computed starts at pbase :
2185
2186 ✗ DO k = 1, nl
2187 ✗ DO i = 1, ncum
2188
2189 ✗ IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
2190
2191 ✗ deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
2192 ✗ cape(i) = cape(i) + rrd*buoy(i, k-1)*deltap/p(i, k-1)
2193 ✗ cape(i) = amax1(0.0, cape(i))
2194 ✗ sigold(i, k) = sig(i, k)
2195
2196 ! dtmin(i,k)=100.0
2197 ! do 97 j=icb(i),k-1 ! mauvaise vectorisation
2198 ! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
2199 ! 97 continue
2200
2201 ✗ sig(i, k) = beta*sig(i, k) + alpha*dtmin(i, k)*abs(dtmin(i,k))
2202 ✗ sig(i, k) = max(sig(i,k), 0.0)
2203 ✗ sig(i, k) = amin1(sig(i,k), 0.01)
2204 ✗ fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
2205 ✗ w = (1.-beta)*fac*sqrt(cape(i)) + beta*w0(i, k)
2206 ✗ amu = 0.5*(sig(i,k)+sigold(i,k))*w
2207 ✗ m(i, k) = amu*0.007*p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
2208 ✗ w0(i, k) = w
2209 END IF
2210
2211 END DO
2212 END DO
2213
2214 ✗ DO i = 1, ncum
2215 ✗ w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
2216 ✗ m(i, icb(i)) = 0.5*m(i, icb(i)+1)*(ph(i,icb(i))-ph(i,icb(i)+1))/(ph(i,icb(i)+1)-ph(i,icb(i)+2))
2217 ✗ sig(i, icb(i)) = sig(i, icb(i)+1)
2218 ✗ sig(i, icb(i)-1) = sig(i, icb(i))
2219 END DO
2220
2221 ! ccc 3. Compute final cloud base mass flux and set iflag to 3 if
2222 ! ccc cloud base mass flux is exceedingly small and is decreasing (i.e. if
2223 ! ccc the final mass flux (cbmflast) is greater than the target mass flux
2224 ! ccc (cbmf) ??).
2225 ! cc
2226 ! c do i = 1,ncum
2227 ! c cbmflast(i) = 0.
2228 ! c enddo
2229 ! cc
2230 ! c do k= 1,nl
2231 ! c do i = 1,ncum
2232 ! c IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
2233 ! c cbmflast(i) = cbmflast(i)+M(i,k)
2234 ! c ENDIF
2235 ! c enddo
2236 ! c enddo
2237 ! cc
2238 ! c do i = 1,ncum
2239 ! c IF (cbmflast(i) .lt. 1.e-6) THEN
2240 ! c iflag(i) = 3
2241 ! c ENDIF
2242 ! c enddo
2243 ! cc
2244 ! c do k= 1,nl
2245 ! c do i = 1,ncum
2246 ! c IF (iflag(i) .ge. 3) THEN
2247 ! c M(i,k) = 0.
2248 ! c sig(i,k) = 0.
2249 ! c w0(i,k) = 0.
2250 ! c ENDIF
2251 ! c enddo
2252 ! c enddo
2253 ! cc
2254 !! cape=0.0
2255 !! do 98 i=icb+1,inb
2256 !! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
2257 !! cape=cape+rrd*buoy(i-1)*deltap/p(i-1)
2258 !! dcape=rrd*buoy(i-1)*deltap/p(i-1)
2259 !! dlnp=deltap/p(i-1)
2260 !! cape=max(0.0,cape)
2261 !! sigold=sig(i)
2262
2263 !! dtmin=100.0
2264 !! do 97 j=icb,i-1
2265 !! dtmin=amin1(dtmin,buoy(j))
2266 !! 97 continue
2267
2268 !! sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin)
2269 !! sig(i)=max(sig(i),0.0)
2270 !! sig(i)=amin1(sig(i),0.01)
2271 !! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
2272 !! w=(1.-beta)*fac*sqrt(cape)+beta*w0(i)
2273 !! amu=0.5*(sig(i)+sigold)*w
2274 !! m(i)=amu*0.007*p(i)*(ph(i)-ph(i+1))/tv(i)
2275 !! w0(i)=w
2276 !! 98 continue
2277 !! w0(icb)=0.5*w0(icb+1)
2278 !! m(icb)=0.5*m(icb+1)*(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
2279 !! sig(icb)=sig(icb+1)
2280 !! sig(icb-1)=sig(icb)
2281
2282 ✗ RETURN
2283 END SUBROUTINE cv3_closure
2284
2285 ✗ SUBROUTINE cv3_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
2286 ✗ ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
2287 ✗ unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
2288 ✗ ment, qent, uent, vent, nent, sij, elij, ments, qents, traent)
2289 IMPLICIT NONE
2290
2291 ! ---------------------------------------------------------------------
2292 ! a faire:
2293 ! - vectorisation de la partie normalisation des flux (do 789...)
2294 ! ---------------------------------------------------------------------
2295
2296 include "cvthermo.h"
2297 include "cv3param.h"
2298 include "cvflag.h"
2299
2300 !inputs:
2301 INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
2302 INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
2303 REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig
2304 REAL, DIMENSION (nloc), INTENT (IN) :: qnk, unk, vnk
2305 REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
2306 REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
2307 REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
2308 REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra ! input of convect3
2309 REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, h, hp
2310 REAL, DIMENSION (nloc, na), INTENT (IN) :: lf, frac
2311 REAL, DIMENSION (nloc, na), INTENT (IN) :: tv, tvp, ep, clw
2312 REAL, DIMENSION (nloc, na), INTENT (IN) :: m ! input of convect3
2313
2314 !outputs:
2315 REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: ment, qent
2316 REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: uent, vent
2317 REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: sij, elij
2318 REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
2319 REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: ments, qents
2320 INTEGER, DIMENSION (nloc, nd), INTENT (OUT) :: nent
2321
2322 !local variables:
2323 INTEGER i, j, k, il, im, jm
2324 INTEGER num1, num2
2325 REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
2326 REAL alt, smid, sjmin, sjmax, delp, delm
2327 ✗ REAL asij(nloc), smax(nloc), scrit(nloc)
2328 ✗ REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
2329 ✗ REAL sigij(nloc, nd, nd)
2330 REAL wgh
2331 ✗ REAL zm(nloc, na)
2332 ✗ LOGICAL lwork(nloc)
2333
2334 ! =====================================================================
2335 ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
2336 ! =====================================================================
2337
2338 ! ori do 360 i=1,ncum*nlp
2339 ✗ DO j = 1, nl
2340 ✗ DO i = 1, ncum
2341 ✗ nent(i, j) = 0
2342 ! in convect3, m is computed in cv3_closure
2343 ! ori m(i,1)=0.0
2344 END DO
2345 END DO
2346
2347 ! ori do 400 k=1,nlp
2348 ! ori do 390 j=1,nlp
2349 ✗ DO j = 1, nl
2350 ✗ DO k = 1, nl
2351 ✗ DO i = 1, ncum
2352 ✗ qent(i, k, j) = rr(i, j)
2353 ✗ uent(i, k, j) = u(i, j)
2354 ✗ vent(i, k, j) = v(i, j)
2355 ✗ elij(i, k, j) = 0.0
2356 !ym ment(i,k,j)=0.0
2357 !ym sij(i,k,j)=0.0
2358 END DO
2359 END DO
2360 END DO
2361
2362 !ym
2363 ✗ ment(1:ncum, 1:nd, 1:nd) = 0.0
2364 ✗ sij(1:ncum, 1:nd, 1:nd) = 0.0
2365
2366 !AC! do k=1,ntra
2367 !AC! do j=1,nd ! instead nlp
2368 !AC! do i=1,nd ! instead nlp
2369 !AC! do il=1,ncum
2370 !AC! traent(il,i,j,k)=tra(il,j,k)
2371 !AC! enddo
2372 !AC! enddo
2373 !AC! enddo
2374 !AC! enddo
2375 ✗ zm(:, :) = 0.
2376
2377 ! =====================================================================
2378 ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
2379 ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
2380 ! --- FRACTION (sij)
2381 ! =====================================================================
2382
2383 ✗ DO i = minorig + 1, nl
2384
2385 ✗ DO j = minorig, nl
2386 ✗ DO il = 1, ncum
2387 ✗ IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) .AND. (j<=inb(il))) THEN
2388
2389 ✗ rti = qnk(il) - ep(il, i)*clw(il, i)
2390 ✗ bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
2391
2392
2393 ✗ IF (cvflag_ice) THEN
2394 ! print*,cvflag_ice,'cvflag_ice dans do 700'
2395 ✗ IF (t(il,j)<=263.15) THEN
2396 bf2 = 1. + (lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)* &
2397 ✗ lf(il,j))*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
2398 END IF
2399 END IF
2400
2401 ✗ anum = h(il, j) - hp(il, i) + (cpv-cpd)*t(il, j)*(rti-rr(il,j))
2402 ✗ denom = h(il, i) - hp(il, i) + (cpd-cpv)*(rr(il,i)-rti)*t(il, j)
2403 dei = denom
2404 ✗ IF (abs(dei)<0.01) dei = 0.01
2405 ✗ sij(il, i, j) = anum/dei
2406 ✗ sij(il, i, i) = 1.0
2407 ✗ altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - rs(il, j)
2408 ✗ altem = altem/bf2
2409 ✗ cwat = clw(il, j)*(1.-ep(il,j))
2410 stemp = sij(il, i, j)
2411 ✗ IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
2412
2413 ✗ IF (cvflag_ice) THEN
2414 ✗ anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat*bf2)
2415 ✗ denom = denom + (lv(il,j)+frac(il,j)*lf(il,j))*(rr(il,i)-rti)
2416 ELSE
2417 ✗ anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
2418 ✗ denom = denom + lv(il, j)*(rr(il,i)-rti)
2419 END IF
2420
2421 ✗ IF (abs(denom)<0.01) denom = 0.01
2422 ✗ sij(il, i, j) = anum/denom
2423 ✗ altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - rs(il, j)
2424 ✗ altem = altem - (bf2-1.)*cwat
2425 END IF
2426 ✗ IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
2427 ✗ qent(il, i, j) = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti
2428 ✗ uent(il, i, j) = sij(il, i, j)*u(il, i) + (1.-sij(il,i,j))*unk(il)
2429 ✗ vent(il, i, j) = sij(il, i, j)*v(il, i) + (1.-sij(il,i,j))*vnk(il)
2430 !!!! do k=1,ntra
2431 !!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
2432 !!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
2433 !!!! end do
2434 elij(il, i, j) = altem
2435 ✗ elij(il, i, j) = max(0.0, elij(il,i,j))
2436 ✗ ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
2437 ✗ nent(il, i) = nent(il, i) + 1
2438 END IF
2439 ✗ sij(il, i, j) = max(0.0, sij(il,i,j))
2440 ✗ sij(il, i, j) = amin1(1.0, sij(il,i,j))
2441 END IF ! new
2442 END DO
2443 END DO
2444
2445 !AC! do k=1,ntra
2446 !AC! do j=minorig,nl
2447 !AC! do il=1,ncum
2448 !AC! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
2449 !AC! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
2450 !AC! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
2451 !AC! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
2452 !AC! endif
2453 !AC! enddo
2454 !AC! enddo
2455 !AC! enddo
2456
2457
2458 ! *** if no air can entrain at level i assume that updraft detrains ***
2459 ! *** at that level and calculate detrained air flux and properties ***
2460
2461
2462 ! @ do 170 i=icb(il),inb(il)
2463
2464 ✗ DO il = 1, ncum
2465 ✗ IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
2466 ! @ if(nent(il,i).eq.0)then
2467 ✗ ment(il, i, i) = m(il, i)
2468 ✗ qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
2469 ✗ uent(il, i, i) = unk(il)
2470 ✗ vent(il, i, i) = vnk(il)
2471 ✗ elij(il, i, i) = clw(il, i)
2472 ! MAF sij(il,i,i)=1.0
2473 ✗ sij(il, i, i) = 0.0
2474 END IF
2475 END DO
2476 END DO
2477
2478 !AC! do j=1,ntra
2479 !AC! do i=minorig+1,nl
2480 !AC! do il=1,ncum
2481 !AC! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
2482 !AC! traent(il,i,i,j)=tra(il,nk(il),j)
2483 !AC! endif
2484 !AC! enddo
2485 !AC! enddo
2486 !AC! enddo
2487
2488 ✗ DO j = minorig, nl
2489 ✗ DO i = minorig, nl
2490 ✗ DO il = 1, ncum
2491 ✗ IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<=inb(il))) THEN
2492 ✗ sigij(il, i, j) = sij(il, i, j)
2493 END IF
2494 END DO
2495 END DO
2496 END DO
2497 ! @ enddo
2498
2499 ! @170 continue
2500
2501 ! =====================================================================
2502 ! --- NORMALIZE ENTRAINED AIR MASS FLUXES
2503 ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
2504 ! =====================================================================
2505
2506 ✗ CALL zilch(asum, nloc*nd)
2507 ✗ CALL zilch(csum, nloc*nd)
2508 ✗ CALL zilch(csum, nloc*nd)
2509
2510 ✗ DO il = 1, ncum
2511 ✗ lwork(il) = .FALSE.
2512 END DO
2513
2514 ✗ DO i = minorig + 1, nl
2515
2516 num1 = 0
2517 ✗ DO il = 1, ncum
2518 ✗ IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
2519 END DO
2520 ✗ IF (num1<=0) GO TO 789
2521
2522
2523 ✗ DO il = 1, ncum
2524