source: trunk/libf/phylmd/cv3p_mixing.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

Import initial LMDZ5
File size: 18.3 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
16c
17#include "cvthermo.h"
18#include "cv3param.h"
19#include "YOMCST2.h"
20
21c 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
35c 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
47c 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, Scrit2
56      integer Jx
57      real smid(nloc), sjmin(nloc), sjmax(nloc)
58      real Sbef(nloc), Sup(nloc), Smin(nloc)
59      real asij(nloc), smax(nloc), scrit(nloc)
60      real csum(nloc,nd)
61      real awat
62      logical lwork(nloc)
63c
64      REAL amxupcrit, df, ff
65      INTEGER nstep
66C
67c --   Mixing probability distribution functions
68c
69      real Qcoef1,Qcoef2,QFF,QFFF,Qmix,Rmix,Qmix1,Rmix1,Qmix2,Rmix2,F
70      Qcoef1(F) = tanh(F/gammas)
71      Qcoef2(F) = ( tanh(F/gammas) + gammas *
72     $            log(cosh((1.- F)/gammas)/cosh(F/gammas)))
73      QFF(F) = Max(Min(F,1.),0.)
74      QFFF(F) = Min(QFF(F),scut)
75      Qmix1(F) = ( tanh((QFF(F) - Fmax)/gammas)+Qcoef1max )/
76     $           Qcoef2max
77      Rmix1(F) = ( gammas*log(cosh((QFF(F)-Fmax)/gammas))
78     1             +QFF(F)*Qcoef1max ) / Qcoef2max
79      Qmix2(F) = -Log(1.-QFFF(F))/scut
80      Rmix2(F) = (QFFF(F)+(1.-QFF(F))*Log(1.-QFFF(F)))/scut
81      Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F)
82      Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F)
83C
84      INTEGER ifrst
85      DATA ifrst/0/
86C
87
88c=====================================================================
89c --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
90c=====================================================================
91c
92c -- Initialize mixing PDF coefficients
93      IF (ifrst .EQ. 0) THEN
94        ifrst = 1
95        Qcoef1max = Qcoef1(Fmax)
96        Qcoef2max = Qcoef2(Fmax)
97c
98      ENDIF
99c
100
101c ori        do 360 i=1,ncum*nlp
102        do 361 j=1,nl
103        do 360 i=1,ncum
104          nent(i,j)=0
105c in convect3, m is computed in cv3_closure
106c ori          m(i,1)=0.0
107 360    continue
108 361    continue
109
110c ori      do 400 k=1,nlp
111c 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
136c=====================================================================
137c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
138c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
139c --- FRACTION (sij)
140c=====================================================================
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)/(rrv*t(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)-cwat*bf2)
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
171ccc                 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
184c
185c   ***   if no air can entrain at level i assume that updraft detrains  ***
186c   ***   at that level and calculate detrained air flux and properties  ***
187c
188
189c@      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
194c@      if(nent(il,i).eq.0)then
195ccc      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
227c@      enddo
228
229c@170   continue
230
231c=====================================================================
232c   ---  NORMALIZE ENTRAINED AIR MASS FLUXES
233c   ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
234c=====================================================================
235
236      call zilch(csum,nloc*nd)
237
238      do il=1,ncum
239       lwork(il) = .FALSE.
240      enddo
241
242      DO 789 i=minorig+1,nl
243
244      num1=0
245      do il=1,ncum
246       if ( i.ge.icb(il) .and. i.le.inb(il) ) num1=num1+1
247      enddo
248      if (num1.le.0) goto 789
249
250
251      do 781 il=1,ncum
252       if ( i.ge.icb(il) .and. i.le.inb(il) ) then
253        lwork(il)=(nent(il,i).ne.0)
254        Signhpmh(il) = sign(1.,hp(il,i)-h(il,i))
255        qp=qnk(il)-ep(il,i)*clw(il,i)
256        anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i))
257     :           +(cpv-cpd)*t(il,i)*(qp-rr(il,i))
258        denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp)
259     :           +(cpd-cpv)*t(il,i)*(rr(il,i)-qp)
260        if(abs(denom).lt.0.01)denom=0.01
261        scrit(il)=min(anum/denom,1.)
262        alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp)
263c
264cjyg1    Find maximum of SIJ for J>I, if any, and new critical value Scrit2
265c        such that : Sij > Scrit2  => mixed draught will detrain at J<I
266c                    Sij < Scrit2  => mixed draught will detrain at J>I
267c
268       Sx = 0.
269       Jx = 0.
270       Sbef(il) = max(0.,signhpmh(il))
271       DO j = i+1,inb(il)
272         IF (Sbef(il) .LT. Sij(il,i,j)) THEN
273           Sx = max(Sij(il,i,j),Sx)
274           Jx = J
275         ENDIF
276         Sbef(il) = Sij(il,i,j)
277       ENDDO
278c
279       Scrit2 = min(Scrit(il),Sx)*max(0.,-signhpmh(il))
280     :         +Scrit(il)*max(0.,signhpmh(il))
281c
282       Scrit(il) = Scrit2
283c
284cjyg    Correction pour la nouvelle logique; la correction pour ALT
285c       est un peu au hazard
286       if(scrit(il).le.0.0)scrit(il)=0.0
287       if(alt.le.0.0) scrit(il)=1.0
288C
289        smax(il)=0.0
290        asij(il)=0.0
291       Sup(il)=0.     ! upper S-value reached by descending draughts
292       endif
293781   continue
294
295      do 175 j=minorig,nl
296
297      num2=0
298      do il=1,ncum
299       if ( i.ge.icb(il) .and. i.le.inb(il) .and.
300     :      j.ge.(icb(il)-1) .and. j.le.inb(il)
301     :      .and. lwork(il) ) num2=num2+1
302      enddo
303      if (num2.le.0) goto 175
304
305c -----------------------------------------------
306         IF (j .GT. i) THEN
307c -----------------------------------------------
308      do 782 il=1,ncum
309      if ( i.ge.icb(il) .and. i.le.inb(il) .and.
310     :      j.ge.(icb(il)-1) .and. j.le.inb(il)
311     :      .and. lwork(il) ) then
312       if(sij(il,i,j).gt.0.0)then
313           Smid(il)=min(Sij(il,i,j),Scrit(il))
314           Sjmax(il)=Smid(il)
315           Sjmin(il)=Smid(il)
316           IF (Smid(il) .LT. Smin(il) .AND.
317     1                         Sij(il,i,j+1) .LT. Smid(il)) THEN
318             Smin(il)=Smid(il)
319             Sjmax(il)=min( (Sij(il,i,j+1)+Sij(il,i,j))/2. ,
320     1                  Sij(il,i,j) ,
321     1                  Scrit(il) )
322             Sjmin(il)=max( (Sbef(il)+Sij(il,i,j))/2. ,
323     1                  Sij(il,i,j) )
324             Sjmin(il)=min(Sjmin(il),Scrit(il))
325             Sbef(il) = Sij(il,i,j)
326           ENDIF
327      endif
328      endif
329782   continue
330c -----------------------------------------------
331         ELSE IF (j .EQ. i) THEN
332c -----------------------------------------------
333      do 783 il=1,ncum
334      if ( i.ge.icb(il) .and. i.le.inb(il) .and.
335     :      j.ge.(icb(il)-1) .and. j.le.inb(il)
336     :      .and. lwork(il) ) then
337       if(sij(il,i,j).gt.0.0)then
338           Smid(il) = 1.
339           Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2.,Scrit(il))
340     1                                         *max(0.,-signhpmh(il))
341     1            +min((Sij(il,i,j+1)+Smid(il))/2.,Scrit(il))
342     1                                         *max(0., signhpmh(il))
343           Sjmin(il) = max(Sjmin(il),Sup(il))
344           Sjmax(il) = 1.
345c
346c-           preparation des variables Scrit, Smin et Sbef pour la partie j>i
347           Scrit(il) = min(Sjmin(il),Sjmax(il),Scrit(il))
348
349           Smin(il) = 1.
350           Sbef(il) = max(0.,signhpmh(il))
351           Supmax(il,i) = sign(Scrit(il),-signhpmh(il))
352      endif
353      endif
354783   continue
355c -----------------------------------------------
356         ELSE IF ( j .LT. i) THEN
357c -----------------------------------------------
358      do 784 il=1,ncum
359      if ( i.ge.icb(il) .and. i.le.inb(il) .and.
360     :      j.ge.(icb(il)-1) .and. j.le.inb(il)
361     :      .and. lwork(il) ) then
362       if(sij(il,i,j).gt.0.0)then
363           Smid(il)=max(Sij(il,i,j),Scrit(il))
364           Sjmax(il) = Smid(il)
365           Sjmin(il) = Smid(il)
366           IF (Smid(il) .GT. Smax(il) .AND.
367     1                          Sij(il,i,j+1) .GT. Smid(il)) THEN
368             Smax(il) = Smid(il)
369             Sjmax(il) = max( (Sij(il,i,j+1)+Sij(il,i,j))/2. ,
370     1                                               Sij(il,i,j) )
371             Sjmax(il) = max(Sjmax(il),Scrit(il))
372             Sjmin(il) = min( (Sbef(il)+Sij(il,i,j))/2. ,
373     1                                               Sij(il,i,j) )
374             Sjmin(il) = max(Sjmin(il),Scrit(il))
375             Sbef(il) = Sij(il,i,j)
376           ENDIF
377          IF (abs(Sjmin(il)-Sjmax(il)) .GT. 1.e-10) Sup(il)=
378     1                            max(Sjmin(il),Sjmax(il),Sup(il))
379      endif
380      endif
381784   continue
382c -----------------------------------------------
383         END IF
384c -----------------------------------------------
385c
386c
387      do il=1,ncum
388      if ( i.ge.icb(il) .and. i.le.inb(il) .and.
389     :      j.ge.(icb(il)-1) .and. j.le.inb(il)
390     :      .and. lwork(il) ) then
391       if(sij(il,i,j).gt.0.0)then
392         rti=qnk(il)-ep(il,i)*clw(il,i)
393         Qmixmax(il)=Qmix(Sjmax(il))
394         Qmixmin(il)=Qmix(Sjmin(il))
395         Rmixmax(il)=Rmix(Sjmax(il))
396         Rmixmin(il)=Rmix(Sjmin(il))
397         SQmRmax(il)= Sjmax(il)*Qmix(Sjmax(il))-Rmix(Sjmax(il))
398         SQmRmin(il)= Sjmin(il)*Qmix(Sjmin(il))-Rmix(Sjmin(il))
399c
400         Ment(il,i,j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il,i,j)
401c
402c    Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
403         IF (abs(Qmixmax(il)-Qmixmin(il)) .GT. 1.e-10) THEN
404           Sigij(il,i,j) =
405     :           (SQmRmax(il)-SQmRmin(il))/(Qmixmax(il)-Qmixmin(il))
406         ELSE
407           Sigij(il,i,j) = 0.
408         ENDIF
409c
410c --    Compute Qent, uent, vent according to the true mixing fraction
411        Qent(il,i,j) = (1.-Sigij(il,i,j))*rti
412     :               + Sigij(il,i,j)*rr(il,i)
413        uent(il,i,j) = (1.-Sigij(il,i,j))*unk(il)
414     :               + Sigij(il,i,j)*u(il,i)
415        vent(il,i,j) = (1.-Sigij(il,i,j))*vnk(il)
416     :               + Sigij(il,i,j)*v(il,i)
417c
418c--     Compute liquid water static energy of mixed draughts
419c       IF (j .GT. i) THEN
420c        awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
421c        awat=amax1(awat,0.0)
422c       ELSE
423c        awat = 0.
424c       ENDIF
425c       Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
426c    :         + Sigij(il,i,j)*H(il,i)
427c    :         + (LV(il,j)+(cpd-cpv)*t(il,j))*awat
428cIM 301008 beg
429        Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
430     :         + Sigij(il,i,j)*H(il,i)
431
432        Elij(il,i,j) = Qent(il,i,j)-rs(il,j)
433        Elij(il,i,j) = Elij(il,i,j)
434     :    + ((h(il,j)-Hent(il,i,j))*rs(il,j)*LV(il,j)
435     :    / ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)
436     :    * rrv*t(il,j)*t(il,j)))
437        Elij(il,i,j) = Elij(il,i,j)
438     :    / (1.+LV(il,j)*LV(il,j)*rs(il,j)
439     :    / ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)
440     :    * rrv*t(il,j)*t(il,j)))
441
442        Elij(il,i,j) = max(elij(il,i,j),0.)
443
444        Elij(il,i,j) = min(elij(il,i,j),Qent(il,i,j))
445
446        IF (j .GT. i) THEN
447         awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
448         awat=amax1(awat,0.0)
449        ELSE
450         awat = 0.
451        ENDIF
452
453c        print *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)*
454c    :         t(il,j))
455
456        Hent(il,i,j) =  Hent(il,i,j)+(LV(il,j)+(cpd-cpv)*t(il,j))
457     :         * awat
458cIM 301008 end
459c
460c      print *,'mix : i,j,hent(il,i,j),sigij(il,i,j) ',
461c     :               i,j,hent(il,i,j),sigij(il,i,j)
462c
463c --      ASij is the integral of P(F) over the relevant F interval
464         ASij(il) = ASij(il)
465     1               + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il)
466     1                    -Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
467c
468      endif
469      endif
470      enddo
471       do k=1,ntra
472         do il=1,ncum
473          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
474     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il))
475     :      .and. lwork(il) ) then
476          if(sij(il,i,j).gt.0.0)then
477            traent(il,i,j,k)=sigij(il,i,j)*tra(il,i,k)
478     :            +(1.-sigij(il,i,j))*tra(il,nk(il),k)
479          endif
480          endif
481         enddo
482       enddo
483c
484c --    If I=J (detrainement and entrainement at the same level), then only the
485c --    adiabatic ascent part of the mixture is considered
486        IF (I .EQ. J) THEN
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          rti=qnk(il)-ep(il,i)*clw(il,i)
493ccc          Ment(il,i,i) = m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il))
494          Ment(il,i,i) = abs(Qmixmax(il)*(1.-Sjmax(il))
495     1                    +Rmixmax(il)
496     1                    -Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
497          Qent(il,i,i) = rti
498          uent(il,i,i) = unk(il)
499          vent(il,i,i) = vnk(il)
500          Hent(il,i,i) = hp(il,i)
501          Elij(il,i,i) = clw(il,i)*(1.-ep(il,i))
502          Sigij(il,i,i) = 0.
503      endif
504      endif
505      enddo
506       do k=1,ntra
507         do il=1,ncum
508          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
509     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il))
510     :      .and. lwork(il) ) then
511          if(sij(il,i,j).gt.0.0)then
512            traent(il,i,i,k)=tra(il,nk(il),k)
513          endif
514          endif
515         enddo
516       enddo
517c
518        ENDIF
519c
520175   continue
521
522      do il=1,ncum
523       if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
524        asij(il)=amax1(1.0e-16,asij(il))
525        asij(il)=1.0/asij(il)
526        csum(il,i)=0.0
527       endif
528      enddo
529
530      do 180 j=minorig,nl
531       do il=1,ncum
532        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
533     :   .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
534         ment(il,i,j)=ment(il,i,j)*asij(il)
535        endif
536       enddo
537180   continue
538
539      do 197 j=minorig,nl
540       do il=1,ncum
541        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
542     :   .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
543         csum(il,i)=csum(il,i)+ment(il,i,j)
544        endif
545       enddo
546197   continue
547
548      do il=1,ncum
549       if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
550     :     .and. csum(il,i).lt.1. ) then
551ccc     :     .and. csum(il,i).lt.m(il,i) ) then
552        nent(il,i)=0
553ccc        ment(il,i,i)=m(il,i)
554        ment(il,i,i)=1.
555        qent(il,i,i)=qnk(il)-ep(il,i)*clw(il,i)
556        uent(il,i,i)=unk(il)
557        vent(il,i,i)=vnk(il)
558        elij(il,i,i)=clw(il,i)*(1.-ep(il,i))
559        sij(il,i,i)=0.0
560       endif
561      enddo ! il
562
563      do j=1,ntra
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
567ccc     :     .and. csum(il,i).lt.m(il,i) ) then
568         traent(il,i,i,j)=tra(il,nk(il),j)
569        endif
570       enddo
571      enddo
572c
573789   continue
574c
575      return
576      end
577
Note: See TracBrowser for help on using the repository browser.