Context Navigation

cv3p_mixing.F90 @ 5445

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

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

Download in other formats:

Original Format