My Project
 All Classes Files Functions Variables Macros
cv3p_mixing.F
Go to the documentation of this file.
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  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) = qqa1*qmix1(f) + qqa2*qmix2(f)
81  rmix(f) = qqa1*rmix1(f) + qqa2*rmix2(f)
82 C
83  INTEGER, SAVE :: 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 !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)/(rrv*t(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)-cwat*bf2)
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  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
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 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  : * rrv*t(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  : * rrv*t(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)/(cpd*rrv*t(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
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