Context Navigation

cv3p_mixing.F90 @ 1992

Last change on this file since 1992 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

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: 19.0 KB

Rev	Line
[1992]	1	SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, ph, t, rr, rs, &
	2	u, v, tra, h, lv, qnk, unk, vnk, hp, tv, tvp, ep, clw, sig, ment, qent, &
	3	hent, uent, vent, nent, sigij, elij, supmax, ments, qents, traent)
	4	! **************************************************************
	5	! *
	6	! CV3P_MIXING : compute mixed draught properties and, *
	7	! within a scaling factor, mixed draught *
	8	! mass fluxes. *
	9	! written by : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15*
	10	! modified by : *
	11	! **************************************************************
[879]	12
[1992]	13	IMPLICIT NONE
[879]	14
[1992]	15	include "cvthermo.h"
	16	include "cv3param.h"
	17	include "YOMCST2.h"
[879]	18
[1992]	19	! inputs:
	20	INTEGER ncum, nd, na, ntra, nloc
	21	INTEGER icb(nloc), inb(nloc), nk(nloc)
	22	REAL sig(nloc, nd)
	23	REAL qnk(nloc), unk(nloc), vnk(nloc)
	24	REAL ph(nloc, nd+1)
	25	REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
	26	REAL u(nloc, nd), v(nloc, nd)
	27	REAL tra(nloc, nd, ntra) ! input of convect3
	28	REAL lv(nloc, na)
	29	REAL h(nloc, na) !liquid water static energy of environment
	30	REAL hp(nloc, na) !liquid water static energy of air shed from adiab. asc.
	31	REAL tv(nloc, na), tvp(nloc, na), ep(nloc, na), clw(nloc, na)
[879]	32
[1992]	33	! outputs:
	34	REAL ment(nloc, na, na), qent(nloc, na, na)
	35	REAL uent(nloc, na, na), vent(nloc, na, na)
	36	REAL sigij(nloc, na, na), elij(nloc, na, na)
	37	REAL supmax(nloc, na) ! Highest mixing fraction of mixed updraughts
	38	! with the sign of (h-hp)
	39	REAL traent(nloc, nd, nd, ntra)
	40	REAL ments(nloc, nd, nd), qents(nloc, nd, nd)
	41	REAL hent(nloc, nd, nd)
	42	INTEGER nent(nloc, nd)
[879]	43
[1992]	44	! local variables:
	45	INTEGER i, j, k, il, im, jm
	46	INTEGER num1, num2
	47	REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
	48	REAL alt, delp, delm
	49	REAL qmixmax(nloc), rmixmax(nloc), sqmrmax(nloc)
	50	REAL qmixmin(nloc), rmixmin(nloc), sqmrmin(nloc)
	51	REAL signhpmh(nloc)
	52	REAL sx(nloc), scrit2
	53	REAL smid(nloc), sjmin(nloc), sjmax(nloc)
	54	REAL sbef(nloc), sup(nloc), smin(nloc)
	55	REAL asij(nloc), smax(nloc), scrit(nloc)
	56	REAL sij(nloc, nd, nd)
	57	REAL csum(nloc, nd)
	58	REAL awat
	59	LOGICAL lwork(nloc)
[879]	60
[1992]	61	REAL amxupcrit, df, ff
	62	INTEGER nstep
[879]	63
[1992]	64	! -- Mixing probability distribution functions
[1650]	65
[1992]	66	REAL qcoef1, qcoef2, qff, qfff, qmix, rmix, qmix1, rmix1, qmix2, rmix2, f
[879]	67
[1992]	68	qcoef1(f) = tanh(f/gammas)
	69	qcoef2(f) = (tanh(f/gammas)+gammas*log(cosh((1.-f)/gammas)/cosh(f/gammas)))
	70	qff(f) = max(min(f,1.), 0.)
	71	qfff(f) = min(qff(f), scut)
	72	qmix1(f) = (tanh((qff(f)-fmax)/gammas)+qcoef1max)/qcoef2max
	73	rmix1(f) = (gammaslog(cosh((qff(f)-fmax)/gammas))+qff(f)qcoef1max)/ &
	74	qcoef2max
	75	qmix2(f) = -log(1.-qfff(f))/scut
	76	rmix2(f) = (qfff(f)+(1.-qff(f))*log(1.-qfff(f)))/scut
	77	qmix(f) = qqa1qmix1(f) + qqa2qmix2(f)
	78	rmix(f) = qqa1rmix1(f) + qqa2rmix2(f)
[879]	79
[1992]	80	INTEGER, SAVE :: ifrst
	81	DATA ifrst/0/
	82	!$OMP THREADPRIVATE(ifrst)
[879]	83
	84
[1992]	85	! =====================================================================
	86	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	87	! =====================================================================
[879]	88
[1992]	89	! -- Initialize mixing PDF coefficients
	90	IF (ifrst==0) THEN
	91	ifrst = 1
	92	qcoef1max = qcoef1(fmax)
	93	qcoef2max = qcoef2(fmax)
[879]	94
[1992]	95	END IF
[879]	96
	97
[1992]	98	! ori do 360 i=1,ncum*nlp
	99	DO j = 1, nl
	100	DO i = 1, ncum
	101	nent(i, j) = 0
	102	! in convect3, m is computed in cv3_closure
	103	! ori m(i,1)=0.0
	104	END DO
	105	END DO
[879]	106
[1992]	107	! ori do 400 k=1,nlp
	108	! ori do 390 j=1,nlp
	109	DO j = 1, nl
	110	DO k = 1, nl
	111	DO i = 1, ncum
	112	qent(i, k, j) = rr(i, j)
	113	uent(i, k, j) = u(i, j)
	114	vent(i, k, j) = v(i, j)
	115	elij(i, k, j) = 0.0
	116	hent(i, k, j) = 0.0
	117	! AC! ment(i,k,j)=0.0
	118	! AC! sij(i,k,j)=0.0
	119	END DO
	120	END DO
	121	END DO
[879]	122
[1992]	123	! AC!
	124	ment(1:ncum, 1:nd, 1:nd) = 0.0
	125	sij(1:ncum, 1:nd, 1:nd) = 0.0
	126	! AC!
[879]	127
[1992]	128	DO k = 1, ntra
	129	DO j = 1, nd ! instead nlp
	130	DO i = 1, nd ! instead nlp
	131	DO il = 1, ncum
	132	traent(il, i, j, k) = tra(il, j, k)
	133	END DO
	134	END DO
	135	END DO
	136	END DO
[879]	137
[1992]	138	! =====================================================================
	139	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
	140	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	141	! --- FRACTION (sij)
	142	! =====================================================================
[879]	143
[1992]	144	DO i = minorig + 1, nl
[879]	145
[1992]	146	DO j = minorig, nl
	147	DO il = 1, ncum
	148	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
	149	1)) .AND. (j<=inb(il))) THEN
[879]	150
[1992]	151	rti = qnk(il) - ep(il, i)*clw(il, i)
	152	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
	153	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
	154	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
	155	dei = denom
	156	IF (abs(dei)<0.01) dei = 0.01
	157	sij(il, i, j) = anum/dei
	158	sij(il, i, i) = 1.0
	159	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
	160	altem = altem/bf2
	161	cwat = clw(il, j)*(1.-ep(il,j))
	162	stemp = sij(il, i, j)
	163	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
	164	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
	165	denom = denom + lv(il, j)*(rr(il,i)-rti)
	166	IF (abs(denom)<0.01) denom = 0.01
	167	sij(il, i, j) = anum/denom
	168	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - &
	169	rs(il, j)
	170	altem = altem - (bf2-1.)*cwat
	171	END IF
	172	IF (sij(il,i,j)>0.0) THEN
	173	! cc ment(il,i,j)=m(il,i)
	174	ment(il, i, j) = 1.
	175	elij(il, i, j) = altem
	176	elij(il, i, j) = amax1(0.0, elij(il,i,j))
	177	nent(il, i) = nent(il, i) + 1
	178	END IF
[879]	179
[1992]	180	sij(il, i, j) = amax1(0.0, sij(il,i,j))
	181	sij(il, i, j) = amin1(1.0, sij(il,i,j))
	182	END IF ! new
	183	END DO
	184	END DO
[879]	185
	186
[1992]	187	! *** if no air can entrain at level i assume that updraft detrains
	188	! ***
	189	! *** at that level and calculate detrained air flux and properties
	190	! ***
[1494]	191
[879]	192
[1992]	193	! @ do 170 i=icb(il),inb(il)
[879]	194
[1992]	195	DO il = 1, ncum
	196	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
	197	! @ if(nent(il,i).eq.0)then
	198	! cc ment(il,i,i)=m(il,i)
	199	ment(il, i, i) = 1.
	200	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
	201	uent(il, i, i) = unk(il)
	202	vent(il, i, i) = vnk(il)
	203	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
	204	sij(il, i, i) = 0.0
	205	END IF
	206	END DO
	207	END DO
[879]	208
[1992]	209	DO j = 1, ntra
	210	DO i = minorig + 1, nl
	211	DO il = 1, ncum
	212	IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
	213	traent(il, i, i, j) = tra(il, nk(il), j)
	214	END IF
	215	END DO
	216	END DO
	217	END DO
[879]	218
[1992]	219	DO j = minorig, nl
	220	DO i = minorig, nl
	221	DO il = 1, ncum
	222	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<= &
	223	inb(il))) THEN
	224	sigij(il, i, j) = sij(il, i, j)
	225	END IF
	226	END DO
	227	END DO
	228	END DO
	229	! @ enddo
[1041]	230
[1992]	231	! @170 continue
[1041]	232
[1992]	233	! =====================================================================
	234	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
	235	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	236	! =====================================================================
[1041]	237
[1992]	238	CALL zilch(csum, nloc*nd)
[1041]	239
[1992]	240	DO il = 1, ncum
	241	lwork(il) = .FALSE.
	242	END DO
[1041]	243
[1992]	244	! ---------------------------------------------------------------
	245	DO i = minorig + 1, nl !Loop on origin level "i"
	246	! ---------------------------------------------------------------
[1041]	247
[1992]	248	num1 = 0
	249	DO il = 1, ncum
	250	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
	251	END DO
	252	IF (num1<=0) GO TO 789
[879]	253
	254
[1992]	255	! jyg1 Find maximum of SIJ for J>I, if any.
[879]	256
[1992]	257	sx(:) = 0.
[879]	258
[1992]	259	DO il = 1, ncum
	260	IF (i>=icb(il) .AND. i<=inb(il)) THEN
	261	signhpmh(il) = sign(1., hp(il,i)-h(il,i))
	262	sbef(il) = max(0., signhpmh(il))
	263	END IF
	264	END DO
[879]	265
[1992]	266	DO j = i + 1, nl
	267	DO il = 1, ncum
	268	IF (i>=icb(il) .AND. i<=inb(il) .AND. j<=inb(il)) THEN
	269	IF (sbef(il)<sij(il,i,j)) THEN
	270	sx(il) = max(sij(il,i,j), sx(il))
	271	END IF
	272	sbef(il) = sij(il, i, j)
	273	END IF
	274	END DO
	275	END DO
[879]	276
[1992]	277
	278	DO il = 1, ncum
	279	IF (i>=icb(il) .AND. i<=inb(il)) THEN
	280	lwork(il) = (nent(il,i)/=0)
	281	qp = qnk(il) - ep(il, i)*clw(il, i)
	282	anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
	283	(cpv-cpd)t(il, i)(qp-rr(il,i))
	284	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
	285	(cpd-cpv)t(il, i)(rr(il,i)-qp)
	286	IF (abs(denom)<0.01) denom = 0.01
	287	scrit(il) = min(anum/denom, 1.)
	288	alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
	289
	290	! jyg1 Find new critical value Scrit2
	291	! such that : Sij > Scrit2 => mixed draught will detrain at J<I
	292	! Sij < Scrit2 => mixed draught will detrain at J>I
	293
	294	scrit2 = min(scrit(il), sx(il))*max(0., -signhpmh(il)) + &
	295	scrit(il)*max(0., signhpmh(il))
	296
	297	scrit(il) = scrit2
	298
	299	! jyg Correction pour la nouvelle logique; la correction pour ALT
	300	! est un peu au hazard
	301	IF (scrit(il)<=0.0) scrit(il) = 0.0
	302	IF (alt<=0.0) scrit(il) = 1.0
	303
	304	smax(il) = 0.0
	305	asij(il) = 0.0
	306	sup(il) = 0. ! upper S-value reached by descending draughts
	307	END IF
	308	END DO
	309
	310	! ---------------------------------------------------------------
	311	DO j = minorig, nl !Loop on destination level "j"
	312	! ---------------------------------------------------------------
	313
	314	num2 = 0
	315	DO il = 1, ncum
	316	IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
	317	il)-1) .AND. j<=inb(il) .AND. lwork(il)) num2 = num2 + 1
	318	END DO
	319	IF (num2<=0) GO TO 175
	320
	321	! -----------------------------------------------
	322	IF (j>i) THEN
	323	! -----------------------------------------------
	324	DO il = 1, ncum
	325	IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
	326	il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
	327	IF (sij(il,i,j)>0.0) THEN
	328	smid(il) = min(sij(il,i,j), scrit(il))
	329	sjmax(il) = smid(il)
	330	sjmin(il) = smid(il)
	331	IF (smid(il)<smin(il) .AND. sij(il,i,j+1)<smid(il)) THEN
	332	smin(il) = smid(il)
	333	sjmax(il) = min((sij(il,i,j+1)+sij(il,i, &
	334	j))/2., sij(il,i,j), scrit(il))
	335	sjmin(il) = max((sbef(il)+sij(il,i,j))/2., sij(il,i,j))
	336	sjmin(il) = min(sjmin(il), scrit(il))
	337	sbef(il) = sij(il, i, j)
	338	END IF
	339	END IF
	340	END IF
	341	END DO
	342	! -----------------------------------------------
	343	ELSE IF (j==i) THEN
	344	! -----------------------------------------------
	345	DO il = 1, ncum
	346	IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
	347	il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
	348	IF (sij(il,i,j)>0.0) THEN
	349	smid(il) = 1.
	350	sjmin(il) = max((sij(il,i,j-1)+smid(il))/2., scrit(il))*max(0., &
	351	-signhpmh(il)) + min((sij(il,i,j+1)+smid(il))/2., scrit(il))* &
	352	max(0., signhpmh(il))
	353	sjmin(il) = max(sjmin(il), sup(il))
	354	sjmax(il) = 1.
	355
	356	! - preparation des variables Scrit, Smin et Sbef
	357	! pour la partie j>i
	358	scrit(il) = min(sjmin(il), sjmax(il), scrit(il))
	359
	360	smin(il) = 1.
	361	sbef(il) = max(0., signhpmh(il))
	362	supmax(il, i) = sign(scrit(il), -signhpmh(il))
	363	END IF
	364	END IF
	365	END DO
	366	! -----------------------------------------------
	367	ELSE IF (j<i) THEN
	368	! -----------------------------------------------
	369	DO il = 1, ncum
	370	IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
	371	il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
	372	IF (sij(il,i,j)>0.0) THEN
	373	smid(il) = max(sij(il,i,j), scrit(il))
	374	sjmax(il) = smid(il)
	375	sjmin(il) = smid(il)
	376	IF (smid(il)>smax(il) .AND. sij(il,i,j+1)>smid(il)) THEN
	377	smax(il) = smid(il)
	378	sjmax(il) = max((sij(il,i,j+1)+sij(il,i,j))/2., sij(il,i,j))
	379	sjmax(il) = max(sjmax(il), scrit(il))
	380	sjmin(il) = min((sbef(il)+sij(il,i,j))/2., sij(il,i,j))
	381	sjmin(il) = max(sjmin(il), scrit(il))
	382	sbef(il) = sij(il, i, j)
	383	END IF
	384	IF (abs(sjmin(il)-sjmax(il))>1.E-10) sup(il) = max(sjmin(il), &
	385	sjmax(il), sup(il))
	386	END IF
	387	END IF
	388	END DO
	389	! -----------------------------------------------
	390	END IF
	391	! -----------------------------------------------
	392
	393
	394	DO il = 1, ncum
	395	IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
	396	il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
	397	IF (sij(il,i,j)>0.0) THEN
	398	rti = qnk(il) - ep(il, i)*clw(il, i)
	399	qmixmax(il) = qmix(sjmax(il))
	400	qmixmin(il) = qmix(sjmin(il))
	401	rmixmax(il) = rmix(sjmax(il))
	402	rmixmin(il) = rmix(sjmin(il))
	403	sqmrmax(il) = sjmax(il)*qmix(sjmax(il)) - rmix(sjmax(il))
	404	sqmrmin(il) = sjmin(il)*qmix(sjmin(il)) - rmix(sjmin(il))
	405
	406	ment(il, i, j) = abs(qmixmax(il)-qmixmin(il))*ment(il, i, j)
	407
	408	! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
	409	IF (abs(qmixmax(il)-qmixmin(il))>1.E-10) THEN
	410	sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/ &
	411	(qmixmax(il)-qmixmin(il))
	412	ELSE
	413	sigij(il, i, j) = 0.
	414	END IF
	415
	416	! -- Compute Qent, uent, vent according to the true mixing
	417	! fraction
	418	qent(il, i, j) = (1.-sigij(il,i,j))*rti + &
	419	sigij(il, i, j)*rr(il, i)
	420	uent(il, i, j) = (1.-sigij(il,i,j))*unk(il) + &
	421	sigij(il, i, j)*u(il, i)
	422	vent(il, i, j) = (1.-sigij(il,i,j))*vnk(il) + &
	423	sigij(il, i, j)*v(il, i)
	424
	425	! -- Compute liquid water static energy of mixed draughts
	426	! IF (j .GT. i) THEN
	427	! awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
	428	! awat=amax1(awat,0.0)
	429	! ELSE
	430	! awat = 0.
	431	! ENDIF
	432	! Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
	433	! : + Sigij(il,i,j)*H(il,i)
	434	! : + (LV(il,j)+(cpd-cpv)t(il,j))awat
	435	! IM 301008 beg
	436	hent(il, i, j) = (1.-sigij(il,i,j))*hp(il, i) + &
	437	sigij(il, i, j)*h(il, i)
	438
	439	elij(il, i, j) = qent(il, i, j) - rs(il, j)
	440	elij(il, i, j) = elij(il, i, j) + ((h(il,j)-hent(il,i, &
	441	j))rs(il,j)lv(il,j)/((cpd*(1.-qent(il,i,j))+ &
	442	qent(il,i,j)cpv)rrvt(il,j)t(il,j)))
	443	elij(il, i, j) = elij(il, i, j)/(1.+lv(il,j)lv(il,j)rs(il,j)/(( &
	444	cpd(1.-qent(il,i,j))+qent(il,i,j)cpv)rrvt(il,j)*t(il,j)))
	445
	446	elij(il, i, j) = max(elij(il,i,j), 0.)
	447
	448	elij(il, i, j) = min(elij(il,i,j), qent(il,i,j))
	449
	450	IF (j>i) THEN
	451	awat = elij(il, i, j) - (1.-ep(il,j))*clw(il, j)
	452	awat = amax1(awat, 0.0)
	453	ELSE
	454	awat = 0.
	455	END IF
	456
	457	! print
	458	! ,h(il,j)-hent(il,i,j),LV(il,j)rs(il,j)/(cpdrrvt(il,j)*
	459	! : t(il,j))
	460
	461	hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)t(il,j)) &
	462	awat
	463	! IM 301008 end
	464
	465	! print *,'mix : i,j,hent(il,i,j),sigij(il,i,j) ',
	466	! : i,j,hent(il,i,j),sigij(il,i,j)
	467
	468	! -- ASij is the integral of P(F) over the relevant F
	469	! interval
	470	asij(il) = asij(il) + abs(qmixmax(il)*(1.-sjmax(il))+rmixmax(il)- &
	471	qmixmin(il)*(1.-sjmin(il))-rmixmin(il))
	472
	473	END IF
	474	END IF
	475	END DO
	476	DO k = 1, ntra
	477	DO il = 1, ncum
	478	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
	479	1)) .AND. (j<=inb(il)) .AND. lwork(il)) THEN
	480	IF (sij(il,i,j)>0.0) THEN
	481	traent(il, i, j, k) = sigij(il, i, j)*tra(il, i, k) + &
	482	(1.-sigij(il,i,j))*tra(il, nk(il), k)
	483	END IF
	484	END IF
	485	END DO
	486	END DO
	487
	488	! -- If I=J (detrainement and entrainement at the same level), then
	489	! only the
	490	! -- adiabatic ascent part of the mixture is considered
	491	IF (i==j) THEN
	492	DO il = 1, ncum
	493	IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
	494	il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
	495	IF (sij(il,i,j)>0.0) THEN
	496	rti = qnk(il) - ep(il, i)*clw(il, i)
	497	! cc Ment(il,i,i) =
	498	! m(il,i)abs(Qmixmax(il)(1.-Sjmax(il))
	499	ment(il, i, i) = abs(qmixmax(il)*(1.-sjmax( &
	500	il))+rmixmax(il)-qmixmin(il)*(1.-sjmin(il))-rmixmin(il))
	501	qent(il, i, i) = rti
	502	uent(il, i, i) = unk(il)
	503	vent(il, i, i) = vnk(il)
	504	hent(il, i, i) = hp(il, i)
	505	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
	506	sigij(il, i, i) = 0.
	507	END IF
	508	END IF
	509	END DO
	510	DO k = 1, ntra
	511	DO il = 1, ncum
	512	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
	513	1)) .AND. (j<=inb(il)) .AND. lwork(il)) THEN
	514	IF (sij(il,i,j)>0.0) THEN
	515	traent(il, i, i, k) = tra(il, nk(il), k)
	516	END IF
	517	END IF
	518	END DO
	519	END DO
	520
	521	END IF
	522
	523	175 END DO
	524
	525	DO il = 1, ncum
	526	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
	527	asij(il) = amax1(1.0E-16, asij(il))
	528	asij(il) = 1.0/asij(il)
	529	csum(il, i) = 0.0
	530	END IF
	531	END DO
	532
	533	DO j = minorig, nl
	534	DO il = 1, ncum
	535	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
	536	il)-1) .AND. j<=inb(il)) THEN
	537	ment(il, i, j) = ment(il, i, j)*asij(il)
	538	END IF
	539	END DO
	540	END DO
	541
	542	DO j = minorig, nl
	543	DO il = 1, ncum
	544	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
	545	il)-1) .AND. j<=inb(il)) THEN
	546	csum(il, i) = csum(il, i) + ment(il, i, j)
	547	END IF
	548	END DO
	549	END DO
	550
	551	DO il = 1, ncum
	552	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) &
	553	THEN
	554	! cc : .and. csum(il,i).lt.m(il,i) ) then
	555	nent(il, i) = 0
	556	! cc ment(il,i,i)=m(il,i)
	557	ment(il, i, i) = 1.
	558	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
	559	uent(il, i, i) = unk(il)
	560	vent(il, i, i) = vnk(il)
	561	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
	562	sij(il, i, i) = 0.0
	563	END IF
	564	END DO ! il
	565
	566	DO j = 1, ntra
	567	DO il = 1, ncum
	568	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) &
	569	THEN
	570	! cc : .and. csum(il,i).lt.m(il,i) ) then
	571	traent(il, i, i, j) = tra(il, nk(il), j)
	572	END IF
	573	END DO
	574	END DO
	575
	576	789 END DO
	577
	578	RETURN
	579	END SUBROUTINE cv3p_mixing
	580

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

Download in other formats:

Original Format