Context Navigation

cv3p1_closure.F90 @ 4100

Last change on this file since 4100 was 2826, checked in by fhourdin, 8 years ago
Nouvelle option (provisoire) de contrôle de la convection par wbmax
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 21.9 KB

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

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

Download in other formats:

Original Format