Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

cv3p_mixing.F90 @ 17

Last change on this file since 17 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 21.4 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	& ,sigij,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 sigij(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 hent(nloc,nd,nd)
44	integer nent(nloc,nd)
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)
54	real Sx(nloc), Scrit2
55	real smid(nloc), sjmin(nloc), sjmax(nloc)
56	real Sbef(nloc), Sup(nloc), Smin(nloc)
57	real asij(nloc), smax(nloc), scrit(nloc)
58	real sij(nloc,nd,nd)
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	& +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, SAVE :: ifrst
84	DATA ifrst/0/
85	!$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	!AC! ment(i,k,j)=0.0
121	!AC! sij(i,k,j)=0.0
122	385 continue
123	390 continue
124	400 continue
125
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
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
247	!c---------------------------------------------------------------
248	DO 789 i=minorig+1,nl !Loop on origin level "i"
249	!c---------------------------------------------------------------
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
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
268
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
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
294	!cjyg1 Find new critical value Scrit2
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
298	Scrit2 = min(Scrit(il),Sx(il))*max(0.,-signhpmh(il)) &
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
314	!c---------------------------------------------------------------
315	do 175 j=minorig,nl !Loop on destination level "j"
316	!c---------------------------------------------------------------
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	& 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	& Sij(il,i,j) , &
342	& Scrit(il) )
343	Sjmin(il)=max( (Sbef(il)+Sij(il,i,j))/2. , &
344	& 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	& *max(0.,-signhpmh(il)) &
362	& +min((Sij(il,i,j+1)+Smid(il))/2.,Scrit(il)) &
363	& *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	& 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	& Sij(il,i,j) )
392	Sjmax(il) = max(Sjmax(il),Scrit(il))
393	Sjmin(il) = min( (Sbef(il)+Sij(il,i,j))/2. , &
394	& 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	& 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
413	rti=qnk(il)-ep(il,i)*clw(il,i)
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
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	!IM 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
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
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	!IM 301008 end
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	& + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) &
487	& -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	& +Rmixmax(il) &
517	& -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