Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	phylmd/cv3_routines.F90	Lines:	1067	1721	62.0 %
Date:	2023-06-30 12:56:34	Branches:	920	1656	55.6 %

Line	Branch	Exec	Source
1
2			! $Id: cv3_routines.F90 4076 2022-02-04 08:37:44Z jyg $
3
4
5
6
7		288	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		144	noff = min(max(nd-k_upper, 1), (nd+1)/2)
60			!! noff = 1
61			!>jyg
62		144	minorig = 1
63		144	nl = nd - noff
64		144	nlp = nl + 1
65		144	nlm = nl - 1
66
67	✓✓	144	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		144	beta = 1.0 - delt/tau
170			alpha1 = 1.5E-3
171			!JYG Correction bug alpha
172		144	alpha1 = alpha1*1.5
173		144	alpha = alpha1*delt/tau
174			!JYG Bug
175			! cc increase alpha to compensate W decrease:
176			! c alpha = alpha*1.5
177
178		144	noconv_stop = max(2.,tau_stop/delt)
179
180		144	RETURN
181			END SUBROUTINE cv3_param
182
183		144	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	✓✗	144	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	✓✓	143280	DO il = 1, len
214		143136	sig(il, nd) = sig(il, nd) + 1.
215		143280	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		144	RETURN
221			END SUBROUTINE cv3_incrcount
222
223		15458976	SUBROUTINE cv3_prelim(len, nd, ndp1, t, q, p, ph, &
224		288	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		576	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	✓✓	8352	DO k = 1, nlp
255
256	✓✓	8023968	DO i = 1, len
257			! debug lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
258		8015616	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		8015616	lf(i, k) = lf0 + clmci*(t(i,k)-273.15)
261		8015616	cpn(i, k) = cpd(1.0-q(i,k)) + cpvq(i, k)
262		8015616	cpx(i, k) = cpd(1.0-q(i,k)) + clq(i, k)
263			! ori tv(i,k)=t(i,k)(1.0+q(i,k)epsim1)
264		8015616	tv(i, k) = t(i, k)*(1.0+q(i,k)/eps-q(i,k))
265		8015616	rdcp = (rrd(1.-q(i,k))+q(i,k)rrv)/cpn(i, k)
266		8023680	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	✓✓✓✓	11176128	gz(:,:) = 0. !jyg: initialization of the whole array
276			! ori do 140 k=2,nlp
277	✓✓	7776	DO k = 2, nl ! convect3
278	✓✓	7450848	DO i = 1, len
279		7443072	tvx = t(i, k)*(1.+q(i,k)/eps-q(i,k)) !convect3
280		7443072	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.5rrd(tvx+tvy)* & !convect3
282		7450560	(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	✓✓	8064	DO k = 1, nl ! convect3
296	✓✓	7737408	DO i = 1, len
297		7729344	h(i, k) = gz(i, k) + cpn(i, k)*t(i, k)
298		7737120	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		288	RETURN
303			END SUBROUTINE cv3_prelim
304
305		144	SUBROUTINE cv3_feed(len, nd, ok_conserv_q, &
306		144	t, q, u, v, p, ph, h, gz, &
307			p1feed, p2feed, wght, &
308			wghti, tnk, thnk, qnk, qsnk, unk, vnk, &
309		144	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		288	REAL pup(len), plo(len), pfeed(len)
359		288	REAL plclup(len), plcllo(len), plclfeed(len)
360		288	REAL pfeedmin(len)
361		288	REAL posit(len)
362		288	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	✓✓	144	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	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	✓✓	143280	DO i = 1, len
393		143136	nk(i) = minorig
394		143280	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	✓✓	143280	DO i = 1, len
408		143280	pup(i) = p2feed(i)
409			END DO
410	✗✓	144	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		144	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	✓✓	143280	DO i = 1, len
421		143280	plo(i) = ph(i, nk(i)+1)
422			END DO
423	✗✓	144	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		144	wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plcllo)
431			ENDIF ! (fl_cor_ebil >=2 )
432			! 2- Iterations
433			niter = 5
434	✓✓	864	DO iter = 1, niter
435	✓✓	716400	DO i = 1, len
436		715680	plcllo(i) = min(plo(i), plcllo(i))
437		715680	plclup(i) = max(pup(i), plclup(i))
438		716400	nocond(i) = plclup(i) <= pup(i)
439			END DO
440	✓✓	716400	DO i = 1, len
441	✓✓	716400	IF (nocond(i)) THEN
442		502765	pfeed(i) = pup(i)
443			ELSE
444			!JYG20140217<
445	✓✗	212915	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		212915	(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	✓✗	720	IF (ok_new_feed) THEN
461	✓✓	720	IF (iter==niter) THEN
462	✓✓	143280	DO i = 1,len !jyg
463		143280	pfeedmin(i) = ph(i,minorig+1) !jyg
464			ENDDO !jyg
465	✓✓	3888	DO k = minorig+1, nl !jyg
466			!! DO k = minorig, nl !jyg
467	✓✓	3725424	DO i = 1, len
468	✓✓	3725280	IF (ph(i,k)>=plclfeed(i)) pfeedmin(i) = ph(i, k)
469			END DO
470			END DO
471	✓✓	143280	DO i = 1, len
472		143280	pfeed(i) = max(pfeedmin(i), pfeed(i))
473			END DO
474			END IF
475			END IF
476			!jyg>
477
478	✗✓	720	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		720	wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclfeed)
486			ENDIF ! (fl_cor_ebil >=2 )
487			!jyg20140217<
488	✓✗	720	IF (ok_new_feed) THEN
489	✓✓	716400	DO i = 1, len
490		715680	posit(i) = (sign(1.,plclfeed(i)-pfeed(i)+dp_lcl_feed)+1.)*0.5
491	✓✓	716400	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	✓✓	716544	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		715680	pup(i) = posit(i)pfeed(i) + (1.-posit(i))pup(i)
506		715680	plo(i) = (1.-posit(i))pfeed(i) + posit(i)plo(i)
507		715680	plclup(i) = posit(i)plclfeed(i) + (1.-posit(i))plclup(i)
508		716400	plcllo(i) = (1.-posit(i))plclfeed(i) + posit(i)plcllo(i)
509			END DO
510			END DO ! iter
511
512	✓✓	143280	DO i = 1, len
513		143136	p2feed(i) = pfeed(i)
514		143280	plcl(i) = plclfeed(i)
515			END DO
516
517	✓✓	143280	DO i = 1, len
518		143136	cpnk(i) = cpd(1.0-qnk(i)) + cpvqnk(i)
519		143280	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	✗✓	144	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	✓✗	144	ELSEIF (cv_flag_feed >= 2) THEN
533			! --- and demand that LCL be high enough
534	✓✓	143280	DO i = 1, len
535			IF (((tnk(i)<250.0) .OR. &
536			(qnk(i)<=0.0) .OR. &
537	✓✓✓✗ ✓✓✓✗	143136	(plcl(i)>min(0.99*ph(i,1),ph(i,3)))) .AND. &
538		55628	(iflag(i)==0)) iflag(i) = 7
539			END DO
540			ENDIF
541	✗✓	144	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	✓✓	143280	DO i = 1, len
566		143280	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	✓✓	3600	DO k = 3, nl - 1 ! modif pour que icb soit sup/egal a 2
573	✓✓	3438864	DO i = 1, len
574	✓✓	3438720	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	✓✓	143280	DO i = 1, len
584			!@ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
585	✗✓✗✗	143280	IF ((icb(i)==nlm) .AND. (iflag(i)==0)) iflag(i) = 9
586			END DO
587
588	✓✓	143280	DO i = 1, len
589		143280	icb(i) = icb(i) - 1 ! icb sup ou egal a 2
590			END DO
591
592			! Compute icbmax.
593
594		144	icbmax = 2
595	✓✓	143280	DO i = 1, len
596			!! icbmax=max(icbmax,icb(i))
597	✓✓	143280	IF (iflag(i)<7) icbmax = max(icbmax, icb(i)) ! sb Jun7th02
598			END DO
599
600		144	RETURN
601			END SUBROUTINE cv3_feed
602
603		1996096	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		288	INTEGER icb1(len), icbsmax2 ! convect3
638			REAL tg, qg, alv, s, ahg, tc, denom, es, rg
639		288	REAL ah0(len), cpp(len)
640		288	REAL ticb(len), gzicb(len)
641		288	REAL qsicb(len) ! convect3
642		144	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			! cptp+Lqp+phi=cptnk+Lqnk+gznk.
649			! -------------------------------------------------------------------
650
651
652			! * Calculate certain parcel quantities, including static energy *
653
654	✓✓	143280	DO i = 1, len
655		143136	ah0(i) = (cpd(1.-qnk(i))+clqnk(i))tnk(i) + qnk(i)(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
656		143136	cpp(i) = cpd(1.-qnk(i)) + qnk(i)cpv
657		143280	cpinv(i) = 1./cpp(i)
658			END DO
659
660			! * Calculate lifted parcel quantities below cloud base *
661
662	✓✓	143280	DO i = 1, len !convect3
663		143136	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		143136	icbs(i) = icb1(i) !convect3
667	✓✓	143280	IF (plcl(i)<p(i,icb1(i))) THEN
668		52252	icbs(i) = min(icbs(i)+1, nl) !convect3
669			END IF
670			END DO !convect3
671
672	✓✓	143280	DO i = 1, len !convect3
673		143136	ticb(i) = t(i, icbs(i)) !convect3
674		143136	gzicb(i) = gz(i, icbs(i)) !convect3
675		143280	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	✓✓	143280	DO i = 1, len !convect3
683		143280	icbsmax2 = max(icbsmax2, icbs(i)) !convect3
684			END DO !convect3
685
686			! initialization outputs:
687
688	✓✓	1720	DO k = 1, icbsmax2 ! convect3
689	✓✓	1568264	DO i = 1, len ! convect3
690		1566544	tp(i, k) = 0.0 ! convect3
691		1566544	tvp(i, k) = 0.0 ! convect3
692		1568120	clw(i, k) = 0.0 ! convect3
693			END DO ! convect3
694			END DO ! convect3
695
696			! tp and tvp below cloud base:
697
698	✓✓	1576	DO k = minorig, icbsmax2 - 1
699	✓✓	1424984	DO i = 1, len
700		1423408	tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))*cpinv(i)
701		1424840	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	✓✓	143280	DO i = 1, len
708		143136	tg = ticb(i)
709			! ori qg=qs(i,icb(i))
710		143136	qg = qsicb(i) ! convect3
711			! debug alv=lv0-clmcpv*(ticb(i)-t0)
712		143136	alv = lv0 - clmcpv*(ticb(i)-273.15)
713
714			! First iteration.
715
716			! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
717			s = cpd(1.-qnk(i)) + clqnk(i) + & ! convect3
718		143136	alvalvqg/(rrvticb(i)ticb(i)) ! convect3
719		143136	s = 1./s
720			! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
721		143136	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
722		143136	tg = tg + s*(ah0(i)-ahg)
723			! ori tg=max(tg,35.0)
724			! debug tc=tg-t0
725		143136	tc = tg - 273.15
726		143136	denom = 243.5 + tc
727		143136	denom = max(denom, 1.0) ! convect3
728			! ori if(tc.ge.0.0)then
729		143136	es = 6.112exp(17.67tc/denom)
730			! ori else
731			! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
732			! ori endif
733			! ori qg=epses/(p(i,icb(i))-es(1.-eps))
734		143136	qg = epses/(p(i,icbs(i))-es(1.-eps))
735
736			! Second iteration.
737
738
739			! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
740			! ori s=1./s
741			! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
742		143136	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
743		143136	tg = tg + s*(ah0(i)-ahg)
744			! ori tg=max(tg,35.0)
745			! debug tc=tg-t0
746		143136	tc = tg - 273.15
747		143136	denom = 243.5 + tc
748		143136	denom = max(denom, 1.0) ! convect3
749			! ori if(tc.ge.0.0)then
750		143136	es = 6.112exp(17.67tc/denom)
751			! ori else
752			! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
753			! ori end if
754			! ori qg=epses/(p(i,icb(i))-es(1.-eps))
755		143136	qg = epses/(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		143136	tp(i, icbs(i)) = (ah0(i)-gz(i,icbs(i))-alvqg)/(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		143136	clw(i, icbs(i)) = qnk(i) - qg
769		143136	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.+rgepsi)
773			! convect3: (qg utilise au lieu du vrai mixing ratio rg)
774		143280	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	✓✓	143280	DO i = 1, len
801		143136	ticb(i) = t(i, icb(i)+1)
802		143136	gzicb(i) = gz(i, icb(i)+1)
803		143280	qsicb(i) = qs(i, icb(i)+1)
804			END DO
805
806	✓✓	143280	DO i = 1, len
807		143136	tg = ticb(i)
808		143136	qg = qsicb(i) ! convect3
809			! debug alv=lv0-clmcpv*(ticb(i)-t0)
810		143136	alv = lv0 - clmcpv*(ticb(i)-273.15)
811
812			! First iteration.
813
814			! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
815			s = cpd(1.-qnk(i)) + clqnk(i) & ! convect3
816		143136	+alvalvqg/(rrvticb(i)ticb(i)) ! convect3
817		143136	s = 1./s
818			! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
819		143136	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
820		143136	tg = tg + s*(ah0(i)-ahg)
821			! ori tg=max(tg,35.0)
822			! debug tc=tg-t0
823		143136	tc = tg - 273.15
824		143136	denom = 243.5 + tc
825		143136	denom = max(denom, 1.0) ! convect3
826			! ori if(tc.ge.0.0)then
827		143136	es = 6.112exp(17.67tc/denom)
828			! ori else
829			! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
830			! ori endif
831			! ori qg=epses/(p(i,icb(i))-es(1.-eps))
832		143136	qg = epses/(p(i,icb(i)+1)-es(1.-eps))
833
834			! Second iteration.
835
836
837			! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
838			! ori s=1./s
839			! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
840		143136	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
841		143136	tg = tg + s*(ah0(i)-ahg)
842			! ori tg=max(tg,35.0)
843			! debug tc=tg-t0
844		143136	tc = tg - 273.15
845		143136	denom = 243.5 + tc
846		143136	denom = max(denom, 1.0) ! convect3
847			! ori if(tc.ge.0.0)then
848		143136	es = 6.112exp(17.67tc/denom)
849			! ori else
850			! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
851			! ori end if
852			! ori qg=epses/(p(i,icb(i))-es(1.-eps))
853		143136	qg = epses/(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		143136	tp(i, icb(i)+1) = (ah0(i)-gz(i,icb(i)+1)-alvqg)/(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		143136	clw(i, icb(i)+1) = qnk(i) - qg
867		143136	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.+rgepsi)
871			! convect3: (qg utilise au lieu du vrai mixing ratio rg)
872		143280	tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
873
874			END DO
875
876		144	RETURN
877			END SUBROUTINE cv3_undilute1
878
879		5249736	SUBROUTINE cv3_trigger(len, nd, icb, plcl, p, th, tv, tvp, thnk, &
880		144	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	✓✓	143280	DO i = 1, len
922		143136	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		143136	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		143136	tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i)) /(p(i,icb(i))-p(i,icb(i)+1))
927		143280	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) = betasig(i,k) - 2.alphatdiftdif
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	✓✓	4032	DO k = 1, nl
960	✓✓	3868704	DO i = 1, len
961
962		3864672	tdif = buoybase(i)
963		3864672	ath1 = thnk(i)
964		3864672	ath = th(i, icb(i)-1) - dttrig
965
966	✓✓✗✓	3868560	IF (tdif<dtcrit .OR. ath>ath1) THEN
967		1241784	sig(i, k) = betasig(i, k) - 2.alphatdiftdif
968		1241784	sig(i, k) = amax1(sig(i,k), 0.0)
969		1241784	w0(i, k) = beta*w0(i, k)
970		1241784	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		144	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		11044093	SUBROUTINE cv3_undilute2(nloc, ncum, nd, iflag, icb, icbs, nk, &
1132		144	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		288	REAL, DIMENSION (nloc,nd) :: qi
1196		288	REAL, DIMENSION (nloc,nd) :: ha ! moist static energy of adiabatic ascents
1197			! taking into account precip ejection
1198		288	REAL, DIMENSION (nloc,nd) :: hla ! liquid water static energy of adiabatic ascents
1199			! taking into account precip ejection
1200		288	REAL, DIMENSION (nloc,nd) :: qcld ! specific cloud water
1201		288	REAL, DIMENSION (nloc,nd) :: qhsat ! specific humidity at saturation
1202			REAL, DIMENSION (nloc,nd) :: dqhsatdT ! dqhsat/dT
1203		288	REAL, DIMENSION (nloc,nd) :: frac ! ice fraction function of envt temperature
1204		288	REAL, DIMENSION (nloc,nd) :: qps ! specific solid precipitation
1205		288	REAL, DIMENSION (nloc,nd) :: qpl ! specific liquid precipitation
1206		288	REAL, DIMENSION (nloc) :: ah0, cape, capem, byp
1207			LOGICAL, DIMENSION (nloc) :: lcape
1208		288	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	✗✓	144	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	✓✓	4032	DO k = 1, nl
1231	✓✓	1866114	DO i = 1, ncum
1232		1865970	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			! cptp+Lqp+phi=cptnk+Lqnk+gznk.
1243
1244			! * Calculate certain parcel quantities, including static energy *
1245
1246
1247	✓✓	69110	DO i = 1, ncum
1248			ah0(i) = (cpd(1.-qnk(i))+clqnk(i))*tnk(i)+ &
1249			! debug qnk(i)(lv0-clmcpv(tnk(i)-t0))+gznk(i)
1250		69110	qnk(i)(lv0-clmcpv(tnk(i)-273.15)) + gznk(i)
1251			END DO
1252			!
1253			! Ice fraction
1254			!
1255	✓✗	144	IF (cvflag_ice) THEN
1256	✓✓	4032	DO k = minorig, nl
1257	✓✓	1866114	DO i = 1, ncum
1258		1862082	frac(i, k) = (Tx - t(i,k))/(Tx - Tm)
1259		1865970	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	✓✓	4032	DO k = 1, nl
1264	✓✓	1866114	DO i = 1, ncum
1265	✓✓	1865970	IF (k<icb(i)) THEN
1266		305964	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	✓✓	4032	DO k = minorig, nl
1280	✓✓	1866114	DO i = 1,ncum
1281		1862082	ha(i,k) = ah0(i)
1282		1862082	hla(i,k) = hnk(i)
1283		1862082	qta(i,k) = qnk(i)
1284		1862082	qpreca(i,k) = 0.
1285		1862082	frac_a(i,k) = 0.
1286		1862082	frac_s(i,k) = frac(i,k)
1287		1862082	qpl(i,k) = 0.
1288		1862082	qps(i,k) = 0.
1289		1862082	qhsat(i,k) = qs(i,k)
1290		1862082	qcld(i,k) = max(qta(i,k)-qhsat(i,k),0.)
1291	✓✓	1865970	IF (k <= icb(i)+1) THEN
1292		443896	qhsat(i,k) = qnk(i)-clw(i,k)
1293		443896	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	✓✓	4032	DO k = 1, nl
1303	✓✓	1866114	DO i = 1, ncum
1304		1862082	ep(i, k) = 0.0
1305		1865970	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	✗✓	144	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.112exp(17.67(tg - 273.15)/(tg + 243.5 - 273.15))
1333			qg = epses/(p(i,k)-es(1.-eps))
1334			ELSE
1335			esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
1336			qg = epsesi/(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/(rrvtg*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/(rrvtgtg)
1353			ddelta = lf(i,k)*(qnk(i) - qg)
1354			bb = aa + ddeltafp + ddfp*(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 +fpddelta)(Tx - Tg)
1364			Zm = ahg + (coefm +fpddelta)(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 = bbbb - 4ddfp(ah0(i) - ahg + ddeltafp(Tx-tg))
1377			Tp(i,k) = tg + &
1378			Um* (ah0(i) - ahg + ddelta) /(aa + dd) + &
1379			U 2(ah0(i) - ahg + ddeltafp(Tx-tg))/(bb + sqrt(discr)) + &
1380			Ux* (ah0(i) - ahg) /aa
1381			ELSEIF (tg .gt. Tm) THEN
1382			discr = bbbb - 4ddfp(ah0(i) - ahg)
1383			Tp(i,k) = tg + &
1384			Um* (ah0(i) - ahg + ddeltafp(tg-Tm))/(aa + dd) + &
1385			U 2(ah0(i) - ahg) /(bb + sqrt(discr)) + &
1386			Ux* (ah0(i) - ahg + ddeltafp(tg-Tx))/aa
1387			ELSE
1388			discr = bbbb - 4ddfp(ah0(i) - ahg + ddeltafp(Tm-tg))
1389			Tp(i,k) = tg + &
1390			Um* (ah0(i) - ahg) /(aa + dd) + &
1391			U 2(ah0(i) - ahg + ddeltafp(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.112exp(17.67(tg - 273.15)/(tg + 243.5 - 273.15))
1405			qg = epses/(p(i,k)-es(1.-eps))
1406			ELSE
1407			esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
1408			qg = epsesi/(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	✗✓	144	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.112exp(17.67(tg - 273.15)/(tg + 243.5 - 273.15))
1441			qg = epses/(p(i,k)-es(1.-eps))
1442			dqgdT = lv(i,k)qg/(rrvtg*tg)
1443			ELSE
1444			esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
1445			qg = epsesi/(p(i,k)-esi(1.-eps))
1446			dqgdT = (lv(i,k)+lf(i,k))qg/(rrvtg*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.112exp(17.67(tg - 273.15)/(tg + 243.5 - 273.15))
1466			qg = epses/(p(i,k)-es(1.-eps))
1467			ELSE
1468			esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
1469			qg = epsesi/(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) = ejectliqqpl(i,k) + ejecticeqps(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)ejecticeqps(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)ejectliqqpl(i,k) + (lv(i,k)+lf(i,k))ejecticeqps(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	✓✗	144	ELSE IF (icvflag_Tpa == 0) THEN! (icvflag_Tpa == 2) ELSE IF(icvflag_Tpa == 1)
1545			!
1546			!----------------------------------------------------------------------------
1547			!
1548	✓✓	3888	DO k = minorig + 1, nl
1549	✓✓	1797004	DO i = 1, ncum
1550			! ori if(k.ge.(icb(i)+1))then
1551	✓✓	1796860	IF (k>=(icbs(i)+1)) THEN ! convect3
1552		1449226	tg = t(i, k)
1553		1449226	qg = qs(i, k)
1554			! debug alv=lv0-clmcpv*(t(i,k)-t0)
1555		1449226	alv = lv0 - clmcpv*(t(i,k)-273.15)
1556
1557			! First iteration.
1558
1559			! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
1560			s = cpd(1.-qnk(i)) + clqnk(i) + & ! convect3
1561		1449226	alvalvqg/(rrvt(i,k)t(i,k)) ! convect3
1562		1449226	s = 1./s
1563			! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
1564		1449226	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
1565		1449226	tg = tg + s*(ah0(i)-ahg)
1566			! ori tg=max(tg,35.0)
1567			! debug tc=tg-t0
1568		1449226	tc = tg - 273.15
1569		1449226	denom = 243.5 + tc
1570		1449226	denom = max(denom, 1.0) ! convect3
1571			! ori if(tc.ge.0.0)then
1572		1449226	es = 6.112exp(17.67tc/denom)
1573			! ori else
1574			! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
1575			! ori endif
1576		1449226	qg = epses/(p(i,k)-es(1.-eps))
1577
1578			! Second iteration.
1579
1580			! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
1581			! ori s=1./s
1582			! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
1583		1449226	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
1584		1449226	tg = tg + s*(ah0(i)-ahg)
1585			! ori tg=max(tg,35.0)
1586			! debug tc=tg-t0
1587		1449226	tc = tg - 273.15
1588		1449226	denom = 243.5 + tc
1589		1449226	denom = max(denom, 1.0) ! convect3
1590			! ori if(tc.ge.0.0)then
1591		1449226	es = 6.112exp(17.67tc/denom)
1592			! ori else
1593			! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
1594			! ori endif
1595		1449226	qg = epses/(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	✓✗	1449226	IF (cvflag_ice) THEN
1608		1449226	tp(i, k) = max(0., (ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i)))
1609			ELSE
1610			tp(i, k) = (ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i))
1611			END IF
1612
1613		1449226	clw(i, k) = qnk(i) - qg
1614		1449226	clw(i, k) = max(0.0, clw(i,k))
1615			rg = qg/(1.-qnk(i))
1616			! ori tvp(i,k)=tp(i,k)(1.+rgepsi)
1617			! convect3: (qg utilise au lieu du vrai mixing ratio rg):
1618		1449226	tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
1619	✓✗	1449226	IF (cvflag_ice) THEN
1620	✓✓	1449226	IF (clw(i,k)<1.E-11) THEN
1621		263	tp(i, k) = tv(i, k)
1622		263	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	✓✗	1449226	IF (cvflag_ice) THEN
1630			!CR:attention boucle en klon dans Icefrac
1631			! Call Icefrac(t,clw,qi,nl,nloc)
1632	✓✓	1449226	IF (t(i,k)>263.15) THEN
1633		322715	qi(i, k) = 0.
1634			ELSE
1635	✓✓	1126511	IF (t(i,k)<243.15) THEN
1636		951380	qi(i, k) = clw(i, k)
1637			ELSE
1638		175131	fracg = (263.15-t(i,k))/20
1639		175131	qi(i, k) = clw(i, k)*fracg
1640			END IF
1641			END IF
1642			!CR: fin test
1643	✓✓	1449226	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=epsesi/(p(i,k)-esi(1.-eps))
1659			! snew=cpd(1.-qnk(i))+clqnk(i)+alvalvqsat_new/ &
1660			! (rrvtbistbis)
1661			! snew=1./snew
1662			! print*,esi,qsat_new,snew,'esi,qsat,snew'
1663			! tp(i,k)=tg+(alfqi(i,k)+alvqg(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		1126511	alf = lf0 + clmci*(t(i,k)-273.15)
1670		1126511	als = alf + alv
1671		1126511	tg = tp(i, k)
1672		1126511	tp(i, k) = t(i, k)
1673	✓✓	4506044	DO j = 1, 3
1674		3379533	esi = exp(23.33086-(6111.72784/tp(i,k))+0.15215*log(tp(i,k)))
1675		3379533	qsat_new = epsesi/(p(i,k)-esi(1.-eps))
1676			snew = cpd(1.-qnk(i)) + clqnk(i) + alvalsqsat_new/ &
1677		3379533	(rrvtp(i,k)tp(i,k))
1678		3379533	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))+clqnk(i))*(tg-tp(i,k)) + &
1682		4506044	alv(qg-qsat_new)+alfqi(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		1126511	clw(i, k) = qnk(i) - qsat_new
1689		1126511	clw(i, k) = max(0.0, clw(i,k))
1690		1126511	tvp(i, k) = max(0., tp(i,k)*(1.+qsat_new/eps-qnk(i)))
1691			! print*,tvp(i,k),'tvp'
1692			END IF
1693	✗✓	1449226	IF (clw(i,k)<1.E-11) THEN
1694			tp(i, k) = tv(i, k)
1695			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	✗✓	144	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	✓✓	4032	DO k = 1, nl
1731	✓✓	1866114	DO i = 1, ncum
1732	✓✓	1865970	IF(k>=icb(i)) THEN
1733			!! IF (k>=(nk(i)+1)) THEN
1734			!>jyg
1735		1556118	tca = tp(i, k) - t0
1736	✓✓	1556118	IF (tca>=0.0) THEN
1737		308269	elacrit = elcrit
1738			ELSE
1739		1247849	elacrit = elcrit*(1.0-tca/tlcrit)
1740			END IF
1741		1556118	elacrit = max(elacrit, 0.0)
1742		1556118	ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8)
1743		1556118	ep(i, k) = max(ep(i,k), 0.0)
1744		1556118	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	✗✓	144	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	✓✓	69110	DO i = 1, ncum ! convect3
1780		69110	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	✓✓	4032	DO k = 1, nl
1793	✓✓	1866114	DO i = 1, ncum
1794		1865970	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	✓✓	4032	DO k = 1, nl
1803	✓✓	1866114	DO i = 1, ncum
1804	✓✓✓✓	1865970	IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i,k)>=pbase(i))) THEN
1805		117979	buoy(i, k) = buoybase(i)
1806			END IF
1807			END DO
1808			END DO
1809	✓✓	69110	DO i = 1, ncum
1810			! buoy(icb(i),k)=buoybase(i)
1811		69110	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	✓✓	69110	DO i = 1, ncum
1824		68966	inb(i) = nl - 1
1825		69110	iposit(i) = nl
1826			END DO
1827
1828
1829			! -- iposit(i) = first level, above icb, with positive buoyancy
1830	✓✓	3888	DO k = 1, nl - 1
1831	✓✓	1797004	DO i = 1, ncum
1832	✓✓✓✓	1796860	IF (k>=icb(i) .AND. buoy(i,k)>0.) THEN
1833		501077	iposit(i) = min(iposit(i), k)
1834			END IF
1835			END DO
1836			END DO
1837
1838	✓✓	69110	DO i = 1, ncum
1839	✓✓	69110	IF (iposit(i)==nl) THEN
1840		4617	iposit(i) = icb(i)
1841			END IF
1842			END DO
1843
1844	✓✓	3888	DO k = 1, nl - 1
1845	✓✓	1797004	DO i = 1, ncum
1846	✓✓✓✓	1796860	IF ((k>=iposit(i)) .AND. (buoy(i,k)<dtovsh)) THEN
1847		935323	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	✗✓	144	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) + rrdbuoy(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	✓✓	69110	DO i = 1,ncum
1974	✗✓	69110	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	✓✓	4032	DO k = 1, nl
1982	✓✓	1866114	DO i = 1, ncum
1983		1865970	hp(i, k) = h(i, k)
1984			END DO
1985			END DO
1986
1987			!jyg : cvflag_ice test outside the loops (07042015)
1988			!
1989	✓✗	144	IF (cvflag_ice) THEN
1990			!
1991	✓✗	144	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	✓✓	3888	DO k = minorig + 1, nl
2002	✓✓	1797004	DO i = 1, ncum
2003	✓✓✓✓	1796860	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		586707	frac_s(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15)
2007		586707	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	✓✓	3888	DO k = minorig + 1, nl
2013	✓✓	1797004	DO i = 1, ncum
2014	✓✓✓✓	1796860	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		586707	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		144	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.alphabuoy(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.alphabuoy(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)=betasig(i)+2.alphabuoy(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)=betasig(i)-2.alphabuoy(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) + rrdbuoy(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) = betasig(i, k) + alphadtmin(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)facsqrt(cape(i)) + beta*w0(i, k)
2206			amu = 0.5(sig(i,k)+sigold(i,k))w
2207			m(i, k) = amu0.007p(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.5m(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+rrdbuoy(i-1)deltap/p(i-1)
2258			!! dcape=rrdbuoy(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)=betasig(i)+alphadtmin*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)facsqrt(cape)+beta*w0(i)
2273			!! amu=0.5(sig(i)+sigold)w
2274			!! m(i)=amu0.007p(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.5m(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)/(rrvt(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)/(rrvt(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)-cwatbf2)
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			IF (i>=icb(il) .AND. i<=inb(il)) THEN
2525			lwork(il) = (nent(il,i)/=0)
2526			qp = qnk(il) - ep(il, i)*clw(il, i)
2527
2528			IF (cvflag_ice) THEN
2529
2530			anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)lf(il,i)) &
2531			(qp-rs(il,i)) + (cpv-cpd)t(il, i)(qp-rr(il,i))
2532			denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)lf(il,i)) &
2533			(rr(il,i)-qp) + (cpd-cpv)t(il, i)(rr(il,i)-qp)
2534			ELSE
2535
2536			anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
2537			(cpv-cpd)t(il, i)(qp-rr(il,i))
2538			denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
2539			(cpd-cpv)t(il, i)(rr(il,i)-qp)
2540			END IF
2541
2542			IF (abs(denom)<0.01) denom = 0.01
2543			scrit(il) = anum/denom
2544			alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
2545			IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
2546			smax(il) = 0.0
2547			asij(il) = 0.0
2548			END IF
2549			END DO
2550
2551			DO j = nl, minorig, -1
2552
2553			num2 = 0
2554			DO il = 1, ncum
2555			IF (i>=icb(il) .AND. i<=inb(il) .AND. &
2556			j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
2557			lwork(il)) num2 = num2 + 1
2558			END DO
2559			IF (num2<=0) GO TO 175
2560
2561			DO il = 1, ncum
2562			IF (i>=icb(il) .AND. i<=inb(il) .AND. &
2563			j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
2564			lwork(il)) THEN
2565
2566			IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
2567			wgh = 1.0
2568			IF (j>i) THEN
2569			sjmax = max(sij(il,i,j+1), smax(il))
2570			sjmax = amin1(sjmax, scrit(il))
2571			smax(il) = max(sij(il,i,j), smax(il))
2572			sjmin = max(sij(il,i,j-1), smax(il))
2573			sjmin = amin1(sjmin, scrit(il))
2574			IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
2575			smid = amin1(sij(il,i,j), scrit(il))
2576			ELSE
2577			sjmax = max(sij(il,i,j+1), scrit(il))
2578			smid = max(sij(il,i,j), scrit(il))
2579			sjmin = 0.0
2580			IF (j>1) sjmin = sij(il, i, j-1)
2581			sjmin = max(sjmin, scrit(il))
2582			END IF
2583			delp = abs(sjmax-smid)
2584			delm = abs(sjmin-smid)
2585			asij(il) = asij(il) + wgh*(delp+delm)
2586			ment(il, i, j) = ment(il, i, j)(delp+delm)wgh
2587			END IF
2588			END IF
2589			END DO
2590
2591			175 END DO
2592
2593			DO il = 1, ncum
2594			IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
2595			asij(il) = max(1.0E-16, asij(il))
2596			asij(il) = 1.0/asij(il)
2597			asum(il, i) = 0.0
2598			bsum(il, i) = 0.0
2599			csum(il, i) = 0.0
2600			END IF
2601			END DO
2602
2603			DO j = minorig, nl
2604			DO il = 1, ncum
2605			IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2606			j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2607			ment(il, i, j) = ment(il, i, j)*asij(il)
2608			END IF
2609			END DO
2610			END DO
2611
2612			DO j = minorig, nl
2613			DO il = 1, ncum
2614			IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2615			j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2616			asum(il, i) = asum(il, i) + ment(il, i, j)
2617			ment(il, i, j) = ment(il, i, j)*sig(il, j)
2618			bsum(il, i) = bsum(il, i) + ment(il, i, j)
2619			END IF
2620			END DO
2621			END DO
2622
2623			DO il = 1, ncum
2624			IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
2625			bsum(il, i) = max(bsum(il,i), 1.0E-16)
2626			bsum(il, i) = 1.0/bsum(il, i)
2627			END IF
2628			END DO
2629
2630			DO j = minorig, nl
2631			DO il = 1, ncum
2632			IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2633			j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2634			ment(il, i, j) = ment(il, i, j)asum(il, i)bsum(il, i)
2635			END IF
2636			END DO
2637			END DO
2638
2639			DO j = minorig, nl
2640			DO il = 1, ncum
2641			IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2642			j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2643			csum(il, i) = csum(il, i) + ment(il, i, j)
2644			END IF
2645			END DO
2646			END DO
2647
2648			DO il = 1, ncum
2649			IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2650			csum(il,i)<m(il,i)) THEN
2651			nent(il, i) = 0
2652			ment(il, i, i) = m(il, i)
2653			qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
2654			uent(il, i, i) = unk(il)
2655			vent(il, i, i) = vnk(il)
2656			elij(il, i, i) = clw(il, i)
2657			! MAF sij(il,i,i)=1.0
2658			sij(il, i, i) = 0.0
2659			END IF
2660			END DO ! il
2661
2662			!AC! do j=1,ntra
2663			!AC! do il=1,ncum
2664			!AC! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
2665			!AC! : .and. csum(il,i).lt.m(il,i) ) then
2666			!AC! traent(il,i,i,j)=tra(il,nk(il),j)
2667			!AC! endif
2668			!AC! enddo
2669			!AC! enddo
2670			789 END DO
2671
2672			! MAF: renormalisation de MENT
2673			CALL zilch(zm, nloc*na)
2674			DO jm = 1, nl
2675			DO im = 1, nl
2676			DO il = 1, ncum
2677			zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
2678			END DO
2679			END DO
2680			END DO
2681
2682			DO jm = 1, nl
2683			DO im = 1, nl
2684			DO il = 1, ncum
2685			IF (zm(il,im)/=0.) THEN
2686			ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
2687			END IF
2688			END DO
2689			END DO
2690			END DO
2691
2692			DO jm = 1, nl
2693			DO im = 1, nl
2694			DO il = 1, ncum
2695			qents(il, im, jm) = qent(il, im, jm)
2696			ments(il, im, jm) = ment(il, im, jm)
2697			END DO
2698			END DO
2699			END DO
2700
2701			RETURN
2702			END SUBROUTINE cv3_mixing
2703
2704		7605067	SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, &
2705			t, rr, rs, gz, u, v, tra, p, ph, &
2706			th, tv, lv, lf, cpn, ep, sigp, clw, frac_s, qpreca, frac_a, qta , & !!jygprl
2707		144	m, ment, elij, delt, plcl, coef_clos, &
2708		144	mp, rp, up, vp, trap, wt, water, evap, fondue, ice, &
2709		144	faci, b, sigd, &
2710			wdtrainA, wdtrainS, wdtrainM) ! RomP
2711			USE print_control_mod, ONLY: prt_level, lunout
2712			IMPLICIT NONE
2713
2714
2715			include "cvthermo.h"
2716			include "cv3param.h"
2717			include "cvflag.h"
2718			include "nuage.h"
2719
2720			!inputs:
2721			INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
2722			INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
2723			REAL, INTENT(IN) :: delt
2724			REAL, DIMENSION (nloc), INTENT (IN) :: plcl
2725			REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
2726			REAL, DIMENSION (nloc, na), INTENT (IN) :: gz
2727			REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
2728			REAL, DIMENSION (nloc, nd, ntra), INTENT(IN) :: tra
2729			REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
2730			REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
2731			REAL, DIMENSION (nloc, na), INTENT (IN) :: ep, sigp, clw !adiab ascent shedding
2732			REAL, DIMENSION (nloc, na), INTENT (IN) :: frac_s !ice fraction in adiab ascent shedding !!jygprl
2733			REAL, DIMENSION (nloc, na), INTENT (IN) :: qpreca !adiab ascent precip !!jygprl
2734			REAL, DIMENSION (nloc, na), INTENT (IN) :: frac_a !ice fraction in adiab ascent precip !!jygprl
2735			REAL, DIMENSION (nloc, na), INTENT (IN) :: qta !adiab ascent specific total water !!jygprl
2736			REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tv, lv, cpn
2737			REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
2738			REAL, DIMENSION (nloc, na), INTENT (IN) :: m
2739			REAL, DIMENSION (nloc, na, na), INTENT (IN) :: ment, elij
2740			REAL, DIMENSION (nloc), INTENT (IN) :: coef_clos
2741
2742			!input/output
2743			INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag(nloc)
2744
2745			!outputs:
2746			REAL, DIMENSION (nloc, na), INTENT (OUT) :: mp, rp, up, vp
2747			REAL, DIMENSION (nloc, na), INTENT (OUT) :: water, evap, wt
2748			REAL, DIMENSION (nloc, na), INTENT (OUT) :: ice, fondue
2749			REAL, DIMENSION (nloc, na), INTENT (OUT) :: faci ! ice fraction in precipitation
2750			REAL, DIMENSION (nloc, na, ntra), INTENT (OUT) :: trap
2751			REAL, DIMENSION (nloc, na), INTENT (OUT) :: b
2752			REAL, DIMENSION (nloc), INTENT (OUT) :: sigd
2753			! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
2754			! de l ascendance adiabatique et des flux melanges Pa et Pm.
2755			! Distinction des wdtrain
2756			! Pa = wdtrainA Pm = wdtrainM
2757			REAL, DIMENSION (nloc, na), INTENT (OUT) :: wdtrainA, wdtrainS, wdtrainM
2758
2759			!local variables
2760			INTEGER i, j, k, il, num1, ndp1
2761			REAL smallestreal
2762			REAL tinv, delti, coef
2763			REAL awat, afac, afac1, afac2, bfac
2764			REAL pr1, pr2, sigt, b6, c6, d6, e6, f6, revap, delth
2765			REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
2766			REAL ampmax, thaw
2767		288	REAL tevap(nloc)
2768		288	REAL, DIMENSION (nloc, na) :: lvcp, lfcp
2769			REAL, DIMENSION (nloc, na) :: h, hm
2770		288	REAL, DIMENSION (nloc, na) :: ma
2771		288	REAL, DIMENSION (nloc, na) :: frac ! ice fraction in precipitation source
2772		288	REAL, DIMENSION (nloc, na) :: fraci ! provisionnal ice fraction in precipitation
2773		288	REAL, DIMENSION (nloc, na) :: prec
2774		288	REAL wdtrain(nloc)
2775		288	LOGICAL lwork(nloc), mplus(nloc)
2776
2777
2778			! ------------------------------------------------------
2779	✗✓	144	IF (prt_level .GE. 10) print *,' ->cv3_unsat, iflag(1) ', iflag(1)
2780
2781			smallestreal=tiny(smallestreal)
2782
2783			! =============================
2784			! --- INITIALIZE OUTPUT ARRAYS
2785			! =============================
2786			! (loops up to nl+1)
2787	✓✓✓✓	5588064	mp(:,:) = 0.
2788	✓✓✓✓	5588064	rp(:,:) = 0.
2789	✓✓✓✓	5588064	up(:,:) = 0.
2790	✓✓✓✓	5588064	vp(:,:) = 0.
2791	✓✓✓✓	5588064	water(:,:) = 0.
2792	✓✓✓✓	5588064	evap(:,:) = 0.
2793	✓✓✓✓	5588064	wt(:,:) = 0.
2794	✓✓✓✓	5588064	ice(:,:) = 0.
2795	✓✓✓✓	5588064	fondue(:,:) = 0.
2796	✓✓✓✓	5588064	faci(:,:) = 0.
2797	✓✓✓✓	5588064	b(:,:) = 0.
2798	✓✓	143280	sigd(:) = 0.
2799			!! RomP >>>
2800	✓✓✓✓	5588064	wdtrainA(:,:) = 0.
2801	✓✓✓✓	5588064	wdtrainS(:,:) = 0.
2802	✓✓✓✓	5588064	wdtrainM(:,:) = 0.
2803			!! RomP <<<
2804
2805	✓✓	4176	DO i = 1, nlp
2806	✓✓	1935224	DO il = 1, ncum
2807		1931048	rp(il, i) = rr(il, i)
2808		1931048	up(il, i) = u(il, i)
2809		1931048	vp(il, i) = v(il, i)
2810		1935080	wt(il, i) = 0.001
2811			END DO
2812			END DO
2813
2814			! *** Set the fractionnal area sigd of precipitating downdraughts
2815	✓✓	69110	DO il = 1, ncum
2816		69110	sigd(il) = sigdz*coef_clos(il)
2817			END DO
2818
2819			! =====================================================================
2820			! --- INITIALIZE VARIOUS ARRAYS AND PARAMETERS USED IN THE COMPUTATIONS
2821			! =====================================================================
2822			! (loops up to nl+1)
2823
2824		144	delti = 1./delt
2825			tinv = 1./3.
2826
2827	✓✓	4176	DO i = 1, nlp
2828	✓✓	1935224	DO il = 1, ncum
2829		1931048	frac(il, i) = 0.0
2830		1931048	fraci(il, i) = 0.0
2831		1931048	prec(il, i) = 0.0
2832		1931048	lvcp(il, i) = lv(il, i)/cpn(il, i)
2833		1935080	lfcp(il, i) = lf(il, i)/cpn(il, i)
2834			END DO
2835			END DO
2836
2837			!AC! do k=1,ntra
2838			!AC! do i=1,nd
2839			!AC! do il=1,ncum
2840			!AC! trap(il,i,k)=tra(il,i,k)
2841			!AC! enddo
2842			!AC! enddo
2843			!AC! enddo
2844
2845			! * check whether ep(inb)=0, if so, skip precipitating *
2846			! * downdraft calculation *
2847
2848
2849	✓✓	69110	DO il = 1, ncum
2850			!! lwork(il)=.TRUE.
2851			!! if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
2852			!jyg<
2853			!! lwork(il) = ep(il, inb(il)) >= 0.0001
2854	✓✓✓✓	115057	lwork(il) = ep(il, inb(il)) >= 0.0001 .AND. iflag(il) <= 2
2855			END DO
2856
2857			!
2858			! Get adiabatic ascent mass flux
2859			!
2860			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2861	✗✓	144	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
2862			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2863			!!! Warning : this option leads to water conservation violation
2864			!!! Expert only
2865			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2866			DO il = 1, ncum
2867			ma(il, nlp) = 0.
2868			ma(il, 1) = 0.
2869			END DO
2870
2871			DO i = nl, 2, -1
2872			DO il = 1, ncum
2873			ma(il, i) = ma(il, i+1)*(1.-qta(il,i))/(1.-qta(il,i-1)) + m(il, i)
2874			END DO
2875			END DO
2876			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2877			ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
2878			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2879	✓✓	69110	DO il = 1, ncum
2880		68966	ma(il, nlp) = 0.
2881		69110	ma(il, 1) = 0.
2882			END DO
2883
2884	✓✓	3888	DO i = nl, 2, -1
2885	✓✓	1797004	DO il = 1, ncum
2886		1796860	ma(il, i) = ma(il, i+1) + m(il, i)
2887			END DO
2888			END DO
2889
2890			ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
2891			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2892
2893			! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2894			!
2895			! * begin downdraft loop *
2896			!
2897			! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2898
2899	✓✓	4176	DO i = nl + 1, 1, -1
2900
2901			num1 = 0
2902	✓✓	1935080	DO il = 1, ncum
2903	✓✓✓✓	1935080	IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
2904			END DO
2905	✓✓	4032	IF (num1<=0) GO TO 400
2906
2907		3027	CALL zilch(wdtrain, ncum)
2908
2909
2910			! * integrate liquid water equation to find condensed water *
2911			! * and condensed water flux *
2912			!
2913			!
2914			! * calculate detrained precipitation *
2915
2916
2917	✓✓	1452914	DO il = 1, ncum
2918	✓✓✓✓	1452914	IF (i<=inb(il) .AND. lwork(il)) THEN
2919		386747	wdtrain(il) = gravep(il, i)m(il, i)*clw(il, i)
2920		386747	wdtrainS(il, i) = wdtrain(il)/grav ! Ps jyg
2921			!! wdtrainA(il, i) = wdtrain(il)/grav ! Ps RomP
2922			END IF
2923			END DO
2924
2925	✓✓	3027	IF (i>1) THEN
2926	✓✓	33186	DO j = 1, i - 1
2927	✓✓	14549667	DO il = 1, ncum
2928	✓✓✓✓	14546784	IF (i<=inb(il) .AND. lwork(il)) THEN
2929		3192455	awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
2930		3192455	awat = max(awat, 0.0)
2931		3192455	wdtrain(il) = wdtrain(il) + gravawatment(il, j, i)
2932		3192455	wdtrainM(il, i) = wdtrain(il)/grav - wdtrainS(il, i) ! Pm jyg
2933			!! wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i) ! Pm RomP
2934			END IF
2935			END DO
2936			END DO
2937			END IF
2938
2939	✗✓	3027	IF (cvflag_prec_eject) THEN
2940			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2941			IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
2942			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2943			!!! Warning : this option leads to water conservation violation
2944			!!! Expert only
2945			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2946			IF ( i > 1) THEN
2947			DO il = 1, ncum
2948			IF (i<=inb(il) .AND. lwork(il)) THEN
2949			wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i))/(1. - qta(il, i-1)) ! Pa jygprl
2950			wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i)
2951			END IF
2952			END DO
2953			ENDIF ! ( i > 1)
2954			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2955			ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
2956			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2957			IF ( i > 1) THEN
2958			DO il = 1, ncum
2959			IF (i<=inb(il) .AND. lwork(il)) THEN
2960			wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i)) ! Pa jygprl
2961			wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i)
2962			END IF
2963			END DO
2964			ENDIF ! ( i > 1)
2965
2966			ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
2967			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2968			ENDIF ! (cvflag_prec_eject)
2969
2970
2971			! * find rain water and evaporation using provisional *
2972			! * estimates of rp(i)and rp(i-1) *
2973
2974
2975	✓✗	3027	IF (cvflag_ice) THEN !!jygprl
2976	✗✓	3027	IF (cvflag_prec_eject) THEN
2977			DO il = 1, ncum !!jygprl
2978			IF (i<=inb(il) .AND. lwork(il)) THEN !!jygprl
2979			frac(il, i) = (frac_a(il,i)wdtrainA(il,i)+frac_s(il,i)(wdtrainS(il,i)+wdtrainM(il,i))) / & !!jygprl
2980			max(wdtrainA(il,i)+wdtrainS(il,i)+wdtrainM(il,i),smallestreal) !!jygprl
2981			fraci(il, i) = frac(il, i) !!jygprl
2982			END IF !!jygprl
2983			END DO !!jygprl
2984			ELSE ! (cvflag_prec_eject)
2985	✓✓	1452914	DO il = 1, ncum !!jygprl
2986	✓✓✓✓	1452914	IF (i<=inb(il) .AND. lwork(il)) THEN !!jygprl
2987			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2988	✗✓	386747	IF (keepbug_ice_frac) THEN
2989			frac(il, i) = frac_s(il, i)
2990			! Ice fraction computed again here as a function of the temperature seen by unsaturated downdraughts
2991			! (i.e. the cold pool temperature) for compatibility with earlier versions.
2992			fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15)
2993			fraci(il, i) = min(max(fraci(il,i),0.0), 1.0)
2994			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2995			ELSE ! (keepbug_ice_frac)
2996			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2997		386747	frac(il, i) = frac_s(il, i)
2998		386747	fraci(il, i) = frac(il, i) !!jygprl
2999			ENDIF ! (keepbug_ice_frac)
3000			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3001			END IF !!jygprl
3002			END DO !!jygprl
3003			ENDIF ! (cvflag_prec_eject)
3004			END IF !!jygprl
3005
3006
3007	✓✓	1452914	DO il = 1, ncum
3008	✓✓✓✓	1452914	IF (i<=inb(il) .AND. lwork(il)) THEN
3009
3010		386747	wt(il, i) = 45.0
3011
3012	✓✓	386747	IF (i<inb(il)) THEN
3013			rp(il, i) = rp(il, i+1) + &
3014		363728	(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il, i)
3015		363728	rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
3016			END IF
3017		386747	rp(il, i) = max(rp(il,i), 0.0)
3018		386747	rp(il, i) = amin1(rp(il,i), rs(il,i))
3019		386747	rp(il, inb(il)) = rr(il, inb(il))
3020
3021	✓✓	386747	IF (i==1) THEN
3022		23019	afac = p(il, 1)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
3023	✓✗	23019	IF (cvflag_ice) THEN
3024		23019	afac1 = p(il, i)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
3025			END IF
3026			ELSE
3027		363728	rp(il, i-1) = rp(il, i) + (cpd*(t(il,i)-t(il,i-1))+gz(il,i)-gz(il,i-1))/lv(il, i)
3028		363728	rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
3029		363728	rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
3030		363728	rp(il, i-1) = max(rp(il,i-1), 0.0)
3031		363728	afac1 = p(il, i)(rs(il,i)-rp(il,i))/(1.0E4+2000.0p(il,i)*rs(il,i))
3032		363728	afac2 = p(il, i-1)(rs(il,i-1)-rp(il,i-1))/(1.0E4+2000.0p(il,i-1)*rs(il,i-1))
3033		363728	afac = 0.5*(afac1+afac2)
3034			END IF
3035	✓✓	386747	IF (i==inb(il)) afac = 0.0
3036		386747	afac = max(afac, 0.0)
3037		386747	bfac = 1./(sigd(il)*wt(il,i))
3038
3039			!
3040	✗✓	386747	IF (prt_level >= 20) THEN
3041			Print*, 'cv3_unsat after provisional rp estimate: rp, afac, bfac ', &
3042			i, rp(1, i), afac,bfac
3043			ENDIF
3044			!
3045			!JYG1
3046			! cc sigt=1.0
3047			! cc if(i.ge.icb)sigt=sigp(i)
3048			! prise en compte de la variation progressive de sigt dans
3049			! les couches icb et icb-1:
3050			! pour plcl<ph(i+1), pr1=0 & pr2=1
3051			! pour plcl>ph(i), pr1=1 & pr2=0
3052			! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
3053			! sur le nuage, et pr2 est la proportion sous la base du
3054			! nuage.
3055		386747	pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
3056		386747	pr1 = max(0., min(1.,pr1))
3057		386747	pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
3058		386747	pr2 = max(0., min(1.,pr2))
3059		386747	sigt = sigp(il, i)*pr1 + pr2
3060			!JYG2
3061
3062			!JYG----
3063			! b6 = bfac100.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
3064			! c6 = water(il,i+1) + wdtrain(il)*bfac
3065			! c6 = prec(il,i+1) + wdtrain(il)*bfac
3066			! revap=0.5(-b6+sqrt(b6b6+4.*c6))
3067			! evap(il,i)=sigtafacrevap
3068			! water(il,i)=revap*revap
3069			! prec(il,i)=revap*revap
3070			!! print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ', &
3071			!! i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
3072			!!---end jyg---
3073
3074			! --------retour � la formulation originale d'Emanuel.
3075	✓✗	386747	IF (cvflag_ice) THEN
3076
3077			! b6=bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
3078			! c6=prec(il,i+1)+bfac*wdtrain(il) &
3079			! -50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
3080			! if(c6.gt.0.0)then
3081			! revap=0.5(-b6+sqrt(b6b6+4.*c6))
3082
3083			!JAM Attention: evap=sigt*E
3084			! Modification: evap devient l'�vaporation en milieu de couche
3085			! car n�cessaire dans cv3_yield
3086			! Du coup, il faut modifier pas mal d'�quations...
3087			! et l'expression de afac qui devient afac1
3088			! revap=sqrt((prec(i+1)+prec(i))/2)
3089
3090		386747	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac1
3091		386747	c6 = prec(il, i+1) + 0.5bfacwdtrain(il)
3092			! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
3093			! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
3094			! print ,'b6,c6,b6b6+4.c6 ',b6,c6,b6b6+4.*c6
3095	✓✓	386747	IF (c6>b6*b6+1.E-20) THEN
3096		381250	revap = 2.c6/(b6+sqrt(b6b6+4.*c6))
3097			ELSE
3098		5497	revap = (-b6+sqrt(b6b6+4.c6))/2.
3099			END IF
3100		386747	prec(il, i) = max(0., 2.revaprevap-prec(il,i+1))
3101			! print*,prec(il,i),'neige'
3102
3103			!JYG Dans sa formulation originale, Emanuel calcule l'evaporation par:
3104			! c evap(il,i)=sigtafacrevap
3105			! ce qui n'est pas correct. Dans cv_routines, la formulation a �t� modifiee.
3106			! Ici,l'evaporation evap est simplement calculee par l'equation de
3107			! conservation.
3108			! prec(il,i)=revap*revap
3109			! else
3110			!JYG---- Correction : si c6 <= 0, water(il,i)=0.
3111			! prec(il,i)=0.
3112			! endif
3113
3114			!JYG--- Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
3115			! moins [tt ce qui sort de la couche i]
3116			! print *, 'evap avec ice'
3117			evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(prec(il,i+1)-prec(il,i))) / &
3118		386747	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
3119			!
3120	✗✓	386747	IF (prt_level >= 20) THEN
3121			Print*, 'cv3_unsat after evap computation: wdtrain, sigd, wt, prec(i+1),prec(i) ', &
3122			i, wdtrain(1), sigd(1), wt(1,i), prec(1,i+1),prec(1,i)
3123			ENDIF
3124			!
3125
3126			!jyg<
3127		386747	d6 = prec(il,i)-prec(il,i+1)
3128
3129			!! d6 = bfacwdtrain(il) - 100.sigd(il)bfac(ph(il,i)-ph(il,i+1))*evap(il, i)
3130			!! e6 = bfac*wdtrain(il)
3131			!! f6 = -100.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i)
3132			!>jyg
3133			!CR:tmax_fonte_cv: T for which ice is totally melted (used to be 275.15)
3134		386747	thaw = (t(il,i)-273.15)/(tmax_fonte_cv-273.15)
3135		386747	thaw = min(max(thaw,0.0), 1.0)
3136			!jyg<
3137		386747	water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
3138		386747	ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
3139		386747	water(il, i) = min(prec(il,i), max(water(il,i), 0.))
3140		386747	ice(il, i) = min(prec(il,i), max(ice(il,i), 0.))
3141
3142			!! water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
3143			!! water(il, i) = max(water(il,i), 0.)
3144			!! ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
3145			!! ice(il, i) = max(ice(il,i), 0.)
3146			!>jyg
3147		386747	fondue(il, i) = ice(il, i)*thaw
3148		386747	water(il, i) = water(il, i) + fondue(il, i)
3149		386747	ice(il, i) = ice(il, i) - fondue(il, i)
3150
3151	✓✓	386747	IF (water(il,i)+ice(il,i)<1.E-30) THEN
3152		4681	faci(il, i) = 0.
3153			ELSE
3154		382066	faci(il, i) = ice(il, i)/(water(il,i)+ice(il,i))
3155			END IF
3156
3157			! water(il,i)=water(il,i+1)+(1.-fraci(il,i))e6+(1.-faci(il,i))f6
3158			! water(il,i)=max(water(il,i),0.)
3159			! ice(il,i)=ice(il,i+1)+fraci(il,i)e6+faci(il,i)f6
3160			! ice(il,i)=max(ice(il,i),0.)
3161			! fondue(il,i)=ice(il,i)*thaw
3162			! water(il,i)=water(il,i)+fondue(il,i)
3163			! ice(il,i)=ice(il,i)-fondue(il,i)
3164
3165			! if((water(il,i)+ice(il,i)).lt.1.e-30)then
3166			! faci(il,i)=0.
3167			! else
3168			! faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
3169			! endif
3170
3171			ELSE
3172			b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
3173			c6 = water(il, i+1) + bfac*wdtrain(il) - &
3174			50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i+1)
3175			IF (c6>0.0) THEN
3176			revap = 0.5(-b6+sqrt(b6b6+4.*c6))
3177			water(il, i) = revap*revap
3178			ELSE
3179			water(il, i) = 0.
3180			END IF
3181			! print *, 'evap sans ice'
3182			evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(water(il,i+1)-water(il,i)))/ &
3183			(sigd(il)(ph(il,i)-ph(il,i+1))100.)
3184
3185			END IF
3186			END IF !(i.le.inb(il) .and. lwork(il))
3187			END DO
3188			! ----------------------------------------------------------------
3189
3190			! cc
3191			! * calculate precipitating downdraft mass flux under *
3192			! * hydrostatic approximation *
3193
3194	✓✓	1452914	DO il = 1, ncum
3195	✓✓✓✓ ✓✓	1449887	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
3196
3197		363728	tevap(il) = max(0.0, evap(il,i))
3198		363728	delth = max(0.001, (th(il,i)-th(il,i-1)))
3199	✓✗	363728	IF (cvflag_ice) THEN
3200	✓✗	363728	IF (cvflag_grav) THEN
3201			mp(il, i) = 100.ginv(lvcp(il,i)sigd(il)tevap(il)* &
3202			(p(il,i-1)-p(il,i))/delth + &
3203			lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
3204			(p(il,i-1)-p(il,i))/delth + &
3205			lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
3206		363728	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
3207			ELSE
3208			mp(il, i) = 10.(lvcp(il,i)sigd(il)tevap(il) &
3209			(p(il,i-1)-p(il,i))/delth + &
3210			lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
3211			(p(il,i-1)-p(il,i))/delth + &
3212			lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
3213			(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
3214
3215			END IF
3216			ELSE
3217			IF (cvflag_grav) THEN
3218			mp(il, i) = 100.ginvlvcp(il, i)sigd(il)tevap(il)* &
3219			(p(il,i-1)-p(il,i))/delth
3220			ELSE
3221			mp(il, i) = 10.lvcp(il, i)sigd(il)tevap(il) &
3222			(p(il,i-1)-p(il,i))/delth
3223			END IF
3224
3225			END IF
3226
3227			END IF !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
3228	✗✓	1452914	IF (prt_level .GE. 20) THEN
3229			PRINT *,'cv3_unsat, mp hydrostatic ', i, mp(il,i)
3230			ENDIF
3231			END DO
3232			! ----------------------------------------------------------------
3233
3234			! * if hydrostatic assumption fails, *
3235			! * solve cubic difference equation for downdraft theta *
3236			! * and mass flux from two simultaneous differential eqns *
3237
3238	✓✓	1452914	DO il = 1, ncum
3239	✓✓✓✓ ✓✓	1452914	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
3240
3241			amfac = sigd(il)sigd(il)70.0ph(il, i)(p(il,i-1)-p(il,i))* &
3242		363728	(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
3243		363728	amp2 = abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
3244
3245	✓✓	363728	IF (amp2>(0.1*amfac)) THEN
3246		97873	xf = 100.0sigd(il)sigd(il)sigd(il)(ph(il,i)-ph(il,i+1))
3247			tf = b(il, i) - 5.0(th(il,i)-th(il,i-1))t(il, i) / &
3248		97873	(lvcp(il,i)sigd(il)th(il,i))
3249		97873	af = xftf + mp(il, i+1)mp(il, i+1)*tinv
3250
3251	✓✗	97873	IF (cvflag_ice) THEN
3252			bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
3253			50.(p(il,i-1)-p(il,i))xf(tevap(il)(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
3254		97873	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i)/(ph(il,i)-ph(il,i+1)))
3255			ELSE
3256
3257			bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
3258			50.(p(il,i-1)-p(il,i))xf*tevap(il)
3259			END IF
3260
3261			fac2 = 1.0
3262	✓✓	97873	IF (bf<0.0) fac2 = -1.0
3263		97873	bf = abs(bf)
3264		97873	ur = 0.25bfbf - afafaftinvtinv*tinv
3265	✓✓	97873	IF (ur>=0.0) THEN
3266		63062	sru = sqrt(ur)
3267			fac = 1.0
3268	✓✓	63062	IF ((0.5*bf-sru)<0.0) fac = -1.0
3269			mp(il, i) = mp(il, i+1)tinv + (0.5bf+sru)**tinv + &
3270		63062	fac(abs(0.5bf-sru))**tinv
3271			ELSE
3272		34811	d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
3273	✗✓	34811	IF (fac2<0.0) d = 3.14159 - d
3274		34811	mp(il, i) = mp(il, i+1)tinv + 2.sqrt(aftinv)cos(d*tinv)
3275			END IF
3276		97873	mp(il, i) = max(0.0, mp(il,i))
3277	✗✓	97873	IF (prt_level .GE. 20) THEN
3278			PRINT *,'cv3_unsat, mp cubic ', i, mp(il,i)
3279			ENDIF
3280
3281	✓✗	97873	IF (cvflag_ice) THEN
3282	✓✗	97873	IF (cvflag_grav) THEN
3283			!JYG : il y a vraisemblablement une erreur dans la ligne 2 suivante:
3284			! il faut diviser par (mp(il,i)sigd(il)grav) et non par (mp(il,i)+sigd(il)*0.1).
3285			! Et il faut bien revoir les facteurs 100.
3286			b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i)) &
3287			(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
3288			(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
3289			(ph(il,i)-ph(il,i+1))) / &
3290			(mp(il,i)+sigd(il)*0.1) - &
3291			10.0(th(il,i)-th(il,i-1))t(il, i) / &
3292		97873	(lvcp(il,i)sigd(il)th(il,i))
3293			ELSE
3294			b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))&
3295			(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
3296			(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
3297			(ph(il,i)-ph(il,i+1))) / &
3298			(mp(il,i)+sigd(il)*0.1) - &
3299			10.0(th(il,i)-th(il,i-1))t(il, i) / &
3300			(lvcp(il,i)sigd(il)th(il,i))
3301			END IF
3302			ELSE
3303			IF (cvflag_grav) THEN
3304			b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
3305			(mp(il,i)+sigd(il)*0.1) - &
3306			10.0(th(il,i)-th(il,i-1))t(il, i) / &
3307			(lvcp(il,i)sigd(il)th(il,i))
3308			ELSE
3309			b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
3310			(mp(il,i)+sigd(il)*0.1) - &
3311			10.0(th(il,i)-th(il,i-1))t(il, i) / &
3312			(lvcp(il,i)sigd(il)th(il,i))
3313			END IF
3314			END IF
3315		97873	b(il, i-1) = max(b(il,i-1), 0.0)
3316
3317			END IF !(amp2.gt.(0.1*amfac))
3318
3319			!jyg< This part shifted 10 lines farther
3320			!!! * limit magnitude of mp(i) to meet cfl condition *
3321			!!
3322			!! ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
3323			!! amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
3324			!! ampmax = min(ampmax, amp2)
3325			!! mp(il, i) = min(mp(il,i), ampmax)
3326			!>jyg
3327
3328			! * force mp to decrease linearly to zero *
3329			! * between cloud base and the surface *
3330
3331
3332			! c if(p(il,i).gt.p(il,icb(il)))then
3333			! c mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
3334			! c endif
3335	✓✓	363728	IF (ph(il,i)>0.9*plcl(il)) THEN
3336		152436	mp(il, i) = mp(il, i)(ph(il,1)-ph(il,i))/(ph(il,1)-0.9plcl(il))
3337			END IF
3338
3339			!jyg< Shifted part
3340			! * limit magnitude of mp(i) to meet cfl condition *
3341
3342		363728	ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
3343		363728	amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
3344		363728	ampmax = min(ampmax, amp2)
3345		363728	mp(il, i) = min(mp(il,i), ampmax)
3346			!>jyg
3347
3348			END IF ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
3349			END DO
3350			! ----------------------------------------------------------------
3351			!
3352	✗✓	3027	IF (prt_level >= 20) THEN
3353			Print*, 'cv3_unsat after mp computation: mp, b(i), b(i-1) ', &
3354			i, mp(1, i), b(1,i), b(1,max(i-1,1))
3355			ENDIF
3356			!
3357
3358			! * find mixing ratio of precipitating downdraft *
3359
3360	✓✓	1452914	DO il = 1, ncum
3361	✓✓✓✓	1452914	IF (i<inb(il) .AND. lwork(il)) THEN
3362		363728	mplus(il) = mp(il, i) > mp(il, i+1)
3363			END IF ! (i.lt.inb(il) .and. lwork(il))
3364			END DO
3365
3366	✓✓	1452914	DO il = 1, ncum
3367	✓✓✓✓	1453919	IF (i<inb(il) .AND. lwork(il)) THEN
3368
3369		363728	rp(il, i) = rr(il, i)
3370
3371	✓✓	363728	IF (mplus(il)) THEN
3372
3373	✓✗	157029	IF (cvflag_grav) THEN
3374			rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
3375		157029	100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
3376			ELSE
3377			rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
3378			5.sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
3379			END IF
3380		157029	rp(il, i) = rp(il, i)/mp(il, i)
3381		157029	up(il, i) = up(il, i+1)mp(il, i+1) + u(il, i)(mp(il,i)-mp(il,i+1))
3382		157029	up(il, i) = up(il, i)/mp(il, i)
3383		157029	vp(il, i) = vp(il, i+1)mp(il, i+1) + v(il, i)(mp(il,i)-mp(il,i+1))
3384		157029	vp(il, i) = vp(il, i)/mp(il, i)
3385
3386			ELSE ! if (mplus(il))
3387
3388	✓✓	206699	IF (mp(il,i+1)>1.0E-16) THEN
3389	✓✗	170418	IF (cvflag_grav) THEN
3390			rp(il, i) = rp(il,i+1) + 100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1)) * &
3391		170418	(evap(il,i+1)+evap(il,i))/mp(il,i+1)
3392			ELSE
3393			rp(il, i) = rp(il,i+1) + 5.sigd(il)(ph(il,i)-ph(il,i+1)) * &
3394			(evap(il,i+1)+evap(il,i))/mp(il, i+1)
3395			END IF
3396		170418	up(il, i) = up(il, i+1)
3397		170418	vp(il, i) = vp(il, i+1)
3398			END IF ! (mp(il,i+1).gt.1.0e-16)
3399			END IF ! (mplus(il)) else if (.not.mplus(il))
3400
3401		363728	rp(il, i) = amin1(rp(il,i), rs(il,i))
3402		363728	rp(il, i) = max(rp(il,i), 0.0)
3403
3404			END IF ! (i.lt.inb(il) .and. lwork(il))
3405			END DO
3406			! ----------------------------------------------------------------
3407
3408			! * find tracer concentrations in precipitating downdraft *
3409
3410			!AC! do j=1,ntra
3411			!AC! do il = 1,ncum
3412			!AC! if (i.lt.inb(il) .and. lwork(il)) then
3413			!AC!c
3414			!AC! if(mplus(il))then
3415			!AC! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
3416			!AC! : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
3417			!AC! trap(il,i,j)=trap(il,i,j)/mp(il,i)
3418			!AC! else ! if (mplus(il))
3419			!AC! if(mp(il,i+1).gt.1.0e-16)then
3420			!AC! trap(il,i,j)=trap(il,i+1,j)
3421			!AC! endif
3422			!AC! endif ! (mplus(il)) else if (.not.mplus(il))
3423			!AC!c
3424			!AC! endif ! (i.lt.inb(il) .and. lwork(il))
3425			!AC! enddo
3426			!AC! end do
3427
3428		144	400 END DO
3429			! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3430
3431			! * end of downdraft loop *
3432
3433			! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3434
3435
3436		144	RETURN
3437
3438			END SUBROUTINE cv3_unsat
3439
3440		19406870	SUBROUTINE cv3_yield(nloc, ncum, nd, na, ntra, ok_conserv_q, &
3441			icb, inb, delt, &
3442			t, rr, t_wake, rr_wake, s_wake, u, v, tra, &
3443		144	gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
3444			ep, clw, qpreca, m, tp, mp, rp, up, vp, trap, &
3445			wt, water, ice, evap, fondue, faci, b, sigd, &
3446		144	ment, qent, hent, iflag_mix, uent, vent, &
3447			nent, elij, traent, sig, &
3448			tv, tvp, wghti, &
3449		144	iflag, precip, Vprecip, Vprecipi, & ! jyg: Vprecipi
3450			ft, fr, fu, fv, ftra, & ! jyg
3451		144	cbmf, upwd, dnwd, dnwd0, ma, mip, &
3452			!! tls, tps, ! useless . jyg
3453			qcondc, wd, &
3454			ftd, fqd, qta, qtc, sigt, tau_cld_cv, coefw_cld_cv)
3455
3456			USE print_control_mod, ONLY: lunout, prt_level
3457			USE add_phys_tend_mod, only : fl_cor_ebil
3458
3459			IMPLICIT NONE
3460
3461			include "cvthermo.h"
3462			include "cv3param.h"
3463			include "cvflag.h"
3464			include "conema3.h"
3465
3466			!inputs:
3467			INTEGER, INTENT (IN) :: iflag_mix
3468			INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
3469			LOGICAL, INTENT (IN) :: ok_conserv_q
3470			INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
3471			REAL, INTENT (IN) :: delt
3472			REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, u, v
3473			REAL, DIMENSION (nloc, nd), INTENT (IN) :: t_wake, rr_wake
3474			REAL, DIMENSION (nloc), INTENT (IN) :: s_wake
3475			REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra
3476			REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
3477			REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
3478			REAL, DIMENSION (nloc, na), INTENT (IN) :: gz, h, hp
3479			REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tp
3480			REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, cpn, ep, clw
3481			REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
3482			REAL, DIMENSION (nloc, na), INTENT (IN) :: rp, up
3483			REAL, DIMENSION (nloc, na), INTENT (IN) :: vp
3484			REAL, DIMENSION (nloc, nd), INTENT (IN) :: wt
3485			REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: trap
3486			REAL, DIMENSION (nloc, na), INTENT (IN) :: water, evap, b
3487			REAL, DIMENSION (nloc, na), INTENT (IN) :: fondue, faci, ice
3488			REAL, DIMENSION (nloc, na, na), INTENT (IN) :: qent, uent
3489			REAL, DIMENSION (nloc, na, na), INTENT (IN) :: hent
3490			REAL, DIMENSION (nloc, na, na), INTENT (IN) :: vent, elij
3491			INTEGER, DIMENSION (nloc, nd), INTENT (IN) :: nent
3492			REAL, DIMENSION (nloc, na, na, ntra), INTENT (IN) :: traent
3493			REAL, DIMENSION (nloc, nd), INTENT (IN) :: tv, tvp, wghti
3494			REAL, DIMENSION (nloc, nd), INTENT (IN) :: qta
3495			REAL, DIMENSION (nloc, na),INTENT(IN) :: qpreca
3496			REAL, INTENT(IN) :: tau_cld_cv, coefw_cld_cv
3497			!
3498			!input/output:
3499			REAL, DIMENSION (nloc, na), INTENT (INOUT) :: m, mp
3500			REAL, DIMENSION (nloc, na, na), INTENT (INOUT) :: ment
3501			INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag
3502			REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: sig
3503			REAL, DIMENSION (nloc), INTENT (INOUT) :: sigd
3504			!
3505			!outputs:
3506			REAL, DIMENSION (nloc), INTENT (OUT) :: precip
3507			REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ft, fr, fu, fv
3508			REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ftd, fqd
3509			REAL, DIMENSION (nloc, nd, ntra), INTENT (OUT) :: ftra
3510			REAL, DIMENSION (nloc, nd), INTENT (OUT) :: upwd, dnwd, ma
3511			REAL, DIMENSION (nloc, nd), INTENT (OUT) :: dnwd0, mip
3512			REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecip
3513			REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecipi
3514			!! REAL tls(nloc, nd), tps(nloc, nd) ! useless . jyg
3515			REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qcondc ! cld
3516			REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qtc, sigt ! cld
3517			REAL, DIMENSION (nloc), INTENT (OUT) :: wd ! gust
3518			REAL, DIMENSION (nloc), INTENT (OUT) :: cbmf
3519			!
3520			!local variables:
3521			INTEGER :: i, k, il, n, j, num1
3522			REAL :: rat, delti
3523			REAL :: ax, bx, cx, dx, ex
3524			REAL :: cpinv, rdcp, dpinv
3525			REAL :: sigaq
3526		288	REAL, DIMENSION (nloc) :: awat
3527		288	REAL, DIMENSION (nloc, nd) :: lvcp, lfcp ! , mke ! unused . jyg
3528		288	REAL, DIMENSION (nloc) :: am, work, ad, amp1
3529			!! real up1(nloc), dn1(nloc)
3530		288	REAL, DIMENSION (nloc, nd, nd) :: up1, dn1
3531			!jyg<
3532		288	REAL, DIMENSION (nloc, nd) :: up_to, up_from
3533		288	REAL, DIMENSION (nloc, nd) :: dn_to, dn_from
3534			!>jyg
3535		288	REAL, DIMENSION (nloc) :: asum, bsum, csum, dsum
3536		288	REAL, DIMENSION (nloc) :: esum, fsum, gsum, hsum
3537			REAL, DIMENSION (nloc, nd) :: th_wake
3538		288	REAL, DIMENSION (nloc) :: alpha_qpos, alpha_qpos1
3539		288	REAL, DIMENSION (nloc, nd) :: qcond, nqcond, wa ! cld
3540		288	REAL, DIMENSION (nloc, nd) :: siga, sax, mac ! cld
3541		288	REAL, DIMENSION (nloc) :: sument
3542		288	REAL, DIMENSION (nloc, nd) :: sigment, qtment ! cld
3543			REAL sumdq !jyg
3544			!
3545			! -------------------------------------------------------------
3546
3547			! initialization:
3548
3549		144	delti = 1.0/delt
3550			! print*,'cv3_yield initialisation delt', delt
3551
3552	✓✓	69110	DO il = 1, ncum
3553		68966	precip(il) = 0.0
3554		69110	wd(il) = 0.0 ! gust
3555			END DO
3556
3557			! Fluxes are on a staggered grid : loops extend up to nl+1
3558	✓✓	4176	DO i = 1, nlp
3559	✓✓	1935224	DO il = 1, ncum
3560		1931048	Vprecip(il, i) = 0.0
3561		1931048	Vprecipi(il, i) = 0.0 ! jyg
3562		1931048	upwd(il, i) = 0.0
3563		1931048	dnwd(il, i) = 0.0
3564		1931048	dnwd0(il, i) = 0.0
3565		1935080	mip(il, i) = 0.0
3566			END DO
3567			END DO
3568	✓✓	4032	DO i = 1, nl
3569	✓✓	1866114	DO il = 1, ncum
3570		1862082	ft(il, i) = 0.0
3571		1862082	fr(il, i) = 0.0
3572		1862082	fu(il, i) = 0.0
3573		1862082	fv(il, i) = 0.0
3574		1862082	ftd(il, i) = 0.0
3575		1862082	fqd(il, i) = 0.0
3576		1862082	qcondc(il, i) = 0.0 ! cld
3577		1862082	qcond(il, i) = 0.0 ! cld
3578		1862082	qtc(il, i) = 0.0 ! cld
3579		1862082	qtment(il, i) = 0.0 ! cld
3580		1862082	sigment(il, i) = 0.0 ! cld
3581		1862082	sigt(il, i) = 0.0 ! cld
3582		1865970	nqcond(il, i) = 0.0 ! cld
3583			END DO
3584			END DO
3585			! print*,'cv3_yield initialisation 2'
3586			!AC! do j=1,ntra
3587			!AC! do i=1,nd
3588			!AC! do il=1,ncum
3589			!AC! ftra(il,i,j)=0.0
3590			!AC! enddo
3591			!AC! enddo
3592			!AC! enddo
3593			! print*,'cv3_yield initialisation 3'
3594	✓✓	4032	DO i = 1, nl
3595	✓✓	1866114	DO il = 1, ncum
3596		1862082	lvcp(il, i) = lv(il, i)/cpn(il, i)
3597		1865970	lfcp(il, i) = lf(il, i)/cpn(il, i)
3598			END DO
3599			END DO
3600
3601
3602
3603			! * calculate surface precipitation in mm/day *
3604
3605	✓✓	69110	DO il = 1, ncum
3606	✓✓✓✓	69110	IF (ep(il,inb(il))>=0.0001 .AND. iflag(il)<=1) THEN
3607	✓✗	23019	IF (cvflag_ice) THEN
3608			precip(il) = wt(il, 1)sigd(il)(water(il,1)+ice(il,1)) &
3609		23019	86400.1000./(rowl*grav)
3610			ELSE
3611			precip(il) = wt(il, 1)sigd(il)water(il, 1) &
3612			86400.1000./(rowl*grav)
3613			END IF
3614			END IF
3615			END DO
3616			! print*,'cv3_yield apres calcul precip'
3617
3618
3619			! === calculate vertical profile of precipitation in kg/m2/s ===
3620
3621	✓✓	4032	DO i = 1, nl
3622	✓✓	1866114	DO il = 1, ncum
3623	✓✓✓✓ ✓✓	1865970	IF (ep(il,inb(il))>=0.0001 .AND. i<=inb(il) .AND. iflag(il)<=1) THEN
3624	✓✗	386747	IF (cvflag_ice) THEN
3625		386747	Vprecip(il, i) = wt(il, i)sigd(il)(water(il,i)+ice(il,i))/grav
3626		386747	Vprecipi(il, i) = wt(il, i)sigd(il)ice(il,i)/grav ! jyg
3627			ELSE
3628			Vprecip(il, i) = wt(il, i)sigd(il)water(il, i)/grav
3629			Vprecipi(il, i) = 0. ! jyg
3630			END IF
3631			END IF
3632			END DO
3633			END DO
3634
3635
3636			! * Calculate downdraft velocity scale *
3637			! * NE PAS UTILISER POUR L'INSTANT *
3638
3639			!! do il=1,ncum
3640			!! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il)) &
3641			!! /(sigd(il)*p(il,icb(il)))
3642			!! enddo
3643
3644
3645			! * calculate tendencies of lowest level potential temperature *
3646			! * and mixing ratio *
3647
3648	✓✓	69110	DO il = 1, ncum
3649		68966	work(il) = 1.0/(ph(il,1)-ph(il,2))
3650		69110	cbmf(il) = 0.0
3651			END DO
3652
3653			! - Adiabatic ascent mass flux "ma" and cloud base mass flux "cbmf"
3654			!-----------------------------------------------------------------
3655			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3656	✗✓	144	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
3657			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3658			!!! Warning : this option leads to water conservation violation
3659			!!! Expert only
3660			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3661			DO il = 1, ncum
3662			ma(il, nlp) = 0.
3663			ma(il, 1) = 0.
3664			END DO
3665			DO k = nl, 2, -1
3666			DO il = 1, ncum
3667			ma(il, k) = ma(il, k+1)*(1.-qta(il, k))/(1.-qta(il, k-1)) + m(il, k)
3668			cbmf(il) = max(cbmf(il), ma(il,k))
3669			END DO
3670			END DO
3671			DO k = 2,nl
3672			DO il = 1, ncum
3673			IF (k <icb(il)) THEN
3674			ma(il, k) = ma(il, k-1) + wghti(il,k-1)*cbmf(il)
3675			ENDIF
3676			END DO
3677			END DO
3678			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3679			ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
3680			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3681			!! Line kept for compatibility with earlier versions
3682	✓✓	3888	DO k = 2, nl
3683	✓✓	1797004	DO il = 1, ncum
3684	✓✓	1796860	IF (k>=icb(il)) THEN
3685		1556118	cbmf(il) = cbmf(il) + m(il, k)
3686			END IF
3687			END DO
3688			END DO
3689
3690	✓✓	69110	DO il = 1, ncum
3691		68966	ma(il, nlp) = 0.
3692		69110	ma(il, 1) = 0.
3693			END DO
3694	✓✓	3888	DO k = nl, 2, -1
3695	✓✓	1797004	DO il = 1, ncum
3696		1796860	ma(il, k) = ma(il, k+1) + m(il, k)
3697			END DO
3698			END DO
3699	✓✓	3888	DO k = 2,nl
3700	✓✓	1797004	DO il = 1, ncum
3701	✓✓	1796860	IF (k <icb(il)) THEN
3702		236998	ma(il, k) = ma(il, k-1) + wghti(il,k-1)*cbmf(il)
3703			ENDIF
3704			END DO
3705			END DO
3706
3707			ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
3708			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3709			!
3710			! print*,'cv3_yield avant ft'
3711			! am is the part of cbmf taken from the first level
3712	✓✓	69110	DO il = 1, ncum
3713		69110	am(il) = cbmf(il)*wghti(il, 1)
3714			END DO
3715
3716	✓✓	69110	DO il = 1, ncum
3717	✓✓	69110	IF (iflag(il)<=1) THEN
3718			! convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
3719			!JYG Correction pour conserver l'eau
3720			! cc ft(il,1)=-0.5lvcp(il,1)sigd(il)*(evap(il,1)+evap(il,2)) !precip
3721	✓✗	23019	IF (cvflag_ice) THEN
3722			ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1) - &
3723			lfcp(il, 1)sigd(il)evap(il, 1)*faci(il, 1) - &
3724			lfcp(il, 1)sigd(il)(fondue(il,1)*wt(il,1)) / &
3725		23019	(100.*(ph(il,1)-ph(il,2))) !precip
3726			ELSE
3727			ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1)
3728			END IF
3729
3730		23019	ft(il, 1) = ft(il, 1) - 0.009gravsigd(il)mp(il, 2)t_wake(il, 1)b(il, 1)work(il)
3731
3732	✓✗	23019	IF (cvflag_ice) THEN
3733			ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
3734			(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
3735			0.01sigd(il)wt(il, 1)(ci-cpd)ice(il, 2) * &
3736		23019	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
3737			ELSE
3738			ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
3739			(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
3740			END IF
3741
3742		23019	ftd(il, 1) = ft(il, 1) ! fin precip
3743
3744	✗✓	23019	IF ((0.01gravwork(il)*am(il))>=delti) iflag(il) = 1 !consist vect
3745			!jyg<
3746	✗✓	23019	IF (fl_cor_ebil >= 2) THEN
3747			ft(il, 1) = ft(il, 1) + 0.01gravwork(il)am(il) &
3748			((t(il,2)-t(il,1))*cpn(il,2)+gz(il,2)-gz(il,1))/cpn(il,1)
3749			ELSE
3750			ft(il, 1) = ft(il, 1) + 0.01gravwork(il)am(il) &
3751		23019	(t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
3752			ENDIF
3753			!>jyg
3754			END IF ! iflag
3755			END DO
3756
3757
3758	✓✓	3888	DO j = 2, nl
3759	✓✗	3888	IF (iflag_mix>0) THEN
3760	✓✓	1796860	DO il = 1, ncum
3761			! FH WARNING a modifier :
3762			cpinv = 0.
3763			! cpinv=1.0/cpn(il,1)
3764	✓✓✓✓	1796860	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
3765			ft(il, 1) = ft(il, 1) + 0.01gravwork(il)ment(il, j, 1) &
3766		363728	(hent(il,j,1)-h(il,1)+t(il,1)(cpv-cpd)(rr(il,1)-qent(il,j,1)))*cpinv
3767			END IF ! j
3768			END DO
3769			END IF
3770			END DO
3771			! fin sature
3772
3773
3774	✓✓	69110	DO il = 1, ncum
3775	✓✓	69110	IF (iflag(il)<=1) THEN
3776			!JYG1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
3777			fr(il, 1) = 0.01gravmp(il, 2)(rp(il,2)-rr_wake(il,1))work(il) + &
3778		23019	sigd(il)*evap(il, 1)
3779			!!! sigd(il)0.5(evap(il,1)+evap(il,2))
3780
3781		23019	fqd(il, 1) = fr(il, 1) !precip
3782
3783		23019	fr(il, 1) = fr(il, 1) + 0.01gravam(il)(rr(il,2)-rr(il,1))work(il) !sature
3784
3785			fu(il, 1) = fu(il, 1) + 0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1)) + &
3786		23019	am(il)*(u(il,2)-u(il,1)))
3787			fv(il, 1) = fv(il, 1) + 0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1)) + &
3788		23019	am(il)*(v(il,2)-v(il,1)))
3789			END IF ! iflag
3790			END DO ! il
3791
3792
3793			!AC! do j=1,ntra
3794			!AC! do il=1,ncum
3795			!AC! if (iflag(il) .le. 1) then
3796			!AC! if (cvflag_grav) then
3797			!AC! ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
3798			!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
3799			!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
3800			!AC! else
3801			!AC! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
3802			!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
3803			!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
3804			!AC! endif
3805			!AC! endif ! iflag
3806			!AC! enddo
3807			!AC! enddo
3808
3809	✓✓	3888	DO j = 2, nl
3810	✓✓	1797004	DO il = 1, ncum
3811	✓✓✓✓	1796860	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
3812		363728	fr(il, 1) = fr(il, 1) + 0.01gravwork(il)ment(il, j, 1)(qent(il,j,1)-rr(il,1))
3813		363728	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)ment(il, j, 1)(uent(il,j,1)-u(il,1))
3814		363728	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)ment(il, j, 1)(vent(il,j,1)-v(il,1))
3815			END IF ! j
3816			END DO
3817			END DO
3818
3819			!AC! do k=1,ntra
3820			!AC! do j=2,nl
3821			!AC! do il=1,ncum
3822			!AC! if (j.le.inb(il) .and. iflag(il) .le. 1) then
3823			!AC!
3824			!AC! if (cvflag_grav) then
3825			!AC! ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
3826			!AC! : *(traent(il,j,1,k)-tra(il,1,k))
3827			!AC! else
3828			!AC! ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
3829			!AC! : *(traent(il,j,1,k)-tra(il,1,k))
3830			!AC! endif
3831			!AC!
3832			!AC! endif
3833			!AC! enddo
3834			!AC! enddo
3835			!AC! enddo
3836			! print*,'cv3_yield apres ft'
3837
3838			!jyg<
3839			!-----------------------------------------------------------
3840	✓✗	144	IF (ok_optim_yield) THEN !\|
3841			!-----------------------------------------------------------
3842			!
3843			!* *
3844			!* Compute convective mass fluxes upwd and dnwd *
3845
3846			!
3847			! =================================================
3848			! upward fluxes \|
3849			! ------------------------------------------------
3850			!
3851	✓✓✓✓	5588064	upwd(:,:) = 0.
3852	✓✓✓✓	5588064	up_to(:,:) = 0.
3853	✓✓✓✓	5588064	up_from(:,:) = 0.
3854			!
3855			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3856	✗✓	144	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
3857			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3858			!! The decrease of the adiabatic ascent mass flux due to ejection of precipitation
3859			!! is taken into account.
3860			!! WARNING : in the present version, taking into account the mass-flux decrease due to
3861			!! precipitation ejection leads to water conservation violation.
3862			!
3863			! - Upward mass flux of mixed draughts
3864			!---------------------------------------
3865			DO i = 2, nl
3866			DO j = 1, i-1
3867			DO il = 1, ncum
3868			IF (i<=inb(il)) THEN
3869			up_to(il,i) = up_to(il,i) + ment(il,j,i)
3870			ENDIF
3871			ENDDO
3872			ENDDO
3873			ENDDO
3874			!
3875			DO j = 3, nl
3876			DO i = 2, j-1
3877			DO il = 1, ncum
3878			IF (j<=inb(il)) THEN
3879			up_from(il,i) = up_from(il,i) + ment(il,i,j)
3880			ENDIF
3881			ENDDO
3882			ENDDO
3883			ENDDO
3884			!
3885			! The difference between upwd(il,i) and upwd(il,i-1) is due to updrafts ending in layer
3886			!(i-1) (theses drafts cross interface (i-1) but not interface(i)) and to updrafts starting
3887			!from layer (i-1) (theses drafts cross interface (i) but not interface(i-1)):
3888			!
3889			DO i = 2, nlp
3890			DO il = 1, ncum
3891			IF (i<=inb(il)+1) THEN
3892			upwd(il,i) = max(0., upwd(il,i-1) - up_to(il,i-1) + up_from(il,i-1))
3893			ENDIF
3894			ENDDO
3895			ENDDO
3896			!
3897			! - Total upward mass flux
3898			!---------------------------
3899			DO i = 2, nlp
3900			DO il = 1, ncum
3901			IF (i<=inb(il)+1) THEN
3902			upwd(il,i) = upwd(il,i) + ma(il,i)
3903			ENDIF
3904			ENDDO
3905			ENDDO
3906			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3907			ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
3908			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3909			!! The decrease of the adiabatic ascent mass flux due to ejection of precipitation
3910			!! is not taken into account.
3911			!
3912			! - Upward mass flux
3913			!-------------------
3914	✓✓	3888	DO i = 2, nl
3915	✓✓	1796860	DO il = 1, ncum
3916	✓✓	1796860	IF (i<=inb(il)) THEN
3917		823705	up_to(il,i) = m(il,i)
3918			ENDIF
3919			ENDDO
3920	✓✓	54432	DO j = 1, i-1
3921	✓✓	24261354	DO il = 1, ncum
3922	✓✓	24257610	IF (i<=inb(il)) THEN
3923		6037065	up_to(il,i) = up_to(il,i) + ment(il,j,i)
3924			ENDIF
3925			ENDDO
3926			ENDDO
3927			ENDDO
3928			!
3929	✓✓	4032	DO i = 1, nl
3930	✓✓	1866114	DO il = 1, ncum
3931	✓✓	1865970	IF (i<=inb(il)) THEN
3932		892671	up_from(il,i) = cbmf(il)*wghti(il,i)
3933			ENDIF
3934			ENDDO
3935			ENDDO
3936			!
3937	✓✓	3744	DO j = 3, nl
3938	✓✓	50544	DO i = 2, j-1
3939	✓✓	22464350	DO il = 1, ncum
3940	✓✓	22460750	IF (j<=inb(il)) THEN
3941		5213360	up_from(il,i) = up_from(il,i) + ment(il,i,j)
3942			ENDIF
3943			ENDDO
3944			ENDDO
3945			ENDDO
3946			!
3947			! The difference between upwd(il,i) and upwd(il,i-1) is due to updrafts ending in layer
3948			!(i-1) (theses drafts cross interface (i-1) but not interface(i)) and to updrafts starting
3949			!from layer (i-1) (theses drafts cross interface (i) but not interface(i-1)):
3950			!
3951	✓✓	4032	DO i = 2, nlp
3952	✓✓	1866114	DO il = 1, ncum
3953	✓✓	1865970	IF (i<=inb(il)+1) THEN
3954		892671	upwd(il,i) = max(0., upwd(il,i-1) - up_to(il,i-1) + up_from(il,i-1))
3955			ENDIF
3956			ENDDO
3957			ENDDO
3958
3959
3960			ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
3961			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3962
3963			!
3964			! =================================================
3965			! downward fluxes \|
3966			! ------------------------------------------------
3967	✓✓✓✓	5588064	dnwd(:,:) = 0.
3968	✓✓✓✓	5588064	dn_to(:,:) = 0.
3969	✓✓✓✓	5588064	dn_from(:,:) = 0.
3970	✓✓	4032	DO i = 1, nl
3971	✓✓	54576	DO j = i+1, nl
3972	✓✓	24261498	DO il = 1, ncum
3973	✓✓	24257610	IF (j<=inb(il)) THEN
3974			!! dn_to(il,i) = dn_to(il,i) + ment(il,j,i) !jyg,20220202
3975		6037065	dn_to(il,i) = dn_to(il,i) - ment(il,j,i)
3976			ENDIF
3977			ENDDO
3978			ENDDO
3979			ENDDO
3980			!
3981	✓✓	4032	DO j = 1, nl
3982	✓✓	54576	DO i = j+1, nl
3983	✓✓	24261498	DO il = 1, ncum
3984	✓✓	24257610	IF (i<=inb(il)) THEN
3985			!! dn_from(il,i) = dn_from(il,i) + ment(il,i,j) !jyg,20220202
3986		6037065	dn_from(il,i) = dn_from(il,i) - ment(il,i,j)
3987			ENDIF
3988			ENDDO
3989			ENDDO
3990			ENDDO
3991			!
3992			! The difference between dnwd(il,i) and dnwd(il,i+1) is due to downdrafts ending in layer
3993			!(i) (theses drafts cross interface (i+1) but not interface(i)) and to downdrafts
3994			!starting from layer (i) (theses drafts cross interface (i) but not interface(i+1)):
3995			!
3996	✓✓	3888	DO i = nl-1, 1, -1
3997	✓✓	1797004	DO il = 1, ncum
3998			!! dnwd(il,i) = max(0., dnwd(il,i+1) - dn_to(il,i) + dn_from(il,i)) !jyg,20220202
3999		1796860	dnwd(il,i) = min(0., dnwd(il,i+1) - dn_to(il,i) + dn_from(il,i))
4000			ENDDO
4001			ENDDO
4002			! =================================================
4003			!
4004			!-----------------------------------------------------------
4005			ENDIF !(ok_optim_yield) !\|
4006			!-----------------------------------------------------------
4007			!>jyg
4008
4009			! * calculate tendencies of potential temperature and mixing ratio *
4010			! * at levels above the lowest level *
4011
4012			! * first find the net saturated updraft and downdraft mass fluxes *
4013			! * through each level *
4014
4015
4016			!jyg<
4017			!! DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
4018	✓✓	3888	DO i = 2, nl
4019			!>jyg
4020
4021			num1 = 0
4022	✓✓	1796860	DO il = 1, ncum
4023	✓✓✓✓	1796860	IF (i<=inb(il) .AND. iflag(il)<=1) num1 = num1 + 1
4024			END DO
4025	✓✓	3744	IF (num1<=0) GO TO 500
4026
4027			!
4028			!jyg<
4029			!-----------------------------------------------------------
4030	✓✗	2883	IF (ok_optim_yield) THEN !\|
4031			!-----------------------------------------------------------
4032	✓✓	1383804	DO il = 1, ncum
4033		1380921	amp1(il) = upwd(il,i+1)
4034		1383804	ad(il) = dnwd(il,i)
4035			ENDDO
4036			!-----------------------------------------------------------
4037			ELSE !(ok_optim_yield) !\|
4038			!-----------------------------------------------------------
4039			!>jyg
4040			DO il = 1,ncum
4041			amp1(il) = 0.
4042			ad(il) = 0.
4043			ENDDO
4044
4045			DO k = 1, nl + 1
4046			DO il = 1, ncum
4047			IF (i>=icb(il)) THEN
4048			IF (k>=i+1 .AND. k<=(inb(il)+1)) THEN
4049			amp1(il) = amp1(il) + m(il, k)
4050			END IF
4051			ELSE
4052			! AMP1 is the part of cbmf taken from layers I and lower
4053			IF (k<=i) THEN
4054			amp1(il) = amp1(il) + cbmf(il)*wghti(il, k)
4055			END IF
4056			END IF
4057			END DO
4058			END DO
4059
4060			DO j = i + 1, nl + 1
4061			DO k = 1, i
4062			!yor! reverted j and k loops
4063			DO il = 1, ncum
4064			!yor! IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN ! the second condition implies the first !
4065			IF (j<=(inb(il)+1)) THEN
4066			amp1(il) = amp1(il) + ment(il, k, j)
4067			END IF
4068			END DO
4069			END DO
4070			END DO
4071
4072			DO k = 1, i - 1
4073			!jyg<
4074			!! DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
4075			DO j = i, nl
4076			!>jyg
4077			DO il = 1, ncum
4078			!yor! IF (i<=inb(il) .AND. j<=inb(il)) THEN ! the second condition implies the 1st !
4079			IF (j<=inb(il)) THEN
4080			ad(il) = ad(il) + ment(il, j, k)
4081			END IF
4082			END DO
4083			END DO
4084			END DO
4085			!
4086			!-----------------------------------------------------------
4087			ENDIF !(ok_optim_yield) !\|
4088			!-----------------------------------------------------------
4089			!
4090			!! print *,'yield, i, amp1, ad', i, amp1(1), ad(1)
4091
4092	✓✓	1383804	DO il = 1, ncum
4093	✓✓✓✓	1383804	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
4094		363728	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4095		363728	cpinv = 1.0/cpn(il, i)
4096
4097			! convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
4098	✓✓	363728	IF ((0.01gravdpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
4099
4100			! precip
4101			! cc ft(il,i)= -0.5sigd(il)lvcp(il,i)*(evap(il,i)+evap(il,i+1))
4102	✓✗	363728	IF (cvflag_ice) THEN
4103			ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i) - &
4104			sigd(il)lfcp(il, i)evap(il, i)*faci(il, i) - &
4105		363728	sigd(il)lfcp(il, i)fondue(il, i)wt(il, i)/(100.(p(il,i-1)-p(il,i)))
4106			ELSE
4107			ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i)
4108			END IF
4109
4110		363728	rat = cpn(il, i-1)*cpinv
4111
4112			ft(il, i) = ft(il, i) - 0.009gravsigd(il) * &
4113		363728	(mp(il,i+1)t_wake(il,i)b(il,i)-mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
4114	✓✗	363728	IF (cvflag_ice) THEN
4115			ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
4116			(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv + &
4117			0.01sigd(il)wt(il, i)(ci-cpd)ice(il, i+1) * &
4118		363728	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv
4119			ELSE
4120			ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
4121			(t_wake(il,i+1)-t_wake(il,i))dpinv &
4122			cpinv
4123			END IF
4124
4125		363728	ftd(il, i) = ft(il, i)
4126			! fin precip
4127
4128			! sature
4129			!jyg<
4130	✗✓	363728	IF (fl_cor_ebil >= 2) THEN
4131			ft(il, i) = ft(il, i) + 0.01gravdpinv * &
4132			( amp1(il)( (t(il,i+1)-t(il,i))cpn(il,i+1) + gz(il,i+1)-gz(il,i))*cpinv - &
4133			ad(il)( (t(il,i)-t(il,i-1))cpn(il,i-1) + gz(il,i)-gz(il,i-1))*cpinv)
4134			ELSE
4135			ft(il, i) = ft(il, i) + 0.01gravdpinv * &
4136			(amp1(il)(t(il,i+1)-t(il,i) + (gz(il,i+1)-gz(il,i))cpinv) - &
4137		363728	ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
4138			ENDIF
4139			!>jyg
4140
4141
4142	✗✓	363728	IF (iflag_mix==0) THEN
4143			ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, i, i)(hp(il,i)-h(il,i) + &
4144			t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
4145			END IF
4146			!
4147			! sb: on ne fait pas encore la correction permettant de mieux
4148			! conserver l'eau:
4149			!JYG: correction permettant de mieux conserver l'eau:
4150			! cc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
4151			fr(il, i) = sigd(il)evap(il, i) + 0.01grav(mp(il,i+1)(rp(il,i+1)-rr_wake(il,i)) - &
4152		363728	mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
4153		363728	fqd(il, i) = fr(il, i) ! precip
4154
4155			fu(il, i) = 0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i)) - &
4156		363728	mp(il,i)(up(il,i)-u(il,i-1)))dpinv
4157			fv(il, i) = 0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i)) - &
4158		363728	mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
4159
4160
4161			fr(il, i) = fr(il, i) + 0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i)) - &
4162		363728	ad(il)*(rr(il,i)-rr(il,i-1)))
4163			fu(il, i) = fu(il, i) + 0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i)) - &
4164		363728	ad(il)*(u(il,i)-u(il,i-1)))
4165			fv(il, i) = fv(il, i) + 0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i)) - &
4166		363728	ad(il)*(v(il,i)-v(il,i-1)))
4167
4168			END IF ! i
4169			END DO
4170
4171			!AC! do k=1,ntra
4172			!AC! do il=1,ncum
4173			!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
4174			!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
4175			!AC! cpinv=1.0/cpn(il,i)
4176			!AC! if (cvflag_grav) then
4177			!AC! ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
4178			!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
4179			!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
4180			!AC! else
4181			!AC! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
4182			!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
4183			!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
4184			!AC! endif
4185			!AC! endif
4186			!AC! enddo
4187			!AC! enddo
4188
4189	✓✓	33186	DO k = 1, i - 1
4190
4191	✓✓	14546784	DO il = 1, ncum
4192		14516481	awat(il) = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
4193		14546784	awat(il) = max(awat(il), 0.0)
4194			END DO
4195
4196	✓✗	30303	IF (iflag_mix/=0) THEN
4197	✓✓	14546784	DO il = 1, ncum
4198	✓✓✓✓	14546784	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
4199		3192455	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4200		3192455	cpinv = 1.0/cpn(il, i)
4201			ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
4202		3192455	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-qent(il,k,i)))*cpinv
4203			!
4204			!
4205			END IF ! i
4206			END DO
4207			END IF
4208
4209	✓✓	14549667	DO il = 1, ncum
4210	✓✓✓✓	14546784	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
4211		3192455	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4212			cpinv = 1.0/cpn(il, i)
4213			fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i) &
4214		3192455	(qent(il,k,i)-awat(il)-rr(il,i))
4215		3192455	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
4216		3192455	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
4217
4218			! (saturated updrafts resulting from mixing) ! cld
4219		3192455	qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il)) ! cld
4220		3192455	qtment(il, i) = qtment(il, i) + qent(il,k,i) ! cld
4221		3192455	nqcond(il, i) = nqcond(il, i) + 1. ! cld
4222			END IF ! i
4223			END DO
4224			END DO
4225
4226			!AC! do j=1,ntra
4227			!AC! do k=1,i-1
4228			!AC! do il=1,ncum
4229			!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
4230			!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
4231			!AC! cpinv=1.0/cpn(il,i)
4232			!AC! if (cvflag_grav) then
4233			!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
4234			!AC! : *(traent(il,k,i,j)-tra(il,i,j))
4235			!AC! else
4236			!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
4237			!AC! : *(traent(il,k,i,j)-tra(il,i,j))
4238			!AC! endif
4239			!AC! endif
4240			!AC! enddo
4241			!AC! enddo
4242			!AC! enddo
4243
4244			!jyg<
4245			!! DO k = i, nl + 1
4246	✓✓	50421	DO k = i, nl
4247			!>jyg
4248
4249	✓✗	47538	IF (iflag_mix/=0) THEN
4250	✓✓	22815924	DO il = 1, ncum
4251	✓✓✓✓ ✓✓	22815924	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
4252		3192455	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4253		3192455	cpinv = 1.0/cpn(il, i)
4254			ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
4255		3192455	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)-qent(il,k,i)))*cpinv
4256
4257
4258			END IF ! i
4259			END DO
4260			END IF
4261
4262	✓✓	22818807	DO il = 1, ncum
4263	✓✓✓✓ ✓✓	22815924	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
4264		3192455	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4265			cpinv = 1.0/cpn(il, i)
4266
4267		3192455	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k,i)-rr(il,i))
4268		3192455	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
4269		3192455	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
4270			END IF ! i and k
4271			END DO
4272			END DO
4273
4274			!AC! do j=1,ntra
4275			!AC! do k=i,nl+1
4276			!AC! do il=1,ncum
4277			!AC! if (i.le.inb(il) .and. k.le.inb(il)
4278			!AC! $ .and. iflag(il) .le. 1) then
4279			!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
4280			!AC! cpinv=1.0/cpn(il,i)
4281			!AC! if (cvflag_grav) then
4282			!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
4283			!AC! : *(traent(il,k,i,j)-tra(il,i,j))
4284			!AC! else
4285			!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
4286			!AC! : *(traent(il,k,i,j)-tra(il,i,j))
4287			!AC! endif
4288			!AC! endif ! i and k
4289			!AC! enddo
4290			!AC! enddo
4291			!AC! enddo
4292
4293			! sb: interface with the cloud parameterization: ! cld
4294
4295	✓✓	47538	DO k = i + 1, nl
4296	✓✓	21435003	DO il = 1, ncum
4297	✓✓✓✗ ✓✓	21432120	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
4298			! (saturated downdrafts resulting from mixing) ! cld
4299		2828727	qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
4300		2828727	qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
4301		2828727	nqcond(il, i) = nqcond(il, i) + 1. ! cld
4302			END IF ! cld
4303			END DO ! cld
4304			END DO ! cld
4305
4306			!ym BIG Warning : it seems that the k loop is missing !!!
4307			!ym Strong advice to check this
4308			!ym add a k loop temporary
4309
4310			! (particular case: no detraining level is found) ! cld
4311			! Verif merge Dynamico<<<<<<< .working
4312	✓✓	1383804	DO il = 1, ncum ! cld
4313	✓✓✓✓ ✓✓	1383804	IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
4314		82931	qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
4315			!jyg< Bug correction 20180620
4316			! PROBLEM: Should not qent(il,i,i) be taken into account even if nent(il,i)/=0?
4317			!! qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
4318		82931	qtment(il, i) = qent(il,i,i) + qtment(il,i) ! cld
4319			!>jyg
4320		82931	nqcond(il, i) = nqcond(il, i) + 1. ! cld
4321			END IF ! cld
4322			END DO ! cld
4323			! Verif merge Dynamico =======
4324			! Verif merge Dynamico DO k = i + 1, nl
4325			! Verif merge Dynamico DO il = 1, ncum !ym k loop added ! cld
4326			! Verif merge Dynamico IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
4327			! Verif merge Dynamico qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
4328			! Verif merge Dynamico qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
4329			! Verif merge Dynamico nqcond(il, i) = nqcond(il, i) + 1. ! cld
4330			! Verif merge Dynamico END IF ! cld
4331			! Verif merge Dynamico END DO
4332			! Verif merge Dynamico ENDDO ! cld
4333			! Verif merge Dynamico >>>>>>> .merge-right.r3413
4334
4335	✓✓	1383804	DO il = 1, ncum ! cld
4336	✓✓✓✓ ✓✗	1384665	IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN ! cld
4337		363728	qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
4338		363728	qtment(il, i) = qtment(il,i)/nqcond(il, i) ! cld
4339			END IF ! cld
4340			END DO
4341
4342			!AC! do j=1,ntra
4343			!AC! do il=1,ncum
4344			!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
4345			!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
4346			!AC! cpinv=1.0/cpn(il,i)
4347			!AC!
4348			!AC! if (cvflag_grav) then
4349			!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
4350			!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
4351			!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
4352			!AC! else
4353			!AC! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
4354			!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
4355			!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
4356			!AC! endif
4357			!AC! endif ! i
4358			!AC! enddo
4359			!AC! enddo
4360
4361
4362		144	500 END DO
4363
4364			!JYG<
4365			!Conservation de l'eau
4366			! sumdq = 0.
4367			! DO k = 1, nl
4368			! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
4369			! END DO
4370			! PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
4371			!JYG>
4372			! * move the detrainment at level inb down to level inb-1 *
4373			! * in such a way as to preserve the vertically *
4374			! * integrated enthalpy and water tendencies *
4375
4376			! Correction bug le 18-03-09
4377	✓✓	69110	DO il = 1, ncum
4378	✓✓	69110	IF (iflag(il)<=1) THEN
4379			ax = 0.01gravment(il, inb(il), inb(il))* &
4380			(hp(il,inb(il))-h(il,inb(il))+t(il,inb(il))(cpv-cpd)(rr(il,inb(il))-qent(il,inb(il),inb(il))))/ &
4381		23019	(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
4382		23019	ft(il, inb(il)) = ft(il, inb(il)) - ax
4383			ft(il, inb(il)-1) = ft(il, inb(il)-1) + axcpn(il, inb(il))(ph(il,inb(il))-ph(il,inb(il)+1))/ &
4384		23019	(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il,inb(il))))
4385
4386			bx = 0.01gravment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb(il)))/ &
4387		23019	(ph(il,inb(il))-ph(il,inb(il)+1))
4388		23019	fr(il, inb(il)) = fr(il, inb(il)) - bx
4389			fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
4390		23019	(ph(il,inb(il)-1)-ph(il,inb(il)))
4391
4392			cx = 0.01gravment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il)))/ &
4393		23019	(ph(il,inb(il))-ph(il,inb(il)+1))
4394		23019	fu(il, inb(il)) = fu(il, inb(il)) - cx
4395			fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
4396		23019	(ph(il,inb(il)-1)-ph(il,inb(il)))
4397
4398			dx = 0.01gravment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il)))/ &
4399		23019	(ph(il,inb(il))-ph(il,inb(il)+1))
4400		23019	fv(il, inb(il)) = fv(il, inb(il)) - dx
4401			fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
4402		23019	(ph(il,inb(il)-1)-ph(il,inb(il)))
4403			END IF !iflag
4404			END DO
4405
4406			!JYG<
4407			!Conservation de l'eau
4408			! sumdq = 0.
4409			! DO k = 1, nl
4410			! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
4411			! END DO
4412			! PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
4413			!JYG>
4414
4415			!AC! do j=1,ntra
4416			!AC! do il=1,ncum
4417			!AC! IF (iflag(il) .le. 1) THEN
4418			!AC! IF (cvflag_grav) then
4419			!AC! ex=0.01gravment(il,inb(il),inb(il))
4420			!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
4421			!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
4422			!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
4423			!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
4424			!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
4425			!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
4426			!AC! else
4427			!AC! ex=0.1*ment(il,inb(il),inb(il))
4428			!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
4429			!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
4430			!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
4431			!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
4432			!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
4433			!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
4434			!AC! ENDIF !cvflag grav
4435			!AC! ENDIF !iflag
4436			!AC! enddo
4437			!AC! enddo
4438
4439
4440			! * homogenize tendencies below cloud base *
4441
4442
4443	✓✓	69110	DO il = 1, ncum
4444		68966	asum(il) = 0.0
4445		68966	bsum(il) = 0.0
4446		68966	csum(il) = 0.0
4447		68966	dsum(il) = 0.0
4448		68966	esum(il) = 0.0
4449		68966	fsum(il) = 0.0
4450		68966	gsum(il) = 0.0
4451		69110	hsum(il) = 0.0
4452			END DO
4453
4454			!do i=1,nl
4455			!do il=1,ncum
4456			!th_wake(il,i)=t_wake(il,i)(1000.0/p(il,i))*rdcp
4457			!enddo
4458			!enddo
4459
4460	✓✓	4032	DO i = 1, nl
4461	✓✓	1866114	DO il = 1, ncum
4462	✓✓✓✓	1865970	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
4463			!jyg Saturated part : use T profile
4464		101189	asum(il) = asum(il) + (ft(il,i)-ftd(il,i))*(ph(il,i)-ph(il,i+1))
4465			!jyg<20140311
4466			!Correction pour conserver l eau
4467	✓✗	101189	IF (ok_conserv_q) THEN
4468		101189	bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(ph(il,i)-ph(il,i+1))
4469		101189	csum(il) = csum(il) + (ph(il,i)-ph(il,i+1))
4470
4471			ELSE
4472			bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))* &
4473			(ph(il,i)-ph(il,i+1))
4474			csum(il)=csum(il)+(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1))) &
4475			(ph(il,i)-ph(il,i+1))
4476			ENDIF ! (ok_conserv_q)
4477			!jyg>
4478		101189	dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
4479			!jyg Unsaturated part : use T_wake profile
4480		101189	esum(il) = esum(il) + ftd(il, i)*(ph(il,i)-ph(il,i+1))
4481			!jyg<20140311
4482			!Correction pour conserver l eau
4483	✓✗	101189	IF (ok_conserv_q) THEN
4484		101189	fsum(il) = fsum(il) + fqd(il, i)*(ph(il,i)-ph(il,i+1))
4485		101189	gsum(il) = gsum(il) + (ph(il,i)-ph(il,i+1))
4486			ELSE
4487			fsum(il)=fsum(il)+fqd(il,i)(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1)))* &
4488			(ph(il,i)-ph(il,i+1))
4489			gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1))) &
4490			(ph(il,i)-ph(il,i+1))
4491			ENDIF ! (ok_conserv_q)
4492			!jyg>
4493		101189	hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, i)
4494			END IF
4495			END DO
4496			END DO
4497
4498			!!!! do 700 i=1,icb(il)-1
4499	✓✗	144	IF (ok_homo_tend) THEN
4500	✓✓	4032	DO i = 1, nl
4501	✓✓	1866114	DO il = 1, ncum
4502	✓✓✓✓	1865970	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
4503		101189	ftd(il, i) = esum(il)t_wake(il, i)/(th_wake(il,i)hsum(il))
4504		101189	fqd(il, i) = fsum(il)/gsum(il)
4505		101189	ft(il, i) = ftd(il, i) + asum(il)t(il, i)/(th(il,i)dsum(il))
4506		101189	fr(il, i) = fqd(il, i) + bsum(il)/csum(il)
4507			END IF
4508			END DO
4509			END DO
4510			ENDIF
4511
4512			!jyg<
4513			!Conservation de l'eau
4514			!! sumdq = 0.
4515			!! DO k = 1, nl
4516			!! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
4517			!! END DO
4518			!! PRINT *, 'cv3_yield, apres hom, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
4519			!jyg>
4520
4521
4522			! *** Check that moisture stays positive. If not, scale tendencies
4523			! in order to ensure moisture positivity
4524	✓✓	69110	DO il = 1, ncum
4525		68966	alpha_qpos(il) = 1.
4526	✓✓	69110	IF (iflag(il)<=1) THEN
4527	✓✓	23019	IF (fr(il,1)<=0.) THEN
4528		4849	alpha_qpos(il) = max(alpha_qpos(il), (-deltfr(il,1))/(s_wake(il)rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
4529			END IF
4530			END IF
4531			END DO
4532	✓✓	3888	DO i = 2, nl
4533	✓✓	1797004	DO il = 1, ncum
4534	✓✓	1796860	IF (iflag(il)<=1) THEN
4535	✓✓	598494	IF (fr(il,i)<=0.) THEN
4536		437433	alpha_qpos1(il) = max(1., (-deltfr(il,i))/(s_wake(il)rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
4537	✓✗	437433	IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) = alpha_qpos1(il)
4538			END IF
4539			END IF
4540			END DO
4541			END DO
4542	✓✓	69110	DO il = 1, ncum
4543	✓✓✗✓	69110	IF (iflag(il)<=1 .AND. alpha_qpos(il)>1.001) THEN
4544			alpha_qpos(il) = alpha_qpos(il)*1.1
4545			END IF
4546			END DO
4547			!
4548	✗✓	144	IF (prt_level .GE. 5) THEN
4549			print *,' CV3_YIELD : alpha_qpos ',alpha_qpos(1)
4550			ENDIF
4551			!
4552	✓✓	69110	DO il = 1, ncum
4553	✓✓	69110	IF (iflag(il)<=1) THEN
4554		23019	sigd(il) = sigd(il)/alpha_qpos(il)
4555		23019	precip(il) = precip(il)/alpha_qpos(il)
4556		23019	cbmf(il) = cbmf(il)/alpha_qpos(il)
4557			END IF
4558			END DO
4559	✓✓	4032	DO i = 1, nl
4560	✓✓	1866114	DO il = 1, ncum
4561	✓✓	1865970	IF (iflag(il)<=1) THEN
4562		621513	fr(il, i) = fr(il, i)/alpha_qpos(il)
4563		621513	ft(il, i) = ft(il, i)/alpha_qpos(il)
4564		621513	fqd(il, i) = fqd(il, i)/alpha_qpos(il)
4565		621513	ftd(il, i) = ftd(il, i)/alpha_qpos(il)
4566		621513	fu(il, i) = fu(il, i)/alpha_qpos(il)
4567		621513	fv(il, i) = fv(il, i)/alpha_qpos(il)
4568		621513	m(il, i) = m(il, i)/alpha_qpos(il)
4569		621513	mp(il, i) = mp(il, i)/alpha_qpos(il)
4570		621513	Vprecip(il, i) = Vprecip(il, i)/alpha_qpos(il)
4571		621513	Vprecipi(il, i) = Vprecipi(il, i)/alpha_qpos(il) ! jyg
4572			END IF
4573			END DO
4574			END DO
4575			!jyg<
4576			!-----------------------------------------------------------
4577	✓✗	144	IF (ok_optim_yield) THEN !\|
4578			!-----------------------------------------------------------
4579	✓✓	4032	DO i = 1, nl
4580	✓✓	1866114	DO il = 1, ncum
4581	✓✓	1865970	IF (iflag(il)<=1) THEN
4582		621513	upwd(il, i) = upwd(il, i)/alpha_qpos(il)
4583		621513	dnwd(il, i) = dnwd(il, i)/alpha_qpos(il)
4584			END IF
4585			END DO
4586			END DO
4587			!-----------------------------------------------------------
4588			ENDIF !(ok_optim_yield) !\|
4589			!-----------------------------------------------------------
4590			!>jyg
4591	✓✓	4032	DO j = 1, nl !yor! inverted i and j loops
4592	✓✓	109008	DO i = 1, nl
4593	✓✓	50385078	DO il = 1, ncum
4594	✓✓	50381190	IF (iflag(il)<=1) THEN
4595		16780851	ment(il, i, j) = ment(il, i, j)/alpha_qpos(il)
4596			END IF
4597			END DO
4598			END DO
4599			END DO
4600
4601			!AC! DO j = 1,ntra
4602			!AC! DO i = 1,nl
4603			!AC! DO il = 1,ncum
4604			!AC! IF (iflag(il) .le. 1) THEN
4605			!AC! ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
4606			!AC! ENDIF
4607			!AC! ENDDO
4608			!AC! ENDDO
4609			!AC! ENDDO
4610
4611
4612			! * reset counter and return *
4613
4614			! Reset counter only for points actually convective (jyg)
4615			! In order take into account the possibility of changing the compression,
4616			! reset m, sig and w0 to zero for non-convecting points.
4617	✓✓	69110	DO il = 1, ncum
4618	✓✓	69110	IF (iflag(il) < 3) THEN
4619		23019	sig(il, nd) = 2.0
4620			ENDIF
4621			END DO
4622
4623
4624	✓✓	4032	DO i = 1, nl
4625	✓✓	1866114	DO il = 1, ncum
4626		1865970	dnwd0(il, i) = -mp(il, i)
4627			END DO
4628			END DO
4629			!jyg< (loops stop at nl)
4630			!! DO i = nl + 1, nd
4631			!! DO il = 1, ncum
4632			!! dnwd0(il, i) = 0.
4633			!! END DO
4634			!! END DO
4635			!>jyg
4636
4637
4638			!jyg<
4639			!-----------------------------------------------------------
4640	✗✓	144	IF (.NOT.ok_optim_yield) THEN !\|
4641			!-----------------------------------------------------------
4642			DO i = 1, nl
4643			DO il = 1, ncum
4644			upwd(il, i) = 0.0
4645			dnwd(il, i) = 0.0
4646			END DO
4647			END DO
4648
4649			!! DO i = 1, nl ! useless; jyg
4650			!! DO il = 1, ncum ! useless; jyg
4651			!! IF (i>=icb(il) .AND. i<=inb(il)) THEN ! useless; jyg
4652			!! upwd(il, i) = 0.0 ! useless; jyg
4653			!! dnwd(il, i) = 0.0 ! useless; jyg
4654			!! END IF ! useless; jyg
4655			!! END DO ! useless; jyg
4656			!! END DO ! useless; jyg
4657
4658			DO i = 1, nl
4659			DO k = 1, nl
4660			DO il = 1, ncum
4661			up1(il, k, i) = 0.0
4662			dn1(il, k, i) = 0.0
4663			END DO
4664			END DO
4665			END DO
4666
4667			!yor! commented original
4668			! DO i = 1, nl
4669			! DO k = i, nl
4670			! DO n = 1, i - 1
4671			! DO il = 1, ncum
4672			! IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
4673			! up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
4674			! dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
4675			! END IF
4676			! END DO
4677			! END DO
4678			! END DO
4679			! END DO
4680			!yor! replaced with
4681			DO i = 1, nl
4682			DO k = i, nl
4683			DO n = 1, i - 1
4684			DO il = 1, ncum
4685			IF (i>=icb(il) .AND. k<=inb(il)) THEN ! yor ! as i always <= k
4686			up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
4687			END IF
4688			END DO
4689			END DO
4690			END DO
4691			END DO
4692			DO i = 1, nl
4693			DO n = 1, i - 1
4694			DO k = i, nl
4695			DO il = 1, ncum
4696			IF (i>=icb(il) .AND. k<=inb(il)) THEN ! yor ! i always <= k
4697			dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
4698			END IF
4699			END DO
4700			END DO
4701			END DO
4702			END DO
4703			!yor! end replace
4704
4705			DO i = 1, nl
4706			DO k = 1, nl
4707			DO il = 1, ncum
4708			IF (i>=icb(il)) THEN
4709			IF (k>=i .AND. k<=(inb(il))) THEN
4710			upwd(il, i) = upwd(il, i) + m(il, k)
4711			END IF
4712			ELSE
4713			IF (k<i) THEN
4714			upwd(il, i) = upwd(il, i) + cbmf(il)*wghti(il, k)
4715			END IF
4716			END IF
4717			! c print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
4718			END DO
4719			END DO
4720			END DO
4721
4722			DO i = 2, nl
4723			DO k = i, nl
4724			DO il = 1, ncum
4725			! test if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
4726			IF (i<=inb(il) .AND. k<=inb(il)) THEN
4727			upwd(il, i) = upwd(il, i) + up1(il, k, i)
4728			dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
4729			END IF
4730			! c print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
4731			END DO
4732			END DO
4733			END DO
4734
4735
4736			!!!! DO il=1,ncum
4737			!!!! do i=icb(il),inb(il)
4738			!!!!
4739			!!!! upwd(il,i)=0.0
4740			!!!! dnwd(il,i)=0.0
4741			!!!! do k=i,inb(il)
4742			!!!! up1=0.0
4743			!!!! dn1=0.0
4744			!!!! do n=1,i-1
4745			!!!! up1=up1+ment(il,n,k)
4746			!!!! dn1=dn1-ment(il,k,n)
4747			!!!! enddo
4748			!!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
4749			!!!! dnwd(il,i)=dnwd(il,i)+dn1
4750			!!!! enddo
4751			!!!! enddo
4752			!!!!
4753			!!!! ENDDO
4754
4755			!! DO i = 1, nlp
4756			!! DO il = 1, ncum
4757			!! ma(il, i) = 0
4758			!! END DO
4759			!! END DO
4760			!!
4761			!! DO i = 1, nl
4762			!! DO j = i, nl
4763			!! DO il = 1, ncum
4764			!! ma(il, i) = ma(il, i) + m(il, j)
4765			!! END DO
4766			!! END DO
4767			!! END DO
4768
4769			!jyg< (loops stop at nl)
4770			!! DO i = nl + 1, nd
4771			!! DO il = 1, ncum
4772			!! ma(il, i) = 0.
4773			!! END DO
4774			!! END DO
4775			!>jyg
4776
4777			!! DO i = 1, nl
4778			!! DO il = 1, ncum
4779			!! IF (i<=(icb(il)-1)) THEN
4780			!! ma(il, i) = 0
4781			!! END IF
4782			!! END DO
4783			!! END DO
4784
4785			!-----------------------------------------------------------
4786			ENDIF !(.NOT.ok_optim_yield) !\|
4787			!-----------------------------------------------------------
4788			!>jyg
4789
4790			! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4791			! determination de la variation de flux ascendant entre
4792			! deux niveau non dilue mip
4793			! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4794
4795	✓✓	4032	DO i = 1, nl
4796	✓✓	1866114	DO il = 1, ncum
4797		1865970	mip(il, i) = m(il, i)
4798			END DO
4799			END DO
4800
4801			!jyg< (loops stop at nl)
4802			!! DO i = nl + 1, nd
4803			!! DO il = 1, ncum
4804			!! mip(il, i) = 0.
4805			!! END DO
4806			!! END DO
4807			!>jyg
4808
4809
4810			! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4811			! icb represente de niveau ou se trouve la
4812			! base du nuage , et inb le top du nuage
4813			! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4814
4815			!! DO i = 1, nd ! unused . jyg
4816			!! DO il = 1, ncum ! unused . jyg
4817			!! mke(il, i) = upwd(il, i) + dnwd(il, i) ! unused . jyg
4818			!! END DO ! unused . jyg
4819			!! END DO ! unused . jyg
4820
4821			!! DO i = 1, nd ! unused . jyg
4822			!! DO il = 1, ncum ! unused . jyg
4823			!! rdcp = (rrd(1.-rr(il,i))-rr(il,i)rrv)/(cpd(1.-rr(il,i))+rr(il,i)cpv) ! unused . jyg
4824			!! tls(il, i) = t(il, i)(1000.0/p(il,i))*rdcp ! unused . jyg
4825			!! tps(il, i) = tp(il, i) ! unused . jyg
4826			!! END DO ! unused . jyg
4827			!! END DO ! unused . jyg
4828
4829
4830			! * diagnose the in-cloud mixing ratio * ! cld
4831			! * of condensed water * ! cld
4832			!! cld
4833
4834	✓✓	4176	DO i = 1, nl+1 ! cld
4835	✓✓	1935224	DO il = 1, ncum ! cld
4836		1931048	mac(il, i) = 0.0 ! cld
4837		1931048	wa(il, i) = 0.0 ! cld
4838		1931048	siga(il, i) = 0.0 ! cld
4839		1935080	sax(il, i) = 0.0 ! cld
4840			END DO ! cld
4841			END DO ! cld
4842
4843	✓✓	4032	DO i = minorig, nl ! cld
4844	✓✓	58464	DO k = i + 1, nl + 1 ! cld
4845	✓✓	26127468	DO il = 1, ncum ! cld
4846	✓✓✓✓ ✓✓	26123580	IF (i<=inb(il) .AND. k<=(inb(il)+1) .AND. iflag(il)<=1) THEN ! cld
4847		3579202	mac(il, i) = mac(il, i) + m(il, k) ! cld
4848			END IF ! cld
4849			END DO ! cld
4850			END DO ! cld
4851			END DO ! cld
4852
4853	✓✓	4032	DO i = 1, nl ! cld
4854	✓✓	58464	DO j = 1, i ! cld
4855	✓✓	26127468	DO il = 1, ncum ! cld
4856			IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
4857	✓✓✓✓ ✓✓✓✓	26123580	.AND. j>=icb(il) .AND. iflag(il)<=1) THEN ! cld
4858			sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
4859		1790910	*(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
4860			END IF ! cld
4861			END DO ! cld
4862			END DO ! cld
4863			END DO ! cld
4864
4865	✓✓	4032	DO i = 1, nl ! cld
4866	✓✓	1866114	DO il = 1, ncum ! cld
4867			IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
4868	✓✓✓✓ ✓✓✓✗	1865970	.AND. sax(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld
4869		230071	wa(il, i) = sqrt(2.*sax(il,i)) ! cld
4870			END IF ! cld
4871			END DO ! cld
4872			END DO
4873			! cld
4874	✓✓	4032	DO i = 1, nl
4875
4876			! 14/01/15 AJ je remets les parties manquantes cf JYG
4877			! Initialize sument to 0
4878
4879	✓✓	1865970	DO il = 1,ncum
4880		1865970	sument(il) = 0.
4881			ENDDO
4882
4883			! Sum mixed mass fluxes in sument
4884
4885	✓✓	108864	DO k = 1,nl
4886	✓✓	50385078	DO il = 1,ncum
4887	✓✓✓✓ ✓✓	50381190	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
4888		6771657	sument(il) =sument(il) + abs(ment(il,k,i))
4889			ENDIF
4890			ENDDO ! il
4891			ENDDO ! k
4892
4893			! 14/01/15 AJ delta n'a rien � faire l�...
4894	✓✓	1866114	DO il = 1, ncum ! cld
4895			!! IF (wa(il,i)>0.0 .AND. iflag(il)<=1) & ! cld
4896			!! siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) & ! cld
4897			!! rrdtvp(il, i)/p(il, i)/100. ! cld
4898			!!
4899			!! siga(il, i) = min(siga(il,i), 1.0) ! cld
4900			sigaq = 0.
4901	✓✓✓✗	1862082	IF (wa(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld
4902			siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) & ! cld
4903		230071	rrdtvp(il, i)/p(il, i)/100. ! cld
4904		230071	siga(il, i) = min(siga(il,i), 1.0) ! cld
4905		230071	sigaq = siga(il,i)*qta(il,i-1) ! cld
4906			ENDIF
4907
4908			! IM cf. FH
4909			! 14/01/15 AJ ne correspond pas � ce qui a �t� cod� par JYG et SB
4910
4911	✓✗	1865970	IF (iflag_clw==0) THEN ! cld
4912			qcondc(il, i) = siga(il, i)clw(il, i)(1.-ep(il,i)) & ! cld
4913		1862082	+(1.-siga(il,i))*qcond(il, i) ! cld
4914
4915
4916		1862082	sigment(il,i)=sument(il)*tau_cld_cv/(ph(il,i)-ph(il,i+1)) ! cld
4917		1862082	sigment(il, i) = min(1.e-4+sigment(il,i), 1.0 - siga(il,i)) ! cld
4918			!! qtc(il, i) = (siga(il,i)qta(il,i-1)+sigment(il,i)qtment(il,i)) & ! cld
4919			qtc(il, i) = (sigaq+sigment(il,i)*qtment(il,i)) & ! cld
4920		1862082	/(siga(il,i)+sigment(il,i)) ! cld
4921		1862082	sigt(il,i) = sigment(il, i) + siga(il, i)
4922
4923			! qtc(il, i) = siga(il,i)qta(il,i-1)+(1.-siga(il,i))qtment(il,i) ! cld
4924			! print*,'BIGAUSSIAN CONV',siga(il,i),sigment(il,i),qtc(il,i)
4925
4926			ELSE IF (iflag_clw==1) THEN ! cld
4927			qcondc(il, i) = qcond(il, i) ! cld
4928			qtc(il,i) = qtment(il,i) ! cld
4929			END IF ! cld
4930
4931			END DO ! cld
4932			END DO
4933			! print*,'cv3_yield fin'
4934
4935		144	RETURN
4936			END SUBROUTINE cv3_yield
4937
4938			!AC! et !RomP >>>
4939		144	SUBROUTINE cv3_tracer(nloc, len, ncum, nd, na, &
4940		144	ment, sigij, da, phi, phi2, d1a, dam, &
4941		144	ep, Vprecip, elij, clw, epmlmMm, eplaMm, &
4942			icb, inb)
4943			IMPLICIT NONE
4944
4945			include "cv3param.h"
4946
4947			!inputs:
4948			INTEGER, INTENT (IN) :: ncum, nd, na, nloc, len
4949			INTEGER, DIMENSION (len), INTENT (IN) :: icb, inb
4950			REAL, DIMENSION (len, na, na), INTENT (IN) :: ment, sigij, elij
4951			REAL, DIMENSION (len, nd), INTENT (IN) :: clw
4952			REAL, DIMENSION (len, na), INTENT (IN) :: ep
4953			REAL, DIMENSION (len, nd+1), INTENT (IN) :: Vprecip
4954			!ouputs:
4955			REAL, DIMENSION (len, na, na), INTENT (OUT) :: phi, phi2, epmlmMm
4956			REAL, DIMENSION (len, na), INTENT (OUT) :: da, d1a, dam, eplaMm
4957			!
4958			! variables pour tracer dans precip de l'AA et des mel
4959			!local variables:
4960			INTEGER i, j, k
4961		144	REAL epm(nloc, na, na)
4962
4963			! variables d'Emanuel : du second indice au troisieme
4964			! ---> tab(i,k,j) -> de l origine k a l arrivee j
4965			! ment, sigij, elij
4966			! variables personnelles : du troisieme au second indice
4967			! ---> tab(i,j,k) -> de k a j
4968			! phi, phi2
4969
4970			! initialisations
4971
4972	✓✓✓✓	5588064	da(:, :) = 0.
4973	✓✓✓✓	5588064	d1a(:, :) = 0.
4974	✓✓✓✓	5588064	dam(:, :) = 0.
4975	✓✓✓✓ ✓✓	217934640	epm(:, :, :) = 0.
4976	✓✓✓✓	5588064	eplaMm(:, :) = 0.
4977	✓✓✓✓ ✓✓	217934640	epmlmMm(:, :, :) = 0.
4978	✓✓✓✓ ✓✓	217934640	phi(:, :, :) = 0.
4979	✓✓✓✓ ✓✓	217934640	phi2(:, :, :) = 0.
4980
4981			! fraction deau condensee dans les melanges convertie en precip : epm
4982			! et eau condens�e pr�cipit�e dans masse d'air satur� : l_m*dM_m/dzdz.dzdz
4983	✓✓	4032	DO j = 1, nl
4984	✓✓	109008	DO k = 1, nl
4985	✓✓	50385078	DO i = 1, ncum
4986			IF (k>=icb(i) .AND. k<=inb(i) .AND. &
4987			!!jyg j.ge.k.and.j.le.inb(i)) then
4988			!!jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
4989	✓✓✓✓ ✓✓✓✓	50381190	j>k .AND. j<=inb(i)) THEN
4990		2960946	epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
4991			!!
4992		2960946	epm(i, j, k) = max(epm(i,j,k), 0.0)
4993			END IF
4994			END DO
4995			END DO
4996			END DO
4997
4998
4999	✓✓	4032	DO j = 1, nl
5000	✓✓	109008	DO k = 1, nl
5001	✓✓	50385078	DO i = 1, ncum
5002	✓✓✓✓	50381190	IF (k>=icb(i) .AND. k<=inb(i)) THEN
5003			eplaMm(i, j) = eplamm(i, j) + &
5004		15841089	ep(i, j)clw(i, j)ment(i, j, k)*(1.-sigij(i,j,k))
5005			END IF
5006			END DO
5007			END DO
5008			END DO
5009
5010	✓✓	4032	DO j = 1, nl
5011	✓✓	54576	DO k = 1, j - 1
5012	✓✓	24261498	DO i = 1, ncum
5013	✓✓✓✓ ✓✓	24257610	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
5014		2960946	epmlmMm(i, j, k) = epm(i, j, k)elij(i, k, j)ment(i, k, j)
5015			END IF
5016			END DO
5017			END DO
5018			END DO
5019
5020			! matrices pour calculer la tendance des concentrations dans cvltr.F90
5021	✓✓	4032	DO j = 1, nl
5022	✓✓	109008	DO k = 1, nl
5023	✓✓	50385078	DO i = 1, ncum
5024		50276214	da(i, j) = da(i, j) + (1.-sigij(i,k,j))*ment(i, k, j)
5025		50276214	phi(i, j, k) = sigij(i, k, j)*ment(i, k, j)
5026		50276214	d1a(i, j) = d1a(i, j) + ment(i, k, j)ep(i, k)(1.-sigij(i,k,j))
5027	✓✓	50381190	IF (k<=j) THEN
5028		26069148	dam(i, j) = dam(i, j) + ment(i, k, j)epm(i, k, j)(1.-ep(i,k))*(1.-sigij(i,k,j))
5029		26069148	phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
5030			END IF
5031			END DO
5032			END DO
5033			END DO
5034
5035		144	RETURN
5036			END SUBROUTINE cv3_tracer
5037			!AC! et !RomP <<<
5038
5039			SUBROUTINE cv3_uncompress(nloc, len, ncum, nd, ntra, idcum, &
5040			iflag, &
5041			precip, sig, w0, &
5042			ft, fq, fu, fv, ftra, &
5043			Ma, upwd, dnwd, dnwd0, qcondc, wd, cape, &
5044			epmax_diag, & ! epmax_cape
5045			iflag1, &
5046			precip1, sig1, w01, &
5047			ft1, fq1, fu1, fv1, ftra1, &
5048			Ma1, upwd1, dnwd1, dnwd01, qcondc1, wd1, cape1, &
5049			epmax_diag1) ! epmax_cape
5050			IMPLICIT NONE
5051
5052			include "cv3param.h"
5053
5054			!inputs:
5055			INTEGER len, ncum, nd, ntra, nloc
5056			INTEGER idcum(nloc)
5057			INTEGER iflag(nloc)
5058			REAL precip(nloc)
5059			REAL sig(nloc, nd), w0(nloc, nd)
5060			REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
5061			REAL ftra(nloc, nd, ntra)
5062			REAL ma(nloc, nd)
5063			REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
5064			REAL qcondc(nloc, nd)
5065			REAL wd(nloc), cape(nloc)
5066			REAL epmax_diag(nloc)
5067
5068			!outputs:
5069			INTEGER iflag1(len)
5070			REAL precip1(len)
5071			REAL sig1(len, nd), w01(len, nd)
5072			REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
5073			REAL ftra1(len, nd, ntra)
5074			REAL ma1(len, nd)
5075			REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
5076			REAL qcondc1(nloc, nd)
5077			REAL wd1(nloc), cape1(nloc)
5078			REAL epmax_diag1(len) ! epmax_cape
5079
5080			!local variables:
5081			INTEGER i, k, j
5082
5083			DO i = 1, ncum
5084			precip1(idcum(i)) = precip(i)
5085			iflag1(idcum(i)) = iflag(i)
5086			wd1(idcum(i)) = wd(i)
5087			cape1(idcum(i)) = cape(i)
5088			epmax_diag1(idcum(i))=epmax_diag(i) ! epmax_cape
5089			END DO
5090
5091			DO k = 1, nl
5092			DO i = 1, ncum
5093			sig1(idcum(i), k) = sig(i, k)
5094			w01(idcum(i), k) = w0(i, k)
5095			ft1(idcum(i), k) = ft(i, k)
5096			fq1(idcum(i), k) = fq(i, k)
5097			fu1(idcum(i), k) = fu(i, k)
5098			fv1(idcum(i), k) = fv(i, k)
5099			ma1(idcum(i), k) = ma(i, k)
5100			upwd1(idcum(i), k) = upwd(i, k)
5101			dnwd1(idcum(i), k) = dnwd(i, k)
5102			dnwd01(idcum(i), k) = dnwd0(i, k)
5103			qcondc1(idcum(i), k) = qcondc(i, k)
5104			END DO
5105			END DO
5106
5107			DO i = 1, ncum
5108			sig1(idcum(i), nd) = sig(i, nd)
5109			END DO
5110
5111
5112			!AC! do 2100 j=1,ntra
5113			!AC!c oct3 do 2110 k=1,nl
5114			!AC! do 2110 k=1,nd ! oct3
5115			!AC! do 2120 i=1,ncum
5116			!AC! ftra1(idcum(i),k,j)=ftra(i,k,j)
5117			!AC! 2120 continue
5118			!AC! 2110 continue
5119			!AC! 2100 continue
5120			!
5121			RETURN
5122			END SUBROUTINE cv3_uncompress
5123
5124
5125		144	subroutine cv3_epmax_fn_cape(nloc,ncum,nd &
5126			, ep,hp,icb,inb,clw,nk,t,h,hnk,lv,lf,frac &
5127		144	, pbase, p, ph, tv, buoy, sig, w0,iflag &
5128			, epmax_diag)
5129			implicit none
5130
5131			! On fait varier epmax en fn de la cape
5132			! Il faut donc recalculer ep, et hp qui a d�j� �t� calcul� et
5133			! qui en d�pend
5134			! Toutes les autres variables fn de ep sont calcul�es plus bas.
5135
5136			include "cvthermo.h"
5137			include "cv3param.h"
5138			include "conema3.h"
5139			include "cvflag.h"
5140
5141			! inputs:
5142			INTEGER, INTENT (IN) :: ncum, nd, nloc
5143			INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
5144			REAL, DIMENSION (nloc), INTENT (IN) :: hnk,pbase
5145			REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, lv, lf, tv, h
5146			REAL, DIMENSION (nloc, nd), INTENT (IN) :: clw, buoy,frac
5147			REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig,w0
5148			INTEGER, DIMENSION (nloc), INTENT (IN) :: iflag(nloc)
5149			REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
5150			REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
5151			! inouts:
5152			REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: ep,hp
5153			! outputs
5154			REAL, DIMENSION (nloc), INTENT (OUT) :: epmax_diag
5155
5156			! local
5157			integer i,k
5158			! real hp_bak(nloc,nd)
5159			! real ep_bak(nloc,nd)
5160		288	real m_loc(nloc,nd)
5161		288	real sig_loc(nloc,nd)
5162		288	real w0_loc(nloc,nd)
5163		288	integer iflag_loc(nloc)
5164		288	real cape(nloc)
5165
5166	✗✓	144	if (coef_epmax_cape.gt.1e-12) then
5167
5168			! il faut calculer la cape: on fait un calcule simple car tant qu'on ne
5169			! connait pas ep, on ne connait pas les m�langes, ddfts etc... qui sont
5170			! necessaires au calcul de la cape dans la nouvelle physique
5171
5172			! write(,) 'cv3_routines check 4303'
5173			do i=1,ncum
5174			do k=1,nd
5175			sig_loc(i,k)=sig(i,k)
5176			w0_loc(i,k)=w0(i,k)
5177			iflag_loc(i)=iflag(i)
5178			! ep_bak(i,k)=ep(i,k)
5179			enddo ! do k=1,nd
5180			enddo !do i=1,ncum
5181
5182			! write(,) 'cv3_routines check 4311'
5183			! write(,) 'nl=',nl
5184			CALL cv3_closure(nloc, ncum, nd, icb, inb, & ! na->nd
5185			pbase, p, ph, tv, buoy, &
5186			sig_loc, w0_loc, cape, m_loc,iflag_loc)
5187
5188			! write(,) 'cv3_routines check 4316'
5189			! write(,) 'ep(1,:)=',ep(1,:)
5190			do i=1,ncum
5191			epmax_diag(i)=epmax-coef_epmax_cape*sqrt(cape(i))
5192			epmax_diag(i)=amax1(epmax_diag(i),0.0)
5193			! write(,) 'i,icb,inb,cape,epmax_diag=', &
5194			! i,icb(i),inb(i),cape(i),epmax_diag(i)
5195			do k=1,nl
5196			ep(i,k)=ep(i,k)/epmax*epmax_diag(i)
5197			ep(i,k)=amax1(ep(i,k),0.0)
5198			ep(i,k)=amin1(ep(i,k),epmax_diag(i))
5199			enddo
5200			enddo
5201			! write(,) 'ep(1,:)=',ep(1,:)
5202
5203			!write(,) 'cv3_routines check 4326'
5204			! On recalcule hp:
5205			! do k=1,nl
5206			! do i=1,ncum
5207			! hp_bak(i,k)=hp(i,k)
5208			! enddo
5209			! enddo
5210			do k=1,nl
5211			do i=1,ncum
5212			hp(i,k)=h(i,k)
5213			enddo
5214			enddo
5215
5216			IF (cvflag_ice) THEN
5217
5218			do k=minorig+1,nl
5219			do i=1,ncum
5220			if((k.ge.icb(i)).and.(k.le.inb(i)))then
5221			hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k)+frac(i,k)lf(i,k))* &
5222			ep(i, k)*clw(i, k)
5223			endif
5224			enddo
5225			enddo !do k=minorig+1,n
5226			ELSE !IF (cvflag_ice) THEN
5227
5228			DO k = minorig + 1, nl
5229			DO i = 1, ncum
5230			IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
5231			hp(i,k)=hnk(i)+(lv(i,k)+(cpd-cpv)t(i,k))ep(i,k)*clw(i,k)
5232			endif
5233			enddo
5234			enddo !do k=minorig+1,n
5235
5236			ENDIF !IF (cvflag_ice) THEN
5237			!write(,) 'cv3_routines check 4345'
5238			! do i=1,ncum
5239			! do k=1,nl
5240			! if ((abs(hp_bak(i,k)-hp(i,k))/hp_bak(i,k).gt.1e-1).or. &
5241			! ((abs(hp_bak(i,k)-hp(i,k))/hp_bak(i,k).gt.1e-4).and. &
5242			! (ep(i,k)-ep_bak(i,k).lt.1e-4))) then
5243			! write(,) 'i,k=',i,k
5244			! write(,) 'coef_epmax_cape=',coef_epmax_cape
5245			! write(,) 'epmax_diag(i)=',epmax_diag(i)
5246			! write(,) 'ep(i,k)=',ep(i,k)
5247			! write(,) 'ep_bak(i,k)=',ep_bak(i,k)
5248			! write(,) 'hp(i,k)=',hp(i,k)
5249			! write(,) 'hp_bak(i,k)=',hp_bak(i,k)
5250			! write(,) 'h(i,k)=',h(i,k)
5251			! write(,) 'nk(i)=',nk(i)
5252			! write(,) 'h(i,nk(i))=',h(i,nk(i))
5253			! write(,) 'lv(i,k)=',lv(i,k)
5254			! write(,) 't(i,k)=',t(i,k)
5255			! write(,) 'clw(i,k)=',clw(i,k)
5256			! write(,) 'cpd,cpv=',cpd,cpv
5257			! stop
5258			! endif
5259			! enddo !do k=1,nl
5260			! enddo !do i=1,ncum
5261			endif !if (coef_epmax_cape.gt.1e-12) then
5262			!write(,) 'cv3_routines check 4367'
5263
5264		144	return
5265			end subroutine cv3_epmax_fn_cape
5266
5267
5268


Generated by: GCOVR (Version 4.2)