Context Navigation

cv3p_mixing.F90 @ 2477

Last change on this file since 2477 was 2397, checked in by jyg, 9 years ago
Implementation of ice thermodynamics in cv3p_mixing.F90 . Needs debugging.
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: 24.2 KB

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

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

Download in other formats:

Original Format