Context Navigation

cv3p_mixing.F90 @ 3605

Last change on this file since 3605 was 3605, checked in by lguez, 5 years ago
Merge revisions 3427:3600 of trunk into branch `Ocean_skin`
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, &
[3605]	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
[3605]	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
[3605]	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))
[3605]	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.
[3605]	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)
[3605]	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
[3605]	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
[3605]	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.
[3605]	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