Context Navigation

cv3p_mixing.F @ 3931

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

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

Download in other formats:

Original Format