Context Navigation

cv3p_mixing.F @ 3241

Last change on this file since 3241 was 1494, checked in by Laurent Fairhead, 14 years ago

Optimization of routines from the new physics

Optimisation de routines de la nouvelle physique

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

File size: 18.9 KB

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

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

Download in other formats:

Original Format