Context Navigation

cv3p_mixing.F @ 1784

Last change on this file since 1784 was 1742, checked in by idelkadi, 12 years ago

1- Inclusion des developpements de la these de Romain Pilon sur le
lessivage des aerosols :

a/ par les pluies convectives (modifs cv30_routines et cv3_routines pour

sortir les champs nécessaires au calcul off-line ; modif cvltr)

b/ par les pluies stratiformes (modifs phytrac et introduction

lsc_scav).

2- Choix entre plusieurs schemas pour les pluies stratiformes, commande
par iflag_lscav.

3- Quelques corrections dans la convection "Nouvelle Physique" pour
assurer la conservation des traceurs (cv3p1_mixing et cva_driver) (travail
de Robin Locatelli).

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 19.0 KB

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

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

Download in other formats:

Original Format