Context Navigation

cv3p2_closure.F90 @ 5097

Last change on this file since 5097 was 5082, checked in by abarral, 4 months ago
(lint) Fix obsolete boolean operators
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:
[3671]	43	INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag
[2374]	44	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: sig, w0
	45	REAL, DIMENSION (nloc), INTENT (INOUT) :: ptop2
	46
	47	! output:
	48	REAL, DIMENSION (nloc), INTENT (OUT) :: cape, cin
	49	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: m
	50	REAL, DIMENSION (nloc), INTENT (OUT) :: plim1, plim2
	51	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: asupmax
	52	REAL, DIMENSION (nloc), INTENT (OUT) :: supmax0
	53	REAL, DIMENSION (nloc), INTENT (OUT) :: asupmaxmin
	54	REAL, DIMENSION (nloc), INTENT (OUT) :: cbmflast, plfc
	55	REAL, DIMENSION (nloc), INTENT (OUT) :: wbeff
	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
[5082]	519	IF (iflag_mix_adiab==1) THEN
[2420]	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
[3571]	611	IF (NINT(flag_wb)==0) THEN
[2374]	612	wbeff(il) = wbmax
[3571]	613	ELSE IF (NINT(flag_wb)==1) THEN
[2374]	614	wbeff(il) = wbmax/(1.+500./(ph(il,1)-plfc(il)))
[3571]	615	ELSE IF (NINT(flag_wb)==2) THEN
[2374]	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
[5082]	627	if ((plfc(il)<ph(il,k)).and.(plfc(il)>=ph(il,k+1))) then
[2374]	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"
[5082]	640	if ((coef_clos_ls>0.).and.(plfc(il)>0.)) then
[2374]	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
[5082]	711	IF (iflag_mix_adiab==1) THEN
[2420]	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
[5082]	741	IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k)>0.)) THEN
[2420]	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
[5082]	822	IF (dtminmax(il) <= 0.) THEN
[2398]	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