Context Navigation

cv3p2_closure.F90 @ 3513

Last change on this file since 3513 was 3496, checked in by jyg, 6 years ago

Implementation of the ejection of liquid precipitation from the adiabatic ascents.
New flags:
+cvflag_prec_eject: logical

n -> old code, y -> new code

+ejectliq: real; possible values 0. & 1.

-> no liquid precipitation is ejected
-> all liquid precipitation is ejected

+ejectice: real; any value between 0. and 1.

fraction of solid precipitation ejected at each level

Note that the adiabatic ascent mass flux decrease due to precipitation ejection is not taken into account.

Attempts to do it led to water conservation violation.

File size: 26.3 KB

Rev	Line
[2374]	1
	2
	3	SUBROUTINE cv3p2_closure(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, &
	4	tvp, buoy, supmax, ok_inhib, ale, alp, omega,sig, w0, ptop2, cape, cin, m, &
	5	iflag, coef, plim1, plim2, asupmax, supmax0, asupmaxmin, cbmflast, plfc, &
	6	wbeff)
	7
	8
	9	! **************************************************************
	10	! *
	11	! CV3P2_CLOSURE *
	12	! Ale & Alp Closure of Convect3 *
	13	! *
	14	! written by : Kerry Emanuel *
	15	! vectorization: S. Bony *
	16	! modified by : Jean-Yves Grandpeix, 18/06/2003, 19.32.10 *
	17	! Julie Frohwirth, 14/10/2005 17.44.22 *
	18	! **************************************************************
	19
	20	USE print_control_mod, ONLY: prt_level, lunout
	21	IMPLICIT NONE
	22
	23	include "cvthermo.h"
	24	include "cv3param.h"
[3496]	25	include "cvflag.h"
[2374]	26	include "YOMCST2.h"
	27	include "YOMCST.h"
	28	include "conema3.h"
	29
	30	! input:
	31	INTEGER, INTENT (IN) :: ncum, nd, nloc
	32	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
	33	REAL, DIMENSION (nloc), INTENT (IN) :: pbase, plcl
	34	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
	35	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
	36	REAL, DIMENSION (nloc, nd), INTENT (IN) :: tv, tvp, buoy
	37	REAL, DIMENSION (nloc, nd), INTENT (IN) :: supmax
	38	LOGICAL, INTENT (IN) :: ok_inhib ! enable convection inhibition by dryness
	39	REAL, DIMENSION (nloc), INTENT (IN) :: ale, alp
	40	REAL, DIMENSION (nloc, nd), INTENT (IN) :: omega
	41
	42	! input/output:
	43	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: sig, w0
	44	REAL, DIMENSION (nloc), INTENT (INOUT) :: ptop2
	45
	46	! output:
	47	REAL, DIMENSION (nloc), INTENT (OUT) :: cape, cin
	48	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: m
	49	REAL, DIMENSION (nloc), INTENT (OUT) :: plim1, plim2
	50	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: asupmax
	51	REAL, DIMENSION (nloc), INTENT (OUT) :: supmax0
	52	REAL, DIMENSION (nloc), INTENT (OUT) :: asupmaxmin
	53	REAL, DIMENSION (nloc), INTENT (OUT) :: cbmflast, plfc
	54	REAL, DIMENSION (nloc), INTENT (OUT) :: wbeff
	55	INTEGER, DIMENSION (nloc), INTENT (OUT) :: iflag
	56
	57	! local variables:
	58	INTEGER :: il, i, j, k, icbmax
	59	INTEGER, DIMENSION (nloc) :: i0, klfc
	60	REAL :: deltap, fac, w, amu
	61	REAL, DIMENSION (nloc, nd) :: rhodp ! Factor such that m=rhodpsigw
[2420]	62	REAL :: dz
[2374]	63	REAL :: pbmxup
	64	REAL, DIMENSION (nloc, nd) :: dtmin, sigold
	65	REAL, DIMENSION (nloc, nd) :: coefmix
[2398]	66	REAL, DIMENSION (nloc) :: dtminmax
[2374]	67	REAL, DIMENSION (nloc) :: pzero, ptop2old
	68	REAL, DIMENSION (nloc) :: cina, cinb
	69	INTEGER, DIMENSION (nloc) :: ibeg
	70	INTEGER, DIMENSION (nloc) :: nsupmax
	71	REAL :: supcrit
	72	REAL, DIMENSION (nloc, nd) :: temp
	73	REAL, DIMENSION (nloc) :: p1, pmin
	74	REAL, DIMENSION (nloc) :: asupmax0
	75	LOGICAL, DIMENSION (nloc) :: ok
	76	REAL, DIMENSION (nloc, nd) :: siglim, wlim, mlim
	77	REAL, DIMENSION (nloc) :: wb2
	78	REAL, DIMENSION (nloc) :: cbmf0 ! initial cloud base mass flux
	79	REAL, DIMENSION (nloc) :: cbmflim ! cbmf given by Cape closure
	80	REAL, DIMENSION (nloc) :: cbmfalp ! cbmf given by Alp closure
	81	REAL, DIMENSION (nloc) :: cbmfalpb ! bounded cbmf given by Alp closure
	82	REAL, DIMENSION (nloc) :: cbmfmax ! upper bound on cbmf
	83	REAL, DIMENSION (nloc) :: coef
	84	REAL, DIMENSION (nloc) :: xp, xq, xr, discr, b3, b4
	85	REAL, DIMENSION (nloc) :: theta, bb
	86	REAL :: term1, term2, term3
	87	REAL, DIMENSION (nloc) :: alp2 ! Alp with offset
	88
[2420]	89	!CR: variables for new erosion of adiabiatic ascent
	90	REAL, DIMENSION (nloc, nd) :: mad, me, betalim, beta_coef
	91	REAL, DIMENSION (nloc, nd) :: med, md
[2458]	92	!jyg<
	93	! coef_peel is now in the common cv3_param
	94	!! REAL :: coef_peel
	95	!! PARAMETER (coef_peel=0.25)
	96	!>jyg
[2420]	97
[2374]	98	REAL :: sigmax
	99	PARAMETER (sigmax=0.1)
	100	!! PARAMETER (sigmax=10.)
	101
	102	CHARACTER (LEN=20) :: modname = 'cv3p2_closure'
	103	CHARACTER (LEN=80) :: abort_message
	104
	105	INTEGER,SAVE :: igout=1
	106	!$OMP THREADPRIVATE(igout)
	107
	108	IF (prt_level>=20) print *,' -> cv3p2_closure, Ale ',ale(igout)
	109
	110
	111	! -------------------------------------------------------
	112	! -- Initialization
	113	! -------------------------------------------------------
	114
	115
	116	DO il = 1, ncum
	117	alp2(il) = max(alp(il), 1.E-5)
	118	! IM
	119	alp2(il) = max(alp(il), 1.E-12)
	120	END DO
	121
	122	pbmxup = 50. ! PBMXUP+PBCRIT = cloud depth above which mixed updraughts
	123	! exist (if any)
	124
	125	IF (prt_level>=20) PRINT *, 'cv3p2_closure nloc ncum nd icb inb nl', nloc, &
	126	ncum, nd, icb(nloc), inb(nloc), nl
	127	DO k = 1, nl
	128	DO il = 1, ncum
	129	rhodp(il,k) = 0.007p(il, k)(ph(il,k)-ph(il,k+1))/tv(il, k)
	130	END DO
	131	END DO
	132
[2502]	133	!CR+jyg: initializations (up to nd) for erosion of adiabatic ascent and of m and wlim
	134	DO k = 1,nd
[2420]	135	DO il = 1, ncum
	136	mad(il,k)=0.
	137	me(il,k)=0.
	138	betalim(il,k)=1.
	139	wlim(il,k)=0.
[2502]	140	m(il, k) = 0.0
[2420]	141	ENDDO
	142	ENDDO
	143
[2374]	144	! -------------------------------------------------------
	145	! -- Reset sig(i) and w0(i) for i>inb and i<icb
	146	! -------------------------------------------------------
	147
	148	! update sig and w0 above LNB:
	149
	150	DO k = 1, nl - 1
	151	DO il = 1, ncum
	152	IF ((inb(il)<(nl-1)) .AND. (k>=(inb(il)+1))) THEN
	153	sig(il, k) = betasig(il, k) + 2.alphabuoy(il, inb(il))abs(buoy(il,inb(il)))
	154	sig(il, k) = amax1(sig(il,k), 0.0)
	155	w0(il, k) = beta*w0(il, k)
	156	END IF
	157	END DO
	158	END DO
	159
	160	! if(prt.level.GE.20) print*,'cv3p2_closure apres 100'
	161	! compute icbmax:
	162
	163	icbmax = 2
	164	DO il = 1, ncum
	165	icbmax = max(icbmax, icb(il))
	166	END DO
	167	! if(prt.level.GE.20) print*,'cv3p2_closure apres 200'
	168
	169	! update sig and w0 below cloud base:
	170
	171	DO k = 1, icbmax
	172	DO il = 1, ncum
	173	IF (k<=icb(il)) THEN
	174	sig(il, k) = betasig(il, k) - 2.alphabuoy(il, icb(il))buoy(il,icb(il))
	175	sig(il, k) = amax1(sig(il,k), 0.0)
	176	w0(il, k) = beta*w0(il, k)
	177	END IF
	178	END DO
	179	END DO
	180	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 300'
	181
	182	! -------------------------------------------------------------
	183	! -- Reset fractional areas of updrafts and w0 at initial time
	184	! -- and after 10 time steps of no convection
	185	! -------------------------------------------------------------
	186
[2398]	187	!jyg<
	188	IF (ok_convstop) THEN
	189	DO k = 1, nl - 1
	190	DO il = 1, ncum
	191	IF (sig(il,nd)<1.5 .OR. sig(il,nd)>noconv_stop) THEN
	192	sig(il, k) = 0.0
	193	w0(il, k) = 0.0
	194	END IF
	195	END DO
	196	END DO
	197	ELSE
[2374]	198	DO k = 1, nl - 1
	199	DO il = 1, ncum
	200	IF (sig(il,nd)<1.5 .OR. sig(il,nd)>12.0) THEN
	201	sig(il, k) = 0.0
	202	w0(il, k) = 0.0
	203	END IF
	204	END DO
	205	END DO
[2398]	206	ENDIF ! (ok_convstop)
	207	!>jyg
[2374]	208	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 400'
	209
	210	! -------------------------------------------------------
	211	! -- Compute initial cloud base mass flux (Cbmf0)
	212	! -------------------------------------------------------
	213	DO il = 1, ncum
	214	cbmf0(il) = 0.0
	215	END DO
	216
	217	DO k = 1, nl
	218	DO il = 1, ncum
	219	IF (k>=icb(il) .AND. k<=inb(il) &
	220	.AND. icb(il)+1<=inb(il)) THEN
	221	cbmf0(il) = cbmf0(il) + sig(il, k)w0(il,k)rhodp(il,k)
	222	END IF
	223	END DO
	224	END DO
	225
	226	! -------------------------------------------------------------
	227	! jyg1
	228	! -- Calculate adiabatic ascent top pressure (ptop)
	229	! -------------------------------------------------------------
	230
	231
	232	! c 1. Start at first level where precipitations form
	233	DO il = 1, ncum
	234	pzero(il) = plcl(il) - pbcrit
	235	END DO
	236
	237	! c 2. Add offset
	238	DO il = 1, ncum
	239	pzero(il) = pzero(il) - pbmxup
	240	END DO
	241	DO il = 1, ncum
	242	ptop2old(il) = ptop2(il)
	243	END DO
	244
	245	DO il = 1, ncum
	246	! CR:c est quoi ce 300??
	247	p1(il) = pzero(il) - 300.
	248	END DO
	249
	250	! compute asupmax=abs(supmax) up to lnm+1
	251
	252	DO il = 1, ncum
	253	ok(il) = .TRUE.
	254	nsupmax(il) = inb(il)
	255	END DO
	256
	257	DO i = 1, nl
	258	DO il = 1, ncum
	259	IF (i>icb(il) .AND. i<=inb(il)) THEN
	260	IF (p(il,i)<=pzero(il) .AND. supmax(il,i)<0 .AND. ok(il)) THEN
	261	nsupmax(il) = i
	262	ok(il) = .FALSE.
	263	END IF ! end IF (P(i) ... )
	264	END IF ! end IF (icb+1 le i le inb)
	265	END DO
	266	END DO
	267
	268	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 2.'
	269	DO i = 1, nl
	270	DO il = 1, ncum
	271	asupmax(il, i) = abs(supmax(il,i))
	272	END DO
	273	END DO
	274
	275
	276	DO il = 1, ncum
	277	asupmaxmin(il) = 10.
	278	pmin(il) = 100.
	279	! IM ??
	280	asupmax0(il) = 0.
	281	END DO
	282
	283	! c 3. Compute in which level is Pzero
	284
	285	! IM bug i0 = 18
	286	DO il = 1, ncum
	287	i0(il) = nl
	288	END DO
	289
	290	DO i = 1, nl
	291	DO il = 1, ncum
	292	IF (i>icb(il) .AND. i<=inb(il)) THEN
	293	IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
	294	IF (pzero(il)>p(il,i) .AND. pzero(il)<p(il,i-1)) THEN
	295	i0(il) = i
	296	END IF
	297	END IF
	298	END IF
	299	END DO
	300	END DO
	301	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 3.'
	302
	303	! c 4. Compute asupmax at Pzero
	304
	305	DO i = 1, nl
	306	DO il = 1, ncum
	307	IF (i>icb(il) .AND. i<=inb(il)) THEN
	308	IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
	309	asupmax0(il) = ((pzero(il)-p(il,i0(il)-1))*asupmax(il,i0(il))- &
	310	(pzero(il)-p(il,i0(il)))*asupmax(il,i0(il)-1))/(p(il,i0(il))-p(il,i0(il)-1))
	311	END IF
	312	END IF
	313	END DO
	314	END DO
	315
	316
	317	DO i = 1, nl
	318	DO il = 1, ncum
	319	IF (p(il,i)==pzero(il)) THEN
	320	asupmax(i, il) = asupmax0(il)
	321	END IF
	322	END DO
	323	END DO
	324	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 4.'
	325
	326	! c 5. Compute asupmaxmin, minimum of asupmax
	327
	328	DO i = 1, nl
	329	DO il = 1, ncum
	330	IF (i>icb(il) .AND. i<=inb(il)) THEN
	331	IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
	332	IF (asupmax(il,i)<asupmaxmin(il)) THEN
	333	asupmaxmin(il) = asupmax(il, i)
	334	pmin(il) = p(il, i)
	335	END IF
	336	END IF
	337	END IF
	338	END DO
	339	END DO
	340
	341	DO il = 1, ncum
	342	! IM
	343	IF (prt_level>=20) THEN
	344	PRINT *, 'cv3p2_closure il asupmax0 asupmaxmin', il, asupmax0(il), &
	345	asupmaxmin(il), pzero(il), pmin(il)
	346	END IF
	347	IF (asupmax0(il)<asupmaxmin(il)) THEN
	348	asupmaxmin(il) = asupmax0(il)
	349	pmin(il) = pzero(il)
	350	END IF
	351	END DO
	352	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 5.'
	353
	354
	355	! Compute Supmax at Pzero
	356
	357	DO i = 1, nl
	358	DO il = 1, ncum
	359	IF (i>icb(il) .AND. i<=inb(il)) THEN
	360	IF (p(il,i)<=pzero(il)) THEN
	361	supmax0(il) = ((p(il,i)-pzero(il))*asupmax(il,i-1)- &
	362	(p(il,i-1)-pzero(il))*asupmax(il,i))/(p(il,i)-p(il,i-1))
	363	GO TO 425
	364	END IF ! end IF (P(i) ... )
	365	END IF ! end IF (icb+1 le i le inb)
	366	END DO
	367	END DO
	368
	369	425 CONTINUE
	370	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 425.'
	371
	372	! c 6. Calculate ptop2
	373
	374	DO il = 1, ncum
	375	IF (asupmaxmin(il)<supcrit1) THEN
	376	ptop2(il) = pmin(il)
	377	END IF
	378
	379	IF (asupmaxmin(il)>supcrit1 .AND. asupmaxmin(il)<supcrit2) THEN
	380	ptop2(il) = ptop2old(il)
	381	END IF
	382
	383	IF (asupmaxmin(il)>supcrit2) THEN
	384	ptop2(il) = ph(il, inb(il))
	385	END IF
	386	END DO
	387
	388	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 6.'
	389
	390	! c 7. Compute multiplying factor for adiabatic updraught mass flux
	391
	392
	393	IF (ok_inhib) THEN
	394
	395	DO i = 1, nl
	396	DO il = 1, ncum
	397	IF (i<=nl) THEN
	398	coefmix(il, i) = (min(ptop2(il),ph(il,i))-ph(il,i))/(ph(il,i+1)-ph(il,i))
	399	coefmix(il, i) = min(coefmix(il,i), 1.)
	400	END IF
	401	END DO
	402	END DO
	403
	404
	405	ELSE ! when inhibition is not taken into account, coefmix=1
	406
	407
	408
	409	DO i = 1, nl
	410	DO il = 1, ncum
	411	IF (i<=nl) THEN
	412	coefmix(il, i) = 1.
	413	END IF
	414	END DO
	415	END DO
	416
	417	END IF ! ok_inhib
	418	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 7.'
	419	! -------------------------------------------------------------------
	420	! -------------------------------------------------------------------
	421
	422
	423	! jyg2
	424
	425	! ==========================================================================
	426
	427
	428	! -------------------------------------------------------------
	429	! -- Calculate convective inhibition (CIN)
	430	! -------------------------------------------------------------
	431
	432	! do i=1,nloc
	433	! print*,'avant cine p',pbase(i),plcl(i)
	434	! enddo
	435	! do j=1,nd
	436	! do i=1,nloc
	437	! print*,'avant cine t',tv(i),tvp(i)
	438	! enddo
	439	! enddo
	440	CALL cv3_cine(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, tvp, cina, &
	441	cinb, plfc)
	442
	443	DO il = 1, ncum
	444	cin(il) = cina(il) + cinb(il)
	445	END DO
	446	IF (prt_level>=20) PRINT *, 'cv3p2_closure after cv3_cine: cina, cinb, cin ', &
	447	cina(igout), cinb(igout), cin(igout)
	448	! -------------------------------------------------------------
	449	! --Update buoyancies to account for Ale
	450	! -------------------------------------------------------------
	451
	452	CALL cv3_buoy(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, ale, cin, tv, &
	453	tvp, buoy)
	454	IF (prt_level>=20) PRINT *, 'cv3p2_closure after cv3_buoy'
	455
	456	! -------------------------------------------------------------
	457	! -- Calculate convective available potential energy (cape),
	458	! -- vertical velocity (w), fractional area covered by
	459	! -- undilute updraft (sig), and updraft mass flux (m)
	460	! -------------------------------------------------------------
	461
	462	DO il = 1, ncum
	463	cape(il) = 0.0
[2398]	464	dtminmax(il) = -100.
[2374]	465	END DO
	466
	467	! compute dtmin (minimum buoyancy between ICB and given level k):
	468
	469	DO k = 1, nl
	470	DO il = 1, ncum
	471	dtmin(il, k) = 100.0
	472	END DO
	473	END DO
	474
	475	DO k = 1, nl
	476	DO j = minorig, nl
	477	DO il = 1, ncum
	478	IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (j>=icb(il)) &
	479	.AND. (j<=(k-1))) THEN
	480	dtmin(il, k) = amin1(dtmin(il,k), buoy(il,j))
	481	END IF
	482	END DO
	483	END DO
	484	END DO
[2398]	485	!jyg<
	486	! Store maximum of dtmin
	487	! C est pas terrible d avoir ce test sur Ale+Cin encore une fois ici.
	488	! A REVOIR !
	489	DO k = 1, nl
	490	DO il = 1, ncum
	491	IF (k>=(icb(il)+1) .AND. k<=inb(il) .AND. ale(il)+cin(il)>0.) THEN
	492	dtminmax(il) = max(dtmin(il,k), dtminmax(il))
	493	ENDIF
	494	END DO
	495	END DO
[2374]	496	!
[2398]	497	! prevent convection when ale+cin <= 0
	498	DO k = 1, nl
	499	DO il = 1, ncum
	500	IF (k>=(icb(il)+1) .AND. k<=inb(il)) THEN
	501	dtmin(il,k) = min(dtmin(il,k), dtminmax(il))
	502	ENDIF
	503	END DO
	504	END DO
	505	!>jyg
	506	!
[2374]	507	IF (prt_level >= 20) THEN
	508	print *,'cv3p2_closure: dtmin ', (k, dtmin(igout,k), k=1,nl)
[2398]	509	print *,'cv3p2_closure: dtminmax ', dtminmax(igout)
[2374]	510	ENDIF
	511	!
	512	! the interval on which cape is computed starts at pbase :
	513
	514	DO k = 1, nl
	515	DO il = 1, ncum
	516
	517	IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
	518
[2420]	519	IF (iflag_mix_adiab.eq.1) THEN
	520	!CR:computation of cape from LCL: keep flag or to modify in all cases?
	521	deltap = min(plcl(il), ph(il,k-1)) - min(plcl(il), ph(il,k))
	522	ELSE
[2374]	523	deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
[2420]	524	ENDIF
[2374]	525	cape(il) = cape(il) + rrdbuoy(il, k-1)deltap/p(il, k-1)
	526	cape(il) = amax1(0.0, cape(il))
	527	sigold(il, k) = sig(il, k)
	528
	529
	530	! jyg Coefficient coefmix limits convection to levels where a
	531	! sufficient
	532	! fraction of mixed draughts are ascending.
	533	siglim(il, k) = coefmix(il, k)alpha1dtmin(il, k)*abs(dtmin(il,k))
	534	siglim(il, k) = amax1(siglim(il,k), 0.0)
	535	siglim(il, k) = amin1(siglim(il,k), 0.01)
	536	! c fac=AMIN1(((dtcrit-dtmin(il,k))/dtcrit),1.0)
	537	fac = 1.
	538	wlim(il, k) = fac*sqrt(cape(il))
	539	amu = siglim(il, k)*wlim(il, k)
	540	!! rhodp(il,k) = 0.007p(il, k)(ph(il,k)-ph(il,k+1))/tv(il, k) !cor jyg : computed earlier
	541	mlim(il, k) = amu*rhodp(il,k)
	542	! print*, 'siglim ', k,siglim(1,k)
	543	END IF
	544
	545	END DO
	546	END DO
	547	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 600'
	548
	549	DO il = 1, ncum
	550	! IM beg
	551	IF (prt_level>=20) THEN
	552	PRINT *, 'cv3p2_closure il icb mlim ph ph+1 ph+2', il, icb(il), &
	553	mlim(il, icb(il)+1), ph(il, icb(il)), ph(il, icb(il)+1), &
	554	ph(il, icb(il)+2)
	555	END IF
	556
	557	IF (icb(il)+1<=inb(il)) THEN
	558	! IM end
	559	mlim(il, icb(il)) = 0.5mlim(il,icb(il)+1)(ph(il,icb(il))-ph(il,icb(il)+1))/ &
	560	(ph(il,icb(il)+1)-ph(il,icb(il)+2))
	561	! IM beg
	562	END IF !(icb(il.le.inb(il))) then
	563	! IM end
	564	END DO
	565	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 700'
	566
	567	!
	568	! ------------------------------------------------------------------------
	569	! c Compute Cloud base mass flux given by Cape closure (cbmflim = cbmf of
	570	! c elementary systems), cbmf given by Alp closure (cbmfalp), cbmf given by Alp
	571	! c closure with an upper bound imposed (cbmfalpb) and cbmf resulting from
	572	! c time integration (cbmflast).
	573	! ------------------------------------------------------------------------
	574
	575	DO il = 1, ncum
	576	cbmflim(il) = 0.
	577	cbmfalp(il) = 0.
	578	cbmfalpb(il) = 0.
	579	cbmflast(il) = 0.
	580	END DO
	581
	582	! c 1. Compute cloud base mass flux of elementary system (Cbmflim)
	583
	584	DO k = 1, nl
	585	DO il = 1, ncum
	586	! old IF (k .ge. icb(il) .and. k .le. inb(il)) THEN
	587	! IM IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN
	588	IF (k>=icb(il) .AND. k<=inb(il) & !cor jyg
	589	.AND. icb(il)+1<=inb(il)) THEN !cor jyg
	590	cbmflim(il) = cbmflim(il) + mlim(il, k)
	591	END IF
	592	END DO
	593	END DO
	594	IF (prt_level>=20) PRINT *, 'cv3p2_closure after cbmflim: cbmflim ', cbmflim(igout)
	595
	596	! 1.5 Compute cloud base mass flux given by Alp closure (Cbmfalp), maximum
	597	! allowed mass flux (Cbmfmax) and bounded mass flux (Cbmfalpb)
	598	! Cbmfalpb is set to zero if Cbmflim (the mass flux of elementary cloud)
	599	! is exceedingly small.
	600
	601	DO il = 1, ncum
	602	wb2(il) = sqrt(2.*max(ale(il)+cin(il),0.))
	603	END DO
	604
	605	DO il = 1, ncum
	606	IF (plfc(il)<100.) THEN
	607	! This is an irealistic value for plfc => no calculation of wbeff
	608	wbeff(il) = 100.1
	609	ELSE
	610	! Calculate wbeff
	611	IF (flag_wb==0) THEN
	612	wbeff(il) = wbmax
	613	ELSE IF (flag_wb==1) THEN
	614	wbeff(il) = wbmax/(1.+500./(ph(il,1)-plfc(il)))
	615	ELSE IF (flag_wb==2) THEN
	616	wbeff(il) = wbmax(0.01(ph(il,1)-plfc(il)))**2
	617	END IF
	618	END IF
	619	END DO
	620
	621	!CR:Compute k at plfc
	622	DO il=1,ncum
	623	klfc(il)=nl
	624	ENDDO
	625	DO k=1,nl
	626	DO il=1,ncum
	627	if ((plfc(il).lt.ph(il,k)).and.(plfc(il).ge.ph(il,k+1))) then
	628	klfc(il)=k
	629	endif
	630	ENDDO
	631	ENDDO
	632	!RC
	633
	634	DO il = 1, ncum
	635	! jyg Modification du coef de wb*wb pour conformite avec papier Wake
	636	! c cbmfalp(il) = alp2(il)/(0.5wbwb-Cin(il))
	637	cbmfalp(il) = alp2(il)/(2.wbeff(il)wbeff(il)-cin(il))
	638	!CR: Add large-scale component to the mass-flux
	639	!encore connu sous le nom "Experience du tube de dentifrice"
	640	if ((coef_clos_ls.gt.0.).and.(plfc(il).gt.0.)) then
	641	cbmfalp(il) = cbmfalp(il) - coef_clos_lsmin(0.,1./RGomega(il,klfc(il)))
	642	endif
	643	!RC
	644	IF (cbmfalp(il)==0 .AND. alp2(il)/=0.) THEN
	645	WRITE (lunout, *) 'cv3p2_closure cbmfalp=0 and alp NE 0 il alp2 alp cin ' , &
	646	il, alp2(il), alp(il), cin(il)
	647	abort_message = ''
	648	CALL abort_physic(modname, abort_message, 1)
	649	END IF
	650	cbmfmax(il) = sigmaxwb2(il)100.p(il, icb(il))/(rrdtv(il,icb(il)))
	651	END DO
	652
[2398]	653	!jyg<
	654	IF (OK_intermittent) THEN
	655	DO il = 1, ncum
	656	IF (cbmflim(il)>1.E-6) THEN
	657	cbmfalpb(il) = min(cbmfalp(il), (cbmfmax(il)-beta*cbmf0(il))/(1.-beta))
	658	! print*,'cbmfalpb',cbmfalpb(il),cbmfmax(il)
	659	END IF
	660	END DO
	661	ELSE
	662	!>jyg
[2374]	663	DO il = 1, ncum
	664	IF (cbmflim(il)>1.E-6) THEN
	665	! ATTENTION TEST CR
	666	! if (cbmfmax(il).lt.1.e-12) then
	667	cbmfalpb(il) = min(cbmfalp(il), cbmfmax(il))
	668	! else
	669	! cbmfalpb(il) = cbmfalp(il)
	670	! endif
	671	! print*,'cbmfalpb',cbmfalp(il),cbmfmax(il)
	672	END IF
	673	END DO
[2398]	674	ENDIF !(OK_intermittent)
[2374]	675	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres cbmfalpb: cbmfalpb ',cbmfalpb(igout)
	676
	677	! c 2. Compute coefficient and apply correction
	678
	679	DO il = 1, ncum
	680	coef(il) = (cbmfalpb(il)+1.E-10)/(cbmflim(il)+1.E-10)
	681	END DO
	682	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres coef_plantePLUS'
	683
	684	DO k = 1, nl
	685	DO il = 1, ncum
	686	IF (k>=icb(il)+1 .AND. k<=inb(il)) THEN
	687	amu = betasig(il, k)w0(il, k) + (1.-beta)coef(il)siglim(il, k)*wlim(il, k)
	688	w0(il, k) = wlim(il, k)
	689	w0(il, k) = max(w0(il,k), 1.E-10)
	690	sig(il, k) = amu/w0(il, k)
	691	sig(il, k) = min(sig(il,k), 1.)
	692	! c amu = 0.5(SIG(il,k)+sigold(il,k))W0(il,k)
	693	!jyg m(il, k) = amu0.007p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
	694	m(il, k) = amu*rhodp(il,k)
	695	END IF
	696	END DO
	697	END DO
	698	! jyg2
	699	DO il = 1, ncum
	700	w0(il, icb(il)) = 0.5*w0(il, icb(il)+1)
	701	m(il, icb(il)) = 0.5m(il, icb(il)+1)(ph(il,icb(il))-ph(il,icb(il)+1))/ &
	702	(ph(il,icb(il)+1)-ph(il,icb(il)+2))
	703	sig(il, icb(il)) = sig(il, icb(il)+1)
	704	sig(il, icb(il)-1) = sig(il, icb(il))
	705	END DO
	706	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres w0_sig_M: w0, sig ', &
	707	(k,w0(igout,k),sig(igout,k), k=icb(igout),inb(igout))
	708
[2420]	709	!CR: new erosion of adiabatic ascent: modification of m
	710	!computation of the sum of ascending fluxes
	711	IF (iflag_mix_adiab.eq.1) THEN
	712
	713	!Verification sum(me)=sum(m)
[2502]	714	DO k = 1,nd
[2420]	715	DO il = 1, ncum
	716	md(il,k)=0.
	717	med(il,k)=0.
	718	ENDDO
	719	ENDDO
	720
	721	DO k = nl,1,-1
	722	DO il = 1, ncum
	723	md(il,k)=md(il,k+1)+m(il,k+1)
	724	ENDDO
	725	ENDDO
	726
	727	DO k = nl,1,-1
	728	DO il = 1, ncum
	729	IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN
	730	mad(il,k)=mad(il,k+1)+m(il,k+1)
	731	ENDIF
	732	! print*,"mad",il,k,mad(il,k)
	733	ENDDO
	734	ENDDO
	735
	736	!CR: erosion of each adiabatic ascent during its ascent
	737
	738	!Computation of erosion coefficient beta_coef
	739	DO k = 1, nl
	740	DO il = 1, ncum
	741	IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN
	742	! print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1)
	743	beta_coef(il,k)=RGcoef_peelbuoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2
	744	ELSE
	745	beta_coef(il,k)=0.
	746	ENDIF
	747	ENDDO
	748	ENDDO
	749
	750	! print*,"apres beta_coef"
	751
	752	DO k = 1, nl
	753	DO il = 1, ncum
	754
	755	IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
	756
	757	! print*,"dz",il,k,tv(il, k-1)
	758	dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrdtv(il, k-1))RG)
	759	betalim(il,k)=betalim(il,k-1)exp(-1.beta_coef(il,k-1)*dz)
	760	! betalim(il,k)=betalim(il,k-1)exp(-RGcoef_peelbuoy(il,k-1)/tv(il,k-1)/5.2dz)
	761	! print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k)
	762	dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrdtv(il, k))RG)
	763	! me(il,k)=betalim(il,k)(m(il,k)+RGcoef_peelbuoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)2dz*mad(il,k))
	764	me(il,k)=betalim(il,k)(m(il,k)+beta_coef(il,k)dz*mad(il,k))
	765	! print,"B/w2",il,k,RGcoef_peelbuoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)2dz
	766
	767	END IF
	768
	769	!Modification of m
	770	m(il,k)=me(il,k)
	771	END DO
	772	END DO
	773
	774	! DO il = 1, ncum
	775	! dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrdtv(il, icb(il)))RG)
	776	! m(il,icb(il))=m(il,icb(il))+RGcoef_peelbuoy(il,icb(il))/tv(il,icb(il)) &
	777	! /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)*2dz*mad(il,icb(il))
	778	! print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il))
	779	! ENDDO
	780
	781	!Verification sum(me)=sum(m)
	782	DO k = nl,1,-1
	783	DO il = 1, ncum
	784	med(il,k)=med(il,k+1)+m(il,k+1)
	785	! print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k)
	786	ENDDO
	787	ENDDO
	788
	789
	790	ENDIF !(iflag_mix_adiab)
	791	!RC
	792
[2374]	793	! c 3. Compute final cloud base mass flux;
	794	! c set iflag to 3 if cloud base mass flux is exceedingly small and is
	795	! c decreasing (i.e. if the final mass flux (cbmflast) is greater than
	796	! c the target mass flux (cbmfalpb)).
[2398]	797	! c If(ok_convstop): set iflag to 4 if no positive buoyancy has been met
[2374]	798
	799	!jyg DO il = 1, ncum
	800	!jyg cbmflast(il) = 0.
	801	!jyg END DO
	802
	803	DO k = 1, nl
	804	DO il = 1, ncum
	805	IF (k>=icb(il) .AND. k<=inb(il)) THEN
	806	!IMpropo?? IF ((k.ge.(icb(il)+1)).and.(k.le.inb(il))) THEN
	807	cbmflast(il) = cbmflast(il) + m(il, k)
	808	END IF
	809	END DO
	810	END DO
	811	IF (prt_level>=20) PRINT *, 'cv3p2_closure apres cbmflast: cbmflast ',cbmflast(igout)
	812
	813	DO il = 1, ncum
	814	IF (cbmflast(il)<1.E-6 .AND. cbmflast(il)>=cbmfalpb(il)) THEN
	815	iflag(il) = 3
	816	END IF
	817	END DO
	818
[2398]	819	!jyg<
	820	IF (ok_convstop) THEN
	821	DO il = 1, ncum
	822	IF (dtminmax(il) .LE. 0.) THEN
	823	iflag(il) = 4
	824	END IF
	825	END DO
	826	ELSE
	827	!>jyg
[2374]	828	DO k = 1, nl
	829	DO il = 1, ncum
	830	IF (iflag(il)>=3) THEN
	831	m(il, k) = 0.
	832	sig(il, k) = 0.
	833	w0(il, k) = 0.
	834	END IF
	835	END DO
	836	END DO
[2398]	837	ENDIF ! (ok_convstop)
[2374]	838	!
	839	IF (prt_level >= 10) THEN
	840	print *,'cv3p2_closure: iflag ',iflag(igout)
	841	ENDIF
	842	!
	843
	844	! c 4. Introduce a correcting factor for coef, in order to obtain an
	845	! effective
	846	! c sigdz larger in the present case (using cv3p2_closure) than in the
	847	! old
	848	! c closure (using cv3_closure).
	849	IF (1==0) THEN
	850	DO il = 1, ncum
	851	! c coef(il) = 2.*coef(il)
	852	coef(il) = 5.*coef(il)
	853	END DO
	854	! version CVS du ..2008
	855	ELSE
	856	IF (iflag_cvl_sigd==0) THEN
	857	! test pour verifier qu on fait la meme chose qu avant: sid constant
	858	coef(1:ncum) = 1.
	859	ELSE
	860	coef(1:ncum) = min(2.*coef(1:ncum), 5.)
	861	coef(1:ncum) = max(2.*coef(1:ncum), 0.2)
	862	END IF
	863	END IF
	864
	865	IF (prt_level>=20) PRINT *, 'cv3p2_closure FIN'
	866	RETURN
	867	END SUBROUTINE cv3p2_closure
	868
	869

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format