source: LMDZ5/branches/LMDZ5-DOFOCO/libf/phylmd/cv3p_mixing.F @ 2045

Last change on this file since 2045 was 1742, checked in by idelkadi, 12 years ago

1- Inclusion des developpements de la these de Romain Pilon sur le
lessivage des aerosols :

a/ par les pluies convectives (modifs cv30_routines et cv3_routines pour

sortir les champs nécessaires au calcul off-line ; modif cvltr)

b/ par les pluies stratiformes (modifs phytrac et introduction

lsc_scav).

2- Choix entre plusieurs schemas pour les pluies stratiformes, commande
par iflag_lscav.

3- Quelques corrections dans la convection "Nouvelle Physique" pour
assurer la conservation des traceurs (cv3p1_mixing et cva_driver) (travail
de Robin Locatelli).
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.0 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
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 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
46c 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)
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, SAVE :: ifrst
84      DATA ifrst/0/
85c$OMP THREADPRIVATE(ifrst)
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!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
141c=====================================================================
142c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
143c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
144c --- FRACTION (sij)
145c=====================================================================
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
176ccc                 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
189c
190c   ***   if no air can entrain at level i assume that updraft detrains  ***
191c   ***   at that level and calculate detrained air flux and properties  ***
192c
193
194c@      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
199c@      if(nent(il,i).eq.0)then
200ccc      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
232c@      enddo
233
234c@170   continue
235
236c=====================================================================
237c   ---  NORMALIZE ENTRAINED AIR MASS FLUXES
238c   ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
239c=====================================================================
240
241      call zilch(csum,nloc*nd)
242
243      do il=1,ncum
244       lwork(il) = .FALSE.
245      enddo
246
247c---------------------------------------------------------------
248      DO 789 i=minorig+1,nl     !Loop on origin level "i"
249c---------------------------------------------------------------
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
257c
258cjyg1    Find maximum of SIJ for J>I, if any.
259c
260       Sx(:) = 0.
261c
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
280c
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)
293c
294cjyg1    Find new critical value Scrit2
295c        such that : Sij > Scrit2  => mixed draught will detrain at J<I
296c                    Sij < Scrit2  => mixed draught will detrain at J>I
297c
298       Scrit2 = min(Scrit(il),Sx(il))*max(0.,-signhpmh(il))
299     :         +Scrit(il)*max(0.,signhpmh(il))
300c
301       Scrit(il) = Scrit2
302c
303cjyg    Correction pour la nouvelle logique; la correction pour ALT
304c       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
307C
308        smax(il)=0.0
309        asij(il)=0.0
310       Sup(il)=0.     ! upper S-value reached by descending draughts
311       endif
312781   continue
313
314c---------------------------------------------------------------
315      do 175 j=minorig,nl   !Loop on destination level "j"
316c---------------------------------------------------------------
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
326c -----------------------------------------------
327         IF (j .GT. i) THEN
328c -----------------------------------------------
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
350782   continue
351c -----------------------------------------------
352         ELSE IF (j .EQ. i) THEN
353c -----------------------------------------------
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.
366c
367c-           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
375783   continue
376c -----------------------------------------------
377         ELSE IF ( j .LT. i) THEN
378c -----------------------------------------------
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
402784   continue
403c -----------------------------------------------
404         END IF
405c -----------------------------------------------
406c
407c
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))
420c
421         Ment(il,i,j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il,i,j)
422c
423c    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
430c
431c --    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)
438c
439c--     Compute liquid water static energy of mixed draughts
440c       IF (j .GT. i) THEN
441c        awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
442c        awat=amax1(awat,0.0)
443c       ELSE
444c        awat = 0.
445c       ENDIF
446c       Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
447c    :         + Sigij(il,i,j)*H(il,i)
448c    :         + (LV(il,j)+(cpd-cpv)*t(il,j))*awat
449cIM 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
474c        print *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)*
475c    :         t(il,j))
476
477        Hent(il,i,j) =  Hent(il,i,j)+(LV(il,j)+(cpd-cpv)*t(il,j))
478     :         * awat
479cIM 301008 end
480c
481c      print *,'mix : i,j,hent(il,i,j),sigij(il,i,j) ',
482c     :               i,j,hent(il,i,j),sigij(il,i,j)
483c
484c --      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))
488c
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
504c
505c --    If I=J (detrainement and entrainement at the same level), then only the
506c --    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)
514ccc          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
538c
539        ENDIF
540c
541175   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
558180   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
567197   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
572ccc     :     .and. csum(il,i).lt.m(il,i) ) then
573        nent(il,i)=0
574ccc        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
588ccc     :     .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
593c
594789   continue
595c
596      return
597      end
598
Note: See TracBrowser for help on using the repository browser.