Context Navigation

cv3p2_closure.F90 @ 4122

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

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

Download in other formats:

Original Format