Context Navigation

cv3p2_closure.f90 @ 5319

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

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

Download in other formats:

Original Format