Context Navigation

cv3p_mixing.F90 @ 5051

Last change on this file since 5051 was 3496, checked in by jyg, 6 years ago

Implementation of the ejection of liquid precipitation from the adiabatic ascents.
New flags:
+cvflag_prec_eject: logical

n -> old code, y -> new code

+ejectliq: real; possible values 0. & 1.

-> no liquid precipitation is ejected
-> all liquid precipitation is ejected

+ejectice: real; any value between 0. and 1.

fraction of solid precipitation ejected at each level

Note that the adiabatic ascent mass flux decrease due to precipitation ejection is not taken into account.

Attempts to do it led to water conservation violation.

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: 25.3 KB

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

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

Download in other formats:

Original Format