Context Navigation

cv3p_mixing.F90 @ 2591

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

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

Download in other formats:

Original Format