Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

cv3p_mixing.F90 @ 5452

Last change on this file since 5452 was 3331, checked in by acozic, 7 years ago
Add modification for isotopes
Property svn:executable set to ``*
File size: 24.2 KB

Rev	Line
[3331]	1	SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
	2	ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
	3	unk, vnk, hp, tv, tvp, ep, clw, sig, &
	4	Ment, Qent, hent, uent, vent, nent, &
	5	Sigij, elij, supmax, Ments, Qents, traent)
	6	! **************************************************************
	7	! *
	8	! CV3P_MIXING : compute mixed draught properties and, *
	9	! within a scaling factor, mixed draught *
	10	! mass fluxes. *
	11	! written by : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15*
	12	! modified by : *
	13	! **************************************************************
	14
	15	USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
	16
	17	IMPLICIT NONE
	18
	19	include "cvthermo.h"
	20	include "cv3param.h"
	21	include "YOMCST2.h"
	22	include "cvflag.h"
	23
	24	!inputs:
	25	INTEGER, INTENT (IN) :: ncum, nd, na
	26	INTEGER, INTENT (IN) :: ntra, nloc
	27	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
	28	REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig
	29	REAL, DIMENSION (nloc), INTENT (IN) :: qnk, unk, vnk
	30	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
	31	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
	32	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
	33	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra ! input of convect3
	34	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv
	35	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
	36	REAL, DIMENSION (nloc, na), INTENT (IN) :: frac !ice fraction in condensate
	37	REAL, DIMENSION (nloc, na), INTENT (IN) :: h !liquid water static energy of environment
	38	REAL, DIMENSION (nloc, na), INTENT (IN) :: hp !liquid water static energy of air shed from adiab. asc.
	39	REAL, DIMENSION (nloc, na), INTENT (IN) :: tv, tvp
	40	REAL, DIMENSION (nloc, na), INTENT (IN) :: ep, clw
	41
	42	!outputs:
	43	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: Ment, Qent
	44	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: uent, vent
	45	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: Sigij, elij
	46	REAL, DIMENSION (nloc, na), INTENT (OUT) :: supmax ! Highest mixing fraction of mixed
	47	! updraughts with the sign of (h-hp)
	48	REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
	49	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: Ments, Qents
	50	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: hent
	51	INTEGER, DIMENSION (nloc, nd), INTENT (OUT) :: nent
	52
	53	!local variables:
	54	INTEGER i, j, k, il, im, jm
	55	INTEGER num1, num2
	56	REAL :: rti, bf2, anum, denom, dei, altem, cwat, stemp
	57	REAL :: alt, delp, delm
	58	REAL, DIMENSION (nloc) :: Qmixmax, Rmixmax, sqmrmax
	59	REAL, DIMENSION (nloc) :: Qmixmin, Rmixmin, sqmrmin
	60	REAL, DIMENSION (nloc) :: signhpmh
	61	REAL, DIMENSION (nloc) :: Sx
	62	REAL :: Scrit2
	63	REAL, DIMENSION (nloc) :: Smid, Sjmin, Sjmax
	64	REAL, DIMENSION (nloc) :: Sbef, sup, smin
	65	REAL, DIMENSION (nloc) :: ASij, ASij_inv, smax, Scrit
	66	REAL, DIMENSION (nloc, nd, nd) :: Sij
	67	REAL, DIMENSION (nloc, nd) :: csum
	68	REAL :: awat
	69	REAL :: cpm !Mixed draught heat capacity
	70	REAL :: Tm !Mixed draught temperature
	71	LOGICAL, DIMENSION (nloc) :: lwork
	72
	73	REAL amxupcrit, df, ff
	74	INTEGER nstep
	75
	76	INTEGER,SAVE :: igout=1
	77	!$OMP THREADPRIVATE(igout)
	78
	79	! -- Mixing probability distribution functions
	80
	81	REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
	82
	83	Qcoef1(F) = tanh(F/gammas)
	84	Qcoef2(F) = (tanh(F/gammas)+gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
	85	QFF(F) = max(min(F,1.), 0.)
	86	QFFf(F) = min(QFF(F), scut)
	87	Qmix1(F) = (tanh((QFF(F)-Fmax)/gammas)+Qcoef1max)/Qcoef2max
	88	Rmix1(F) = (gammaslog(cosh((QFF(F)-Fmax)/gammas))+QFF(F)Qcoef1max)/Qcoef2max
	89	Qmix2(F) = -log(1.-QFFf(F))/scut
	90	Rmix2(F) = (QFFf(F)+(1.-QFF(F))*log(1.-QFFf(F)))/scut
	91	Qmix(F) = qqa1Qmix1(F) + qqa2Qmix2(F)
	92	Rmix(F) = qqa1Rmix1(F) + qqa2Rmix2(F)
	93
	94	INTEGER, SAVE :: ifrst
	95	DATA ifrst/0/
	96	!$OMP THREADPRIVATE(ifrst)
	97
	98
	99	! =====================================================================
	100	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	101	! =====================================================================
	102
	103	! -- Initialize mixing PDF coefficients
	104	IF (ifrst==0) THEN
	105	ifrst = 1
	106	Qcoef1max = Qcoef1(Fmax)
	107	Qcoef2max = Qcoef2(Fmax)
	108	!<jyg
	109	print*, 'fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max ', &
	110	fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max
	111	!>jyg
	112
	113	END IF
	114
	115
	116	! ori do 360 i=1,ncum*nlp
	117	DO j = 1, nl
	118	DO i = 1, ncum
	119	nent(i, j) = 0
	120	! in convect3, m is computed in cv3_closure
	121	! ori m(i,1)=0.0
	122	END DO
	123	END DO
	124
	125	! ori do 400 k=1,nlp
	126	! ori do 390 j=1,nlp
	127	DO j = 1, nl
	128	DO k = 1, nl
	129	DO i = 1, ncum
	130	Qent(i, k, j) = rr(i, j)
	131	uent(i, k, j) = u(i, j)
	132	vent(i, k, j) = v(i, j)
	133	elij(i, k, j) = 0.0
	134	hent(i, k, j) = 0.0
	135	!AC! Ment(i,k,j)=0.0
	136	!AC! Sij(i,k,j)=0.0
	137	END DO
	138	END DO
	139	END DO
	140
	141	!AC!
	142	Ment(1:ncum, 1:nd, 1:nd) = 0.0
	143	Sij(1:ncum, 1:nd, 1:nd) = 0.0
	144	!AC!
	145	!ym
	146	Sigij(1:ncum, 1:nd, 1:nd) = 0.0
	147	!ym
	148
	149	!jyg! DO k = 1, ntra
	150	!jyg! DO j = 1, nd ! instead nlp
	151	!jyg! DO i = 1, nd ! instead nlp
	152	!jyg! DO il = 1, ncum
	153	!jyg! traent(il, i, j, k) = tra(il, j, k)
	154	!jyg! END DO
	155	!jyg! END DO
	156	!jyg! END DO
	157	!jyg! END DO
	158
	159	! =====================================================================
	160	! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
	161	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	162	! --- FRACTION (Sij)
	163	! =====================================================================
	164
	165	DO i = minorig + 1, nl
	166
	167	DO j = minorig, nl
	168	DO il = 1, ncum
	169	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
	170	.AND. (j<=inb(il))) THEN
	171
	172	rti = qnk(il) - ep(il, i)*clw(il, i)
	173	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
	174	!jyg(from aj)<
	175	IF (cvflag_ice) THEN
	176	! print*,cvflag_ice,'cvflag_ice dans do 700'
	177	IF (t(il,j)<=263.15) THEN
	178	bf2 = 1. + (lf(il,j)+lv(il,j))(lv(il,j)+frac(il,j) &
	179	lf(il,j))rs(il, j)/(rrvt(il,j)t(il,j)cpd)
	180	END IF
	181	END IF
	182	!>jyg
	183	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
	184	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
	185	dei = denom
	186	IF (abs(dei)<0.01) dei = 0.01
	187	Sij(il, i, j) = anum/dei
	188	Sij(il, i, i) = 1.0
	189	altem = Sij(il, i, j)rr(il, i) + (1.-Sij(il,i,j))rti - rs(il, j)
	190	altem = altem/bf2
	191	cwat = clw(il, j)*(1.-ep(il,j))
	192	stemp = Sij(il, i, j)
	193	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
	194	!jyg(from aj)<
	195	IF (cvflag_ice) THEN
	196	anum = anum - (lv(il,j)+frac(il,j)lf(il,j))(rti-rs(il,j)-cwat*bf2)
	197	denom = denom + (lv(il,j)+frac(il,j)lf(il,j))(rr(il,i)-rti)
	198	ELSE
	199	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
	200	denom = denom + lv(il, j)*(rr(il,i)-rti)
	201	END IF
	202	!>jyg
	203	IF (abs(denom)<0.01) denom = 0.01
	204	Sij(il, i, j) = anum/denom
	205	altem = Sij(il, i, j)rr(il, i) + (1.-Sij(il,i,j))rti - rs(il, j)
	206	altem = altem - (bf2-1.)*cwat
	207	END IF
	208	IF (Sij(il,i,j)>0.0) THEN
	209	!!! Ment(il,i,j)=m(il,i)
	210	Ment(il, i, j) = 1.
	211	elij(il, i, j) = altem
	212	elij(il, i, j) = amax1(0.0, elij(il,i,j))
	213	nent(il, i) = nent(il, i) + 1
	214	END IF
	215
	216	Sij(il, i, j) = amax1(0.0, Sij(il,i,j))
	217	Sij(il, i, j) = amin1(1.0, Sij(il,i,j))
	218	END IF ! new
	219	END DO
	220	END DO
	221
	222	!jygdebug<
	223	IF (prt_level >= 10) THEN
	224	print *,'cv3p_mixing i, nent(i), icb, inb ',i, nent(igout,i), icb(igout), inb(igout)
	225	IF (nent(igout,i) .gt. 0) THEN
	226	print *,'i,(j,Sij(i,j),j=icb-1,inb) ',i,(j,Sij(igout,i,j),j=icb(igout)-1,inb(igout))
	227	ENDIF
	228	ENDIF
	229	!>jygdebug
	230
	231	! * if no air can entrain at level i assume that updraft detrains *
	232	! * at that level and calculate detrained air flux and properties *
	233
	234
	235	! @ do 170 i=icb(il),inb(il)
	236
	237	DO il = 1, ncum
	238	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
	239	! @ if(nent(il,i).eq.0)then
	240	!!! Ment(il,i,i)=m(il,i)
	241	Ment(il, i, i) = 1.
	242	Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
	243	uent(il, i, i) = unk(il)
	244	vent(il, i, i) = vnk(il)
	245	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
	246	Sij(il, i, i) = 0.0
	247	END IF
	248	END DO
	249	END DO
	250
	251	!jyg! DO j = 1, ntra
	252	!jyg! DO i = minorig + 1, nl
	253	!jyg! DO il = 1, ncum
	254	!jyg! IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
	255	!jyg! traent(il, i, i, j) = tra(il, nk(il), j)
	256	!jyg! END IF
	257	!jyg! END DO
	258	!jyg! END DO
	259	!jyg! END DO
	260
	261	DO j = minorig, nl
	262	DO i = minorig, nl
	263	DO il = 1, ncum
	264	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
	265	(i>=icb(il)) .AND. (i<=inb(il))) THEN
	266	Sigij(il, i, j) = Sij(il, i, j)
	267	END IF
	268	END DO
	269	END DO
	270	END DO
	271	! @ enddo
	272
	273	! @170 continue
	274
	275	! =====================================================================
	276	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
	277	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	278	! =====================================================================
	279
	280	CALL zilch(csum, nloc*nd)
	281
	282	DO il = 1, ncum
	283	lwork(il) = .FALSE.
	284	END DO
	285
	286	! ---------------------------------------------------------------
	287	DO i = minorig + 1, nl !Loop on origin level "i"
	288	! ---------------------------------------------------------------
	289
	290	num1 = 0
	291	DO il = 1, ncum
	292	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
	293	END DO
	294	IF (num1<=0) GO TO 789
	295
	296
	297	!JYG1 Find maximum of SIJ for J>I, if any.
	298
	299	Sx(:) = 0.
	300
	301	DO il = 1, ncum
	302	IF (i>=icb(il) .AND. i<=inb(il)) THEN
	303	signhpmh(il) = sign(1., hp(il,i)-h(il,i))
	304	Sbef(il) = max(0., signhpmh(il))
	305	END IF
	306	END DO
	307
	308	DO j = i + 1, nl
	309	DO il = 1, ncum
	310	IF (i>=icb(il) .AND. i<=inb(il) .AND. j<=inb(il)) THEN
	311	IF (Sbef(il)<Sij(il,i,j)) THEN
	312	Sx(il) = max(Sij(il,i,j), Sx(il))
	313	END IF
	314	Sbef(il) = Sij(il, i, j)
	315	END IF
	316	END DO
	317	END DO
	318
	319
	320	DO il = 1, ncum
	321	IF (i>=icb(il) .AND. i<=inb(il)) THEN
	322	lwork(il) = (nent(il,i)/=0)
	323	rti = qnk(il) - ep(il, i)*clw(il, i)
	324	!jyg<
	325	IF (cvflag_ice) THEN
	326
	327	anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)lf(il,i)) &
	328	(rti-rs(il,i)) + (cpv-cpd)t(il, i)(rti-rr(il,i))
	329	denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)lf(il,i)) &
	330	(rr(il,i)-rti) + (cpd-cpv)t(il, i)(rr(il,i)-rti)
	331	ELSE
	332
	333	anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + &
	334	(cpv-cpd)t(il, i)(rti-rr(il,i))
	335	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + &
	336	(cpd-cpv)t(il, i)(rr(il,i)-rti)
	337	END IF
	338	!>jyg
	339	IF (abs(denom)<0.01) denom = 0.01
	340	Scrit(il) = min(anum/denom, 1.)
	341	alt = rti - rs(il, i) + Scrit(il)*(rr(il,i)-rti)
	342
	343	!JYG1 Find new critical value Scrit2
	344	! such that : Sij > Scrit2 => mixed draught will detrain at J<I
	345	! Sij < Scrit2 => mixed draught will detrain at J>I
	346
	347	Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + &
	348	Scrit(il)*max(0., signhpmh(il))
	349
	350	Scrit(il) = Scrit2
	351
	352	!JYG Correction pour la nouvelle logique; la correction pour ALT
	353	! est un peu au hazard
	354	IF (Scrit(il)<=0.0) Scrit(il) = 0.0
	355	IF (alt<=0.0) Scrit(il) = 1.0
	356
	357	smax(il) = 0.0
	358	ASij(il) = 0.0
	359	sup(il) = 0. ! upper S-value reached by descending draughts
	360	END IF
	361	END DO
	362
	363	! ---------------------------------------------------------------
	364	DO j = minorig, nl !Loop on destination level "j"
	365	! ---------------------------------------------------------------
	366
	367	num2 = 0
	368	DO il = 1, ncum
	369	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
	370	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
	371	lwork(il)) num2 = num2 + 1
	372	END DO
	373	IF (num2<=0) GO TO 175
	374
	375	! -----------------------------------------------
	376	IF (j>i) THEN
	377	! -----------------------------------------------
	378	DO il = 1, ncum
	379	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
	380	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
	381	lwork(il)) THEN
	382	IF (Sij(il,i,j)>0.0) THEN
	383	Smid(il) = min(Sij(il,i,j), Scrit(il))
	384	Sjmax(il) = Smid(il)
	385	Sjmin(il) = Smid(il)
	386	IF (Smid(il)<smin(il) .AND. Sij(il,i,j+1)<Smid(il)) THEN
	387	smin(il) = Smid(il)
	388	Sjmax(il) = min((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j), Scrit(il))
	389	Sjmin(il) = max((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
	390	Sjmin(il) = min(Sjmin(il), Scrit(il))
	391	Sbef(il) = Sij(il, i, j)
	392	END IF
	393	END IF
	394	END IF
	395	END DO
	396	! -----------------------------------------------
	397	ELSE IF (j==i) THEN
	398	! -----------------------------------------------
	399	DO il = 1, ncum
	400	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
	401	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
	402	lwork(il)) THEN
	403	IF (Sij(il,i,j)>0.0) THEN
	404	Smid(il) = 1.
	405	Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + &
	406	min((Sij(il,i,j+1)+Smid(il))/2., Scrit(il))*max(0., signhpmh(il))
	407	Sjmin(il) = max(Sjmin(il), sup(il))
	408	Sjmax(il) = 1.
	409
	410	! - preparation des variables Scrit, Smin et Sbef pour la partie j>i
	411	Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
	412
	413	smin(il) = 1.
	414	Sbef(il) = max(0., signhpmh(il))
	415	supmax(il, i) = sign(Scrit(il), -signhpmh(il))
	416	END IF
	417	END IF
	418	END DO
	419	! -----------------------------------------------
	420	ELSE IF (j<i) THEN
	421	! -----------------------------------------------
	422	DO il = 1, ncum
	423	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
	424	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
	425	lwork(il)) THEN
	426	IF (Sij(il,i,j)>0.0) THEN
	427	Smid(il) = max(Sij(il,i,j), Scrit(il))
	428	Sjmax(il) = Smid(il)
	429	Sjmin(il) = Smid(il)
	430	IF (Smid(il)>smax(il) .AND. Sij(il,i,j+1)>Smid(il)) THEN
	431	smax(il) = Smid(il)
	432	Sjmax(il) = max((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j))
	433	Sjmax(il) = max(Sjmax(il), Scrit(il))
	434	Sjmin(il) = min((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
	435	Sjmin(il) = max(Sjmin(il), Scrit(il))
	436	Sbef(il) = Sij(il, i, j)
	437	END IF
	438	IF (abs(Sjmin(il)-Sjmax(il))>1.E-10) &
	439	sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
	440	END IF
	441	END IF
	442	END DO
	443	! -----------------------------------------------
	444	END IF
	445	! -----------------------------------------------
	446
	447
	448	DO il = 1, ncum
	449	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
	450	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
	451	lwork(il)) THEN
	452	IF (Sij(il,i,j)>0.0) THEN
	453	rti = qnk(il) - ep(il, i)*clw(il, i)
	454	Qmixmax(il) = Qmix(Sjmax(il))
	455	Qmixmin(il) = Qmix(Sjmin(il))
	456	Rmixmax(il) = Rmix(Sjmax(il))
	457	Rmixmin(il) = Rmix(Sjmin(il))
	458	sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
	459	sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
	460
	461	Ment(il, i, j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il, i, j)
	462
	463	! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
	464	IF (abs(Qmixmax(il)-Qmixmin(il))>1.E-10) THEN
	465	Sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/(Qmixmax(il)-Qmixmin(il))
	466	ELSE
	467	Sigij(il, i, j) = 0.
	468	END IF
	469
	470	! -- Compute Qent, uent, vent according to the true mixing fraction
	471	Qent(il, i, j) = (1.-Sigij(il,i,j))rti + Sigij(il, i, j)rr(il, i)
	472	uent(il, i, j) = (1.-Sigij(il,i,j))unk(il) + Sigij(il, i, j)u(il, i)
	473	vent(il, i, j) = (1.-Sigij(il,i,j))vnk(il) + Sigij(il, i, j)v(il, i)
	474
	475	! -- Compute liquid water static energy of mixed draughts
	476	! IF (j .GT. i) THEN
	477	! awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
	478	! awat=amax1(awat,0.0)
	479	! ELSE
	480	! awat = 0.
	481	! ENDIF
	482	! Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
	483	! : + Sigij(il,i,j)*H(il,i)
	484	! : + (LV(il,j)+(cpd-cpv)t(il,j))awat
	485	!IM 301008 beg
	486	hent(il, i, j) = (1.-Sigij(il,i,j))hp(il, i) + Sigij(il, i, j)h(il, i)
	487
	488	!jyg<
	489	! elij(il, i, j) = Qent(il, i, j) - rs(il, j)
	490	! elij(il, i, j) = elij(il, i, j) + &
	491	! ((h(il,j)-hent(il,i,j))rs(il,j)lv(il,j) / &
	492	! ((cpd(1.-Qent(il,i,j))+Qent(il,i,j)cpv)rrvt(il,j)*t(il,j)))
	493	! elij(il, i, j) = elij(il, i, j) / &
	494	! (1.+lv(il,j)lv(il,j)rs(il,j) / &
	495	! ((cpd(1.-Qent(il,i,j))+Qent(il,i,j)cpv)rrvt(il,j)*t(il,j)))
	496	!
	497	! Computation of condensate amount Elij, taking into account the ice fraction frac
	498	! Warning : the same saturation humidity rs is used over both liquid water and ice; this
	499	! should be corrected.
	500	!
	501	! Heat capacity of mixed draught
	502	cpm = cpd+Qent(il,i,j)*(cpv-cpd)
	503	!
	504	IF (cvflag_ice .and. frac(il,j) .gt. 0.) THEN
	505	elij(il, i, j) = Qent(il, i, j) - rs(il, j)
	506	elij(il, i, j) = elij(il, i, j) + &
	507	(h(il,j)-hent(il,i,j)+(cpv-cpd)(Qent(il,i,j)-rr(il,j))t(il,j))* &
	508	rs(il,j)lv(il,j) / (cpmrrvt(il,j)t(il,j))
	509	elij(il, i, j) = elij(il, i, j) / &
	510	(1.+(lv(il,j)+frac(il,j)lf(il,j))lv(il,j)*rs(il,j) / &
	511	(cpmrrvt(il,j)*t(il,j)))
	512	ELSE
	513	elij(il, i, j) = Qent(il, i, j) - rs(il, j)
	514	elij(il, i, j) = elij(il, i, j) + &
	515	(h(il,j)-hent(il,i,j)+(cpv-cpd)(Qent(il,i,j)-rr(il,j))t(il,j))* &
	516	rs(il,j)lv(il,j) / (cpmrrvt(il,j)t(il,j))
	517	elij(il, i, j) = elij(il, i, j) / &
	518	(1.+lv(il,j)lv(il,j)rs(il,j) / &
	519	(cpmrrvt(il,j)*t(il,j)))
	520	ENDIF
	521	!>jyg
	522	elij(il, i, j) = max(elij(il,i,j), 0.)
	523
	524	elij(il, i, j) = min(elij(il,i,j), Qent(il,i,j))
	525
	526	IF (j>i) THEN
	527	awat = elij(il, i, j) - (1.-ep(il,j))*clw(il, j)
	528	awat = amax1(awat, 0.0)
	529	ELSE
	530	awat = 0.
	531	END IF
	532
	533	! print ,h(il,j)-hent(il,i,j),LV(il,j)rs(il,j)/(cpdrrvt(il,j)*
	534	! : t(il,j))
	535
	536	!jyg<
	537	! hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)t(il,j))awat
	538	! Mixed draught temperature at level j
	539	IF (cvflag_ice .and. frac(il,j) .gt. 0.) THEN
	540	Tm = t(il,j) + (Qent(il,i,j)-elij(il,i,j)-rs(il,j))rrvt(il,j)t(il,j)/(lv(il,j)rs(il,j))
	541	hent(il, i, j) = hent(il, i, j) + (lv(il,j)+frac(il,j)lf(il,j)+(cpd-cpv)Tm)*awat
	542	ELSE
	543	Tm = t(il,j) + (Qent(il,i,j)-elij(il,i,j)-rs(il,j))rrvt(il,j)t(il,j)/(lv(il,j)rs(il,j))
	544	hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)Tm)awat
	545	ENDIF
	546	!>jyg
	547
	548	!IM 301008 end
	549
	550	! print *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ',
	551	! : i,j,hent(il,i,j),Sigij(il,i,j)
	552
	553	! -- ASij is the integral of P(F) over the relevant F interval
	554	ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
	555	Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
	556
	557	END IF
	558	END IF
	559	END DO
	560	!jyg! DO k = 1, ntra
	561	!jyg! DO il = 1, ncum
	562	!jyg! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
	563	!jyg! (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
	564	!jyg! lwork(il)) THEN
	565	!jyg! IF (Sij(il,i,j)>0.0) THEN
	566	!jyg! traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + &
	567	!jyg! (1.-Sigij(il,i,j))*tra(il, nk(il), k)
	568	!jyg! END IF
	569	!jyg! END IF
	570	!jyg! END DO
	571	!jyg! END DO
	572
	573	! -- If I=J (detrainement and entrainement at the same level), then only the
	574	! -- adiabatic ascent part of the mixture is considered
	575	IF (i==j) THEN
	576	DO il = 1, ncum
	577	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
	578	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
	579	lwork(il)) THEN
	580	IF (Sij(il,i,j)>0.0) THEN
	581	rti = qnk(il) - ep(il, i)*clw(il, i)
	582	!!! Ment(il,i,i) = m(il,i)abs(Qmixmax(il)(1.-Sjmax(il))
	583	Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
	584	Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
	585	Qent(il, i, i) = rti
	586	uent(il, i, i) = unk(il)
	587	vent(il, i, i) = vnk(il)
	588	hent(il, i, i) = hp(il, i)
	589	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
	590	Sigij(il, i, i) = 0.
	591	END IF
	592	END IF
	593	END DO
	594	!jyg! DO k = 1, ntra
	595	!jyg! DO il = 1, ncum
	596	!jyg! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
	597	!jyg! (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
	598	!jyg! lwork(il)) THEN
	599	!jyg! IF (Sij(il,i,j)>0.0) THEN
	600	!jyg! traent(il, i, i, k) = tra(il, nk(il), k)
	601	!jyg! END IF
	602	!jyg! END IF
	603	!jyg! END DO
	604	!jyg! END DO
	605
	606	END IF
	607
	608	! ---------------------------------------------------------------
	609	175 END DO ! End loop on destination level "j"
	610	! ---------------------------------------------------------------
	611
	612	DO il = 1, ncum
	613	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
	614	ASij(il) = amax1(1.0E-16, ASij(il))
	615	!jyg+lluis<
	616	!! ASij(il) = 1.0/ASij(il)
	617	ASij_inv(il) = 1.0/ASij(il)
	618	! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
	619	IF (ASij_inv(il) > 100.) ASij_inv(il) = 0.
	620	!>jyg+lluis
	621	csum(il, i) = 0.0
	622	END IF
	623	END DO
	624
	625	DO j = minorig, nl
	626	DO il = 1, ncum
	627	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	628	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
	629	!jyg Ment(il, i, j) = Ment(il, i, j)*ASij(il)
	630	Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il)
	631	END IF
	632	END DO
	633	END DO
	634
	635	DO j = minorig, nl
	636	DO il = 1, ncum
	637	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	638	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
	639	csum(il, i) = csum(il, i) + Ment(il, i, j)
	640	END IF
	641	END DO
	642	END DO
	643
	644	DO il = 1, ncum
	645	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
	646	! cc : .and. csum(il,i).lt.m(il,i) ) then
	647	nent(il, i) = 0
	648	! cc Ment(il,i,i)=m(il,i)
	649	Ment(il, i, i) = 1.
	650	Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
	651	uent(il, i, i) = unk(il)
	652	vent(il, i, i) = vnk(il)
	653	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
	654	Sij(il, i, i) = 0.0
	655	END IF
	656	END DO ! il
	657
	658	!jyg! DO j = 1, ntra
	659	!jyg! DO il = 1, ncum
	660	!jyg! IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
	661	!jyg!! cc : .and. csum(il,i).lt.m(il,i) ) then
	662	!jyg! traent(il, i, i, j) = tra(il, nk(il), j)
	663	!jyg! END IF
	664	!jyg! END DO
	665	!jyg! END DO
	666
	667	! ---------------------------------------------------------------
	668	789 END DO ! End loop on origin level "i"
	669	! ---------------------------------------------------------------
	670
	671
	672	RETURN
	673	END SUBROUTINE cv3p_mixing
	674

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

Download in other formats:

Original Format