Context Navigation

cv3p_mixing.f90 @ 5625

Last change on this file since 5625 was 5346, checked in by fhourdin, 12 months ago

Debut de replaysation de la convection profonde.

Regroupement de cvparam, cv3param et cvthermo (récemment
passés de statut de .h à module, dans un unique module
lmdz_cv_ini.f90

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

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

Download in other formats:

Original Format