Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

cv3p_mixing.f90

Last change on this file was 5703, checked in by yann meurdesoif, 2 months ago

Convection GPU porting : change some loop name indices to avoid GPUmorphosis confusion (solve in future)

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.6 KB

Rev	Line
[5692]	1	MODULE cv3p_mixing_mod
	2	PRIVATE
	3	PUBLIC cv3p_mixing
	4	CONTAINS
	5
[5697]	6	SUBROUTINE cv3p_mixing_pre
[5692]	7
[5697]	8	END SUBROUTINE cv3p_mixing_pre
	9
[2007]	10	SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
[3496]	11	ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, &
[2007]	12	unk, vnk, hp, tv, tvp, ep, clw, sig, &
	13	Ment, Qent, hent, uent, vent, nent, &
	14	Sigij, elij, supmax, Ments, Qents, traent)
	15	! **************************************************************
	16	! *
	17	! CV3P_MIXING : compute mixed draught properties and, *
	18	! within a scaling factor, mixed draught *
	19	! mass fluxes. *
	20	! written by : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15*
	21	! modified by : *
	22	! **************************************************************
[879]	23
[5697]	24	USE yomcst2_mod_h, ONLY : Fmax, gammas, scut, qqa1, qqa2
[5346]	25	USE lmdz_cv_ini, ONLY : cpd,cpv,minorig,nl,rrv
[5285]	26	USE cvflag_mod_h
[2393]	27	USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
[2900]	28	USE ioipsl_getin_p_mod, ONLY: getin_p
	29	USE add_phys_tend_mod, ONLY: fl_cor_ebil
[2393]	30
[1992]	31	IMPLICIT NONE
[879]	32
	33
[2007]	34	!inputs:
	35	INTEGER, INTENT (IN) :: ncum, nd, na
	36	INTEGER, INTENT (IN) :: ntra, nloc
	37	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
	38	REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig
[3496]	39	REAL, DIMENSION (nloc), INTENT (IN) :: unk, vnk
	40	REAL, DIMENSION (nloc, nd), INTENT (IN) :: qta
[2007]	41	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
	42	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
	43	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
	44	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra ! input of convect3
	45	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv
[2397]	46	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
	47	REAL, DIMENSION (nloc, na), INTENT (IN) :: frac !ice fraction in condensate
[2393]	48	REAL, DIMENSION (nloc, na), INTENT (IN) :: h !liquid water static energy of environment
[2007]	49	REAL, DIMENSION (nloc, na), INTENT (IN) :: hp !liquid water static energy of air shed from adiab. asc.
	50	REAL, DIMENSION (nloc, na), INTENT (IN) :: tv, tvp
	51	REAL, DIMENSION (nloc, na), INTENT (IN) :: ep, clw
[879]	52
[2007]	53	!outputs:
	54	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: Ment, Qent
	55	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: uent, vent
	56	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: Sigij, elij
[2393]	57	REAL, DIMENSION (nloc, na), INTENT (OUT) :: supmax ! Highest mixing fraction of mixed
[2007]	58	! updraughts with the sign of (h-hp)
	59	REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
	60	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: Ments, Qents
	61	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: hent
	62	INTEGER, DIMENSION (nloc, nd), INTENT (OUT) :: nent
[879]	63
[2007]	64	!local variables:
[1992]	65	INTEGER i, j, k, il, im, jm
	66	INTEGER num1, num2
[2397]	67	REAL :: rti, bf2, anum, denom, dei, altem, cwat, stemp
[2007]	68	REAL :: alt, delp, delm
	69	REAL, DIMENSION (nloc) :: Qmixmax, Rmixmax, sqmrmax
	70	REAL, DIMENSION (nloc) :: Qmixmin, Rmixmin, sqmrmin
	71	REAL, DIMENSION (nloc) :: signhpmh
	72	REAL, DIMENSION (nloc) :: Sx
	73	REAL :: Scrit2
	74	REAL, DIMENSION (nloc) :: Smid, Sjmin, Sjmax
	75	REAL, DIMENSION (nloc) :: Sbef, sup, smin
[2226]	76	REAL, DIMENSION (nloc) :: ASij, ASij_inv, smax, Scrit
[2007]	77	REAL, DIMENSION (nloc, nd, nd) :: Sij
	78	REAL, DIMENSION (nloc, nd) :: csum
	79	REAL :: awat
[2397]	80	REAL :: cpm !Mixed draught heat capacity
	81	REAL :: Tm !Mixed draught temperature
[2007]	82	LOGICAL, DIMENSION (nloc) :: lwork
[879]	83
[1992]	84	REAL amxupcrit, df, ff
	85	INTEGER nstep
[879]	86
[5695]	87	INTEGER,PARAMETER :: igout=1
[2393]	88
[2007]	89	! -- Mixing probability distribution functions
[1650]	90
[2007]	91	REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
[5697]	92	REAL :: Qcoef1max,Qcoef2max !ym WARNING
	93	!ym redefine local variable instead use module variable
	94	!ym to check when refactoring deep convection
	95	!ym => eliminate "first" SAVE variable
	96	!ym probably all these folowing lines will be removed
[2007]	97	Qcoef1(F) = tanh(F/gammas)
	98	Qcoef2(F) = (tanh(F/gammas)+gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
	99	QFF(F) = max(min(F,1.), 0.)
	100	QFFf(F) = min(QFF(F), scut)
	101	Qmix1(F) = (tanh((QFF(F)-Fmax)/gammas)+Qcoef1max)/Qcoef2max
	102	Rmix1(F) = (gammaslog(cosh((QFF(F)-Fmax)/gammas))+QFF(F)Qcoef1max)/Qcoef2max
	103	Qmix2(F) = -log(1.-QFFf(F))/scut
	104	Rmix2(F) = (QFFf(F)+(1.-QFF(F))*log(1.-QFFf(F)))/scut
	105	Qmix(F) = qqa1Qmix1(F) + qqa2Qmix2(F)
	106	Rmix(F) = qqa1Rmix1(F) + qqa2Rmix2(F)
[879]	107
	108
	109
[2007]	110	! =====================================================================
	111	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	112	! =====================================================================
[879]	113
[5697]	114	Qcoef1max = Qcoef1(Fmax)
	115	Qcoef2max = Qcoef2(Fmax)
[879]	116
[2007]	117	! ori do 360 i=1,ncum*nlp
[1992]	118	DO j = 1, nl
[5703]	119	DO il = 1, ncum
	120	nent(il, j) = 0
[2007]	121	! in convect3, m is computed in cv3_closure
	122	! ori m(i,1)=0.0
[1992]	123	END DO
	124	END DO
[879]	125
[2007]	126	! ori do 400 k=1,nlp
	127	! ori do 390 j=1,nlp
[1992]	128	DO j = 1, nl
	129	DO k = 1, nl
[5703]	130	DO il = 1, ncum
	131	Qent(il, k, j) = rr(il, j)
	132	uent(il, k, j) = u(il, j)
	133	vent(il, k, j) = v(il, j)
	134	elij(il, k, j) = 0.0
	135	hent(il, k, j) = 0.0
[2007]	136	!AC! Ment(i,k,j)=0.0
	137	!AC! Sij(i,k,j)=0.0
[1992]	138	END DO
	139	END DO
	140	END DO
[879]	141
[2007]	142	!AC!
	143	Ment(1:ncum, 1:nd, 1:nd) = 0.0
	144	Sij(1:ncum, 1:nd, 1:nd) = 0.0
	145	!AC!
[2478]	146	!ym
	147	Sigij(1:ncum, 1:nd, 1:nd) = 0.0
	148	!ym
[879]	149
[2393]	150	!jyg! DO k = 1, ntra
	151	!jyg! DO j = 1, nd ! instead nlp
	152	!jyg! DO i = 1, nd ! instead nlp
	153	!jyg! DO il = 1, ncum
	154	!jyg! traent(il, i, j, k) = tra(il, j, k)
	155	!jyg! END DO
	156	!jyg! END DO
	157	!jyg! END DO
	158	!jyg! END DO
[879]	159
[2007]	160	! =====================================================================
	161	! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
	162	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	163	! --- FRACTION (Sij)
	164	! =====================================================================
[879]	165
[1992]	166	DO i = minorig + 1, nl
[879]	167
[2900]	168	IF (ok_entrain) THEN
	169	DO j = minorig, nl
	170	DO il = 1, ncum
	171	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
	172	.AND. (j<=inb(il))) THEN
[879]	173
[3496]	174	!! rti = qnk(il) - ep(il, i)*clw(il, i)
	175	rti = qta(il,i-1) - ep(il, i)*clw(il, i)
[2900]	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))
[3496]	221	ELSE IF (j > i) THEN
	222	IF (prt_level >= 10) THEN
	223	print *,'cv3p_mixing i, j, Sij given by the no-precip eq. ', i, j, Sij(il,i,j)
	224	ENDIF
[2900]	225	END IF ! new
	226	END DO
[1992]	227	END DO
[2900]	228	ELSE ! (ok_entrain)
	229	DO il = 1,ncum
	230	nent(il,i) = 0
	231	ENDDO
	232	ENDIF ! (ok_entrain)
[879]	233
[2393]	234	!jygdebug<
	235	IF (prt_level >= 10) THEN
	236	print *,'cv3p_mixing i, nent(i), icb, inb ',i, nent(igout,i), icb(igout), inb(igout)
	237	IF (nent(igout,i) .gt. 0) THEN
	238	print *,'i,(j,Sij(i,j),j=icb-1,inb) ',i,(j,Sij(igout,i,j),j=icb(igout)-1,inb(igout))
	239	ENDIF
	240	ENDIF
	241	!>jygdebug
[879]	242
[2007]	243	! * if no air can entrain at level i assume that updraft detrains *
	244	! * at that level and calculate detrained air flux and properties *
[1494]	245
[879]	246
[2007]	247	! @ do 170 i=icb(il),inb(il)
[879]	248
[1992]	249	DO il = 1, ncum
	250	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
[2007]	251	! @ if(nent(il,i).eq.0)then
	252	!!! Ment(il,i,i)=m(il,i)
	253	Ment(il, i, i) = 1.
[3496]	254	!! Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
	255	Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
[1992]	256	uent(il, i, i) = unk(il)
	257	vent(il, i, i) = vnk(il)
[2902]	258	IF (fl_cor_ebil .GE. 2) THEN
	259	hent(il, i, i) = hp(il,i)
	260	ENDIF
[1992]	261	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
[2007]	262	Sij(il, i, i) = 0.0
[1992]	263	END IF
	264	END DO
[2902]	265	END DO ! i = minorig + 1, nl
[879]	266
[2393]	267	!jyg! DO j = 1, ntra
	268	!jyg! DO i = minorig + 1, nl
	269	!jyg! DO il = 1, ncum
	270	!jyg! IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
	271	!jyg! traent(il, i, i, j) = tra(il, nk(il), j)
	272	!jyg! END IF
	273	!jyg! END DO
	274	!jyg! END DO
	275	!jyg! END DO
[879]	276
[1992]	277	DO j = minorig, nl
	278	DO i = minorig, nl
	279	DO il = 1, ncum
[2007]	280	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
	281	(i>=icb(il)) .AND. (i<=inb(il))) THEN
	282	Sigij(il, i, j) = Sij(il, i, j)
[1992]	283	END IF
	284	END DO
	285	END DO
	286	END DO
[2007]	287	! @ enddo
[1041]	288
[2007]	289	! @170 continue
[1041]	290
[2007]	291	! =====================================================================
	292	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
	293	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	294	! =====================================================================
[1041]	295
[5699]	296	csum(:,:) = 0.
	297
[1992]	298	DO il = 1, ncum
	299	lwork(il) = .FALSE.
	300	END DO
[1041]	301
[2007]	302	! ---------------------------------------------------------------
	303	DO i = minorig + 1, nl !Loop on origin level "i"
	304	! ---------------------------------------------------------------
[1041]	305
[1992]	306	num1 = 0
	307	DO il = 1, ncum
	308	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
	309	END DO
[5701]	310	!ym IF (num1<=0) GO TO 789
	311	IF (num1<=0) CYCLE
[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
[5701]	391	!ym IF (num2<=0) GO TO 175
	392	IF (num2<=0) CYCLE
[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	! ---------------------------------------------------------------
[5701]	629	END DO !ym label 175 ! End loop on destination level "j"
[2007]	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	! ---------------------------------------------------------------
[5701]	694	END DO !ym label 789 ! End loop on origin level "i"
	695
[2007]	696	! ---------------------------------------------------------------
[1992]	697
[2007]	698
[1992]	699	RETURN
	700	END SUBROUTINE cv3p_mixing
	701
[5692]	702	END MODULE cv3p_mixing_mod

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

Download in other formats:

Original Format