source: LMDZ5/trunk/libf/phylmd/cv3p_mixing.F @ 1513

Last change on this file since 1513 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
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(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)
62c
63      REAL amxupcrit, df, ff
64      INTEGER nstep
65C
66c --   Mixing probability distribution functions
67c
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)
82C
83      INTEGER ifrst
84      DATA ifrst/0/
85C
86
87c=====================================================================
88c --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
89c=====================================================================
90c
91c -- Initialize mixing PDF coefficients
92      IF (ifrst .EQ. 0) THEN
93        ifrst = 1
94        Qcoef1max = Qcoef1(Fmax)
95        Qcoef2max = Qcoef2(Fmax)
96c
97      ENDIF
98c
99
100c ori        do 360 i=1,ncum*nlp
101        do 361 j=1,nl
102        do 360 i=1,ncum
103          nent(i,j)=0
104c in convect3, m is computed in cv3_closure
105c ori          m(i,1)=0.0
106 360    continue
107 361    continue
108
109c ori      do 400 k=1,nlp
110c 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
135c=====================================================================
136c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
137c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
138c --- FRACTION (sij)
139c=====================================================================
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)/(rrv*t(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)-cwat*bf2)
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
170ccc                 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
183c
184c   ***   if no air can entrain at level i assume that updraft detrains  ***
185c   ***   at that level and calculate detrained air flux and properties  ***
186c
187
188c@      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
193c@      if(nent(il,i).eq.0)then
194ccc      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
226c@      enddo
227
228c@170   continue
229
230c=====================================================================
231c   ---  NORMALIZE ENTRAINED AIR MASS FLUXES
232c   ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
233c=====================================================================
234
235      call zilch(csum,nloc*nd)
236
237      do il=1,ncum
238       lwork(il) = .FALSE.
239      enddo
240
241c---------------------------------------------------------------
242      DO 789 i=minorig+1,nl     !Loop on origin level "i"
243c---------------------------------------------------------------
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
251c
252cjyg1    Find maximum of SIJ for J>I, if any.
253c
254       Sx(:) = 0.
255c
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
274c
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)
287c
288cjyg1    Find new critical value Scrit2
289c        such that : Sij > Scrit2  => mixed draught will detrain at J<I
290c                    Sij < Scrit2  => mixed draught will detrain at J>I
291c
292       Scrit2 = min(Scrit(il),Sx(il))*max(0.,-signhpmh(il))
293     :         +Scrit(il)*max(0.,signhpmh(il))
294c
295       Scrit(il) = Scrit2
296c
297cjyg    Correction pour la nouvelle logique; la correction pour ALT
298c       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
301C
302        smax(il)=0.0
303        asij(il)=0.0
304       Sup(il)=0.     ! upper S-value reached by descending draughts
305       endif
306781   continue
307
308c---------------------------------------------------------------
309      do 175 j=minorig,nl   !Loop on destination level "j"
310c---------------------------------------------------------------
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
320c -----------------------------------------------
321         IF (j .GT. i) THEN
322c -----------------------------------------------
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
344782   continue
345c -----------------------------------------------
346         ELSE IF (j .EQ. i) THEN
347c -----------------------------------------------
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.
360c
361c-           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
369783   continue
370c -----------------------------------------------
371         ELSE IF ( j .LT. i) THEN
372c -----------------------------------------------
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
396784   continue
397c -----------------------------------------------
398         END IF
399c -----------------------------------------------
400c
401c
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))
414c
415         Ment(il,i,j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il,i,j)
416c
417c    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
424c
425c --    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)
432c
433c--     Compute liquid water static energy of mixed draughts
434c       IF (j .GT. i) THEN
435c        awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
436c        awat=amax1(awat,0.0)
437c       ELSE
438c        awat = 0.
439c       ENDIF
440c       Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
441c    :         + Sigij(il,i,j)*H(il,i)
442c    :         + (LV(il,j)+(cpd-cpv)*t(il,j))*awat
443cIM 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     :    * rrv*t(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     :    * rrv*t(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
468c        print *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)*
469c    :         t(il,j))
470
471        Hent(il,i,j) =  Hent(il,i,j)+(LV(il,j)+(cpd-cpv)*t(il,j))
472     :         * awat
473cIM 301008 end
474c
475c      print *,'mix : i,j,hent(il,i,j),sigij(il,i,j) ',
476c     :               i,j,hent(il,i,j),sigij(il,i,j)
477c
478c --      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))
482c
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
498c
499c --    If I=J (detrainement and entrainement at the same level), then only the
500c --    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)
508ccc          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
532c
533        ENDIF
534c
535175   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
552180   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
561197   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
566ccc     :     .and. csum(il,i).lt.m(il,i) ) then
567        nent(il,i)=0
568ccc        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
582ccc     :     .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
587c
588789   continue
589c
590      return
591      end
592
Note: See TracBrowser for help on using the repository browser.