Context Navigation

cv3p1_closure.F90 @ 2461

Last change on this file since 2461 was 2435, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes r2396:2434 into testing branch
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.6 KB

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

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

Download in other formats:

Original Format