/[lmdze]/trunk/Sources/phylmd/CV30_routines/cv30_mixing.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/CV30_routines/cv30_mixing.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/phylmd/CV3_routines/cv3_mixing.f revision 82 by guez, Wed Mar 5 14:57:53 2014 UTC trunk/Sources/phylmd/CV3_routines/cv3_mixing.f revision 145 by guez, Tue Jun 16 15:23:29 2015 UTC
# Line 1  Line 1 
1    module cv3_mixing_m
2    
3        SUBROUTINE cv3_mixing(nloc,ncum,nd,na,ntra,icb,nk,inb &    implicit none
4                            ,ph,t,rr,rs,u,v,tra,h,lv,qnk &  
5                            ,hp,tv,tvp,ep,clw,m,sig &  contains
6           ,ment,qent,uent,vent, nent, sij,elij,ments,qents,traent)  
7              use cv3_param_m    SUBROUTINE cv3_mixing(nloc, ncum, nd, na, icb, nk, inb, t, rr, rs, u, v, h, &
8              use cvthermo         lv, hp, ep, clw, m, sig, ment, qent, uent, vent, nent, sij, elij, &
9        implicit none         ments, qents)
10        use cv3_param_m
11  !---------------------------------------------------------------------      use cvthermo
12  ! a faire:  
13  !     - changer rr(il,1) -> qnk(il)      !---------------------------------------------------------------------
14  !   - vectorisation de la partie normalisation des flux (do 789...)      ! a faire:
15  !---------------------------------------------------------------------      ! - changer rr(il, 1) -> qnk(il)
16        ! - vectorisation de la partie normalisation des flux (do 789)
17        !---------------------------------------------------------------------
18  ! inputs:  
19        integer, intent(in):: ncum, nd, na, ntra, nloc      ! inputs:
20        integer icb(nloc), inb(nloc), nk(nloc)      integer, intent(in):: ncum, nd, na, nloc
21        real sig(nloc,nd)      integer icb(nloc), inb(nloc), nk(nloc)
22        real qnk(nloc)      real sig(nloc, nd)
23        real ph(nloc,nd+1)      real t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
24        real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)      real u(nloc, nd), v(nloc, nd)
25        real u(nloc,nd), v(nloc,nd)      real lv(nloc, na), h(nloc, na), hp(nloc, na)
26        real tra(nloc,nd,ntra) ! input of convect3      real ep(nloc, na), clw(nloc, na)
27        real lv(nloc,na), h(nloc,na), hp(nloc,na)      real m(nloc, na) ! input of convect3
28        real tv(nloc,na), tvp(nloc,na), ep(nloc,na), clw(nloc,na)  
29        real m(nloc,na)        ! input of convect3      ! outputs:
30        real ment(nloc, na, na), qent(nloc, na, na)
31  ! outputs:      real uent(nloc, na, na), vent(nloc, na, na)
32        real ment(nloc,na,na), qent(nloc,na,na)      real sij(nloc, na, na), elij(nloc, na, na)
33        real uent(nloc,na,na), vent(nloc,na,na)      real ments(nloc, nd, nd), qents(nloc, nd, nd)
34        real sij(nloc,na,na), elij(nloc,na,na)      real sigij(nloc, nd, nd)
35        real traent(nloc,nd,nd,ntra)      integer nent(nloc, nd)
36        real ments(nloc,nd,nd), qents(nloc,nd,nd)  
37        real sigij(nloc,nd,nd)      ! local variables:
38        integer nent(nloc,nd)      integer i, j, k, il, im, jm
39        integer num1, num2
40  ! local variables:      real rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
41        integer i, j, k, il, im, jm      real alt, smid, sjmin, sjmax, delp, delm
42        integer num1, num2      real asij(nloc), smax(nloc), scrit(nloc)
43        real rti, bf2, anum, denom, dei, altem, cwat, stemp, qp      real asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
44        real alt, smid, sjmin, sjmax, delp, delm      real wgh
45        real asij(nloc), smax(nloc), scrit(nloc)      real zm(nloc, na)
46        real asum(nloc,nd),bsum(nloc,nd),csum(nloc,nd)      logical lwork(nloc)
47        real wgh  
48        real zm(nloc,na)      ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
49        logical lwork(nloc)  
50        do j=1, nl
51  !=====================================================================         do i=1, ncum
52  ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS            nent(i, j)=0
53  !=====================================================================            ! in convect3, m is computed in cv3_closure
54              ! ori m(i, 1)=0.0
55          do 361 j=1,nl         end do
56          do 360 i=1,ncum      end do
57            nent(i,j)=0  
58  ! in convect3, m is computed in cv3_closure      do j=1, nl
59  ! ori          m(i,1)=0.0         do k=1, nl
60   360    continue            do i=1, ncum
61   361    continue               qent(i, k, j)=rr(i, j)
62                 uent(i, k, j)=u(i, j)
63  ! ori      do 400 k=1,nlp               vent(i, k, j)=v(i, j)
64  ! ori       do 390 j=1,nlp               elij(i, k, j)=0.0
65        do 400 j=1,nl               !ym ment(i, k, j)=0.0
66         do 390 k=1,nl               !ym sij(i, k, j)=0.0
67            do 385 i=1,ncum            end do
68              qent(i,k,j)=rr(i,j)         end do
69              uent(i,k,j)=u(i,j)      end do
70              vent(i,k,j)=v(i,j)  
71              elij(i,k,j)=0.0      !ym
72  !ym            ment(i,k,j)=0.0      ment(1:ncum, 1:nd, 1:nd)=0.0
73  !ym            sij(i,k,j)=0.0      sij(1:ncum, 1:nd, 1:nd)=0.0
74   385      continue  
75   390    continue      zm(:, :)=0.
76   400  continue  
77        ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
78  !ym      ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
79        ment(1:ncum,1:nd,1:nd)=0.0      ! --- FRACTION (sij)
80        sij(1:ncum,1:nd,1:nd)=0.0  
81        do i=minorig+1, nl
82        zm(:,:)=0.  
83           do j=minorig, nl
84  !=====================================================================            do il=1, ncum
85  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING               if((i >= icb(il)).and.(i <= inb(il)).and. &
86  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING                    (j >= (icb(il)-1)).and.(j <= inb(il)))then
87  ! --- FRACTION (sij)  
88  !=====================================================================                  rti=rr(il, 1)-ep(il, i)*clw(il, i)
89                    bf2=1.+lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd)
90        do 750 i=minorig+1, nl                  anum=h(il, j)-hp(il, i)+(cpv-cpd)*t(il, j)*(rti-rr(il, j))
91                    denom=h(il, i)-hp(il, i)+(cpd-cpv)*(rr(il, i)-rti)*t(il, j)
92         do 710 j=minorig,nl                  dei=denom
93          do 700 il=1,ncum                  if(abs(dei) < 0.01)dei=0.01
94           if( (i.ge.icb(il)).and.(i.le.inb(il)).and. &                  sij(il, i, j)=anum/dei
95              (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then                  sij(il, i, i)=1.0
96                    altem=sij(il, i, j)*rr(il, i)+(1.-sij(il, i, j))*rti-rs(il, j)
97            rti=rr(il,1)-ep(il,i)*clw(il,i)                  altem=altem/bf2
98            bf2=1.+lv(il,j)*lv(il,j)*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd)                  cwat=clw(il, j)*(1.-ep(il, j))
99            anum=h(il,j)-hp(il,i)+(cpv-cpd)*t(il,j)*(rti-rr(il,j))                  stemp=sij(il, i, j)
100            denom=h(il,i)-hp(il,i)+(cpd-cpv)*(rr(il,i)-rti)*t(il,j)                  if((stemp < 0.0.or.stemp > 1.0.or.altem > cwat) &
101            dei=denom                       .and.j > i)then
102            if(abs(dei).lt.0.01)dei=0.01                     anum=anum-lv(il, j)*(rti-rs(il, j)-cwat*bf2)
103            sij(il,i,j)=anum/dei                     denom=denom+lv(il, j)*(rr(il, i)-rti)
104            sij(il,i,i)=1.0                     if(abs(denom) < 0.01)denom=0.01
105            altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j)                     sij(il, i, j)=anum/denom
106            altem=altem/bf2                     altem=sij(il, i, j)*rr(il, i)+(1.-sij(il, i, j))*rti-rs(il, j)
107            cwat=clw(il,j)*(1.-ep(il,j))                     altem=altem-(bf2-1.)*cwat
108            stemp=sij(il,i,j)                  end if
109            if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat) &                  if(sij(il, i, j) > 0.0.and.sij(il, i, j) < 0.95)then
110                         .and.j.gt.i)then                     qent(il, i, j)=sij(il, i, j)*rr(il, i)+(1.-sij(il, i, j))*rti
111             anum=anum-lv(il,j)*(rti-rs(il,j)-cwat*bf2)                     uent(il, i, j)=sij(il, i, j)*u(il, i)+(1.-sij(il, i, j))*u(il, nk(il))
112             denom=denom+lv(il,j)*(rr(il,i)-rti)                     vent(il, i, j)=sij(il, i, j)*v(il, i)+(1.-sij(il, i, j))*v(il, nk(il))
113             if(abs(denom).lt.0.01)denom=0.01                     elij(il, i, j)=altem
114             sij(il,i,j)=anum/denom                     elij(il, i, j)=amax1(0.0, elij(il, i, j))
115             altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j)                     ment(il, i, j)=m(il, i)/(1.-sij(il, i, j))
116             altem=altem-(bf2-1.)*cwat                     nent(il, i)=nent(il, i)+1
117                    end if
118                    sij(il, i, j)=amax1(0.0, sij(il, i, j))
119                    sij(il, i, j)=amin1(1.0, sij(il, i, j))
120                 endif ! new
121              end do
122           end do
123    
124           ! *** if no air can entrain at level i assume that updraft detrains ***
125           ! *** at that level and calculate detrained air flux and properties ***
126    
127           !@ do 170 i=icb(il), inb(il)
128    
129           do il=1, ncum
130              if ((i >= icb(il)).and.(i <= inb(il)).and.(nent(il, i) == 0)) then
131                 !@ if(nent(il, i) == 0)then
132                 ment(il, i, i)=m(il, i)
133                 qent(il, i, i)=rr(il, nk(il))-ep(il, i)*clw(il, i)
134                 uent(il, i, i)=u(il, nk(il))
135                 vent(il, i, i)=v(il, nk(il))
136                 elij(il, i, i)=clw(il, i)
137                 !MAF sij(il, i, i)=1.0
138                 sij(il, i, i)=0.0
139            end if            end if
140           if(sij(il,i,j).gt.0.0.and.sij(il,i,j).lt.0.95)then         end do
141            qent(il,i,j)=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti      end do
142            uent(il,i,j)=sij(il,i,j)*u(il,i)+(1.-sij(il,i,j))*u(il,nk(il))  
143            vent(il,i,j)=sij(il,i,j)*v(il,i)+(1.-sij(il,i,j))*v(il,nk(il))      do j=minorig, nl
144  !!!!      do k=1,ntra         do i=minorig, nl
145  !!!!      traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)            do il=1, ncum
146  !!!!     :      +(1.-sij(il,i,j))*tra(il,nk(il),k)               if ((j >= (icb(il)-1)).and.(j <= inb(il)) &
147  !!!!      end do                    .and.(i >= icb(il)).and.(i <= inb(il)))then
148            elij(il,i,j)=altem                  sigij(il, i, j)=sij(il, i, j)
149            elij(il,i,j)=amax1(0.0,elij(il,i,j))               endif
150            ment(il,i,j)=m(il,i)/(1.-sij(il,i,j))            end do
151            nent(il,i)=nent(il,i)+1         end do
152           end if      end do
153           sij(il,i,j)=amax1(0.0,sij(il,i,j))      !@ enddo
154           sij(il,i,j)=amin1(1.0,sij(il,i,j))  
155           endif ! new      !@170 continue
156   700   continue  
157   710  continue      ! --- NORMALIZE ENTRAINED AIR MASS FLUXES
158        ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
 !  
 !   ***   if no air can entrain at level i assume that updraft detrains  ***  
 !   ***   at that level and calculate detrained air flux and properties  ***  
 !  
   
 !@      do 170 i=icb(il),inb(il)  
   
       do 740 il=1,ncum  
       if ((i.ge.icb(il)).and.(i.le.inb(il)).and.(nent(il,i).eq.0)) then  
 !@      if(nent(il,i).eq.0)then  
       ment(il,i,i)=m(il,i)  
       qent(il,i,i)=rr(il,nk(il))-ep(il,i)*clw(il,i)  
       uent(il,i,i)=u(il,nk(il))  
       vent(il,i,i)=v(il,nk(il))  
       elij(il,i,i)=clw(il,i)  
 !MAF      sij(il,i,i)=1.0  
       sij(il,i,i)=0.0  
       end if  
  740  continue  
  750  continue  
   
       do 100 j=minorig,nl  
       do 101 i=minorig,nl  
       do 102 il=1,ncum  
       if ((j.ge.(icb(il)-1)).and.(j.le.inb(il)) &  
           .and.(i.ge.icb(il)).and.(i.le.inb(il)))then  
        sigij(il,i,j)=sij(il,i,j)  
       endif  
  102  continue  
  101  continue  
  100  continue  
 !@      enddo  
   
 !@170   continue  
   
 !=====================================================================  
 !   ---  NORMALIZE ENTRAINED AIR MASS FLUXES  
 !   ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING  
 !=====================================================================  
   
       call zilch(asum,nloc*nd)  
       call zilch(csum,nloc*nd)  
       call zilch(csum,nloc*nd)  
159    
160        do il=1,ncum      call zilch(asum, nloc*nd)
161        call zilch(csum, nloc*nd)
162        call zilch(csum, nloc*nd)
163    
164        do il=1, ncum
165         lwork(il) = .FALSE.         lwork(il) = .FALSE.
166        enddo      enddo
167    
168        DO 789 i=minorig+1,nl      DO i=minorig+1, nl
169    
170        num1=0         num1=0
171        do il=1,ncum         do il=1, ncum
172         if ( i.ge.icb(il) .and. i.le.inb(il) ) num1=num1+1            if (i >= icb(il) .and. i <= inb(il)) num1=num1+1
       enddo  
       if (num1.le.0) goto 789  
   
   
       do 781 il=1,ncum  
        if ( i.ge.icb(il) .and. i.le.inb(il) ) then  
         lwork(il)=(nent(il,i).ne.0)  
         qp=rr(il,1)-ep(il,i)*clw(il,i)  
         anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i)) &  
                  +(cpv-cpd)*t(il,i)*(qp-rr(il,i))  
         denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp) &  
                  +(cpd-cpv)*t(il,i)*(rr(il,i)-qp)  
         if(abs(denom).lt.0.01)denom=0.01  
         scrit(il)=anum/denom  
         alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp)  
         if(scrit(il).le.0.0.or.alt.le.0.0)scrit(il)=1.0  
         smax(il)=0.0  
         asij(il)=0.0  
        endif  
 781   continue  
   
       do 175 j=nl,minorig,-1  
   
       num2=0  
       do il=1,ncum  
        if ( i.ge.icb(il) .and. i.le.inb(il) .and. &  
             j.ge.(icb(il)-1) .and. j.le.inb(il)  &  
             .and. lwork(il) ) num2=num2+1  
       enddo  
       if (num2.le.0) goto 175  
   
       do 782 il=1,ncum  
       if ( i.ge.icb(il) .and. i.le.inb(il) .and. &  
             j.ge.(icb(il)-1) .and. j.le.inb(il)  &  
             .and. lwork(il) ) then  
   
        if(sij(il,i,j).gt.1.0e-16.and.sij(il,i,j).lt.0.95)then  
         wgh=1.0  
         if(j.gt.i)then  
          sjmax=amax1(sij(il,i,j+1),smax(il))  
          sjmax=amin1(sjmax,scrit(il))  
          smax(il)=amax1(sij(il,i,j),smax(il))  
          sjmin=amax1(sij(il,i,j-1),smax(il))  
          sjmin=amin1(sjmin,scrit(il))  
          if(sij(il,i,j).lt.(smax(il)-1.0e-16))wgh=0.0  
          smid=amin1(sij(il,i,j),scrit(il))  
         else  
          sjmax=amax1(sij(il,i,j+1),scrit(il))  
          smid=amax1(sij(il,i,j),scrit(il))  
          sjmin=0.0  
          if(j.gt.1)sjmin=sij(il,i,j-1)  
          sjmin=amax1(sjmin,scrit(il))  
         endif  
         delp=abs(sjmax-smid)  
         delm=abs(sjmin-smid)  
         asij(il)=asij(il)+wgh*(delp+delm)  
         ment(il,i,j)=ment(il,i,j)*(delp+delm)*wgh  
        endif  
       endif  
 782   continue  
   
 175   continue  
   
       do il=1,ncum  
        if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then  
         asij(il)=amax1(1.0e-16,asij(il))  
         asij(il)=1.0/asij(il)  
         asum(il,i)=0.0  
         bsum(il,i)=0.0  
         csum(il,i)=0.0  
        endif  
       enddo  
   
       do 180 j=minorig,nl  
        do il=1,ncum  
         if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &  
          .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then  
          ment(il,i,j)=ment(il,i,j)*asij(il)  
         endif  
173         enddo         enddo
174  180   continue         if (num1 <= 0) cycle
175    
176        do 190 j=minorig,nl         do il=1, ncum
177         do il=1,ncum            if (i >= icb(il) .and. i <= inb(il)) then
178          if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &               lwork(il)=(nent(il, i) /= 0)
179           .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then               qp=rr(il, 1)-ep(il, i)*clw(il, i)
180           asum(il,i)=asum(il,i)+ment(il,i,j)               anum=h(il, i)-hp(il, i)-lv(il, i)*(qp-rs(il, i)) &
181           ment(il,i,j)=ment(il,i,j)*sig(il,j)                    +(cpv-cpd)*t(il, i)*(qp-rr(il, i))
182           bsum(il,i)=bsum(il,i)+ment(il,i,j)               denom=h(il, i)-hp(il, i)+lv(il, i)*(rr(il, i)-qp) &
183          endif                    +(cpd-cpv)*t(il, i)*(rr(il, i)-qp)
184         enddo               if(abs(denom) < 0.01)denom=0.01
185  190   continue               scrit(il)=anum/denom
186                 alt=qp-rs(il, i)+scrit(il)*(rr(il, i)-qp)
187                 if(scrit(il) <= 0.0.or.alt <= 0.0)scrit(il)=1.0
188                 smax(il)=0.0
189                 asij(il)=0.0
190              endif
191           end do
192    
193           do j=nl, minorig, -1
194    
195        do il=1,ncum            num2=0
196         if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then            do il=1, ncum
197          bsum(il,i)=amax1(bsum(il,i),1.0e-16)               if (i >= icb(il) .and. i <= inb(il) .and. &
198          bsum(il,i)=1.0/bsum(il,i)                    j >= (icb(il)-1) .and. j <= inb(il) &
199         endif                    .and. lwork(il)) num2=num2+1
200        enddo            enddo
201              if (num2 <= 0) cycle
202        do 195 j=minorig,nl  
203         do il=1,ncum            do il=1, ncum
204          if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &               if (i >= icb(il) .and. i <= inb(il) .and. &
205           .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then                    j >= (icb(il)-1) .and. j <= inb(il) &
206           ment(il,i,j)=ment(il,i,j)*asum(il,i)*bsum(il,i)                    .and. lwork(il)) then
207          endif  
208                    if(sij(il, i, j) > 1.0e-16.and.sij(il, i, j) < 0.95)then
209                       wgh=1.0
210                       if(j > i)then
211                          sjmax=amax1(sij(il, i, j+1), smax(il))
212                          sjmax=amin1(sjmax, scrit(il))
213                          smax(il)=amax1(sij(il, i, j), smax(il))
214                          sjmin=amax1(sij(il, i, j-1), smax(il))
215                          sjmin=amin1(sjmin, scrit(il))
216                          if(sij(il, i, j) < (smax(il)-1.0e-16))wgh=0.0
217                          smid=amin1(sij(il, i, j), scrit(il))
218                       else
219                          sjmax=amax1(sij(il, i, j+1), scrit(il))
220                          smid=amax1(sij(il, i, j), scrit(il))
221                          sjmin=0.0
222                          if(j > 1)sjmin=sij(il, i, j-1)
223                          sjmin=amax1(sjmin, scrit(il))
224                       endif
225                       delp=abs(sjmax-smid)
226                       delm=abs(sjmin-smid)
227                       asij(il)=asij(il)+wgh*(delp+delm)
228                       ment(il, i, j)=ment(il, i, j)*(delp+delm)*wgh
229                    endif
230                 endif
231              end do
232    
233           end do
234    
235           do il=1, ncum
236              if (i >= icb(il).and.i <= inb(il).and.lwork(il)) then
237                 asij(il)=amax1(1.0e-16, asij(il))
238                 asij(il)=1.0/asij(il)
239                 asum(il, i)=0.0
240                 bsum(il, i)=0.0
241                 csum(il, i)=0.0
242              endif
243         enddo         enddo
 195   continue  
244    
245        do 197 j=minorig,nl         do j=minorig, nl
246         do il=1,ncum            do il=1, ncum
247          if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &               if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
248           .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then                    .and. j >= (icb(il)-1) .and. j <= inb(il)) then
249           csum(il,i)=csum(il,i)+ment(il,i,j)                  ment(il, i, j)=ment(il, i, j)*asij(il)
250          endif               endif
251              enddo
252           end do
253    
254           do j=minorig, nl
255              do il=1, ncum
256                 if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
257                      .and. j >= (icb(il)-1) .and. j <= inb(il)) then
258                    asum(il, i)=asum(il, i)+ment(il, i, j)
259                    ment(il, i, j)=ment(il, i, j)*sig(il, j)
260                    bsum(il, i)=bsum(il, i)+ment(il, i, j)
261                 endif
262              enddo
263           end do
264    
265           do il=1, ncum
266              if (i >= icb(il).and.i <= inb(il).and.lwork(il)) then
267                 bsum(il, i)=amax1(bsum(il, i), 1.0e-16)
268                 bsum(il, i)=1.0/bsum(il, i)
269              endif
270         enddo         enddo
 197   continue  
271    
272        do il=1,ncum         do j=minorig, nl
273         if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &            do il=1, ncum
274             .and. csum(il,i).lt.m(il,i) ) then               if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
275          nent(il,i)=0                    .and. j >= (icb(il)-1) .and. j <= inb(il)) then
276          ment(il,i,i)=m(il,i)                  ment(il, i, j)=ment(il, i, j)*asum(il, i)*bsum(il, i)
277          qent(il,i,i)=rr(il,1)-ep(il,i)*clw(il,i)               endif
278          uent(il,i,i)=u(il,nk(il))            enddo
279          vent(il,i,i)=v(il,nk(il))         end do
280          elij(il,i,i)=clw(il,i)  
281  !MAF        sij(il,i,i)=1.0         do j=minorig, nl
282          sij(il,i,i)=0.0            do il=1, ncum
283         endif               if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
284        enddo ! il                    .and. j >= (icb(il)-1) .and. j <= inb(il)) then
285                    csum(il, i)=csum(il, i)+ment(il, i, j)
286  789   continue               endif
287  !            enddo
288  ! MAF: renormalisation de MENT         end do
289        do jm=1,nd  
290          do im=1,nd         do il=1, ncum
291            do il=1,ncum            if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
292            zm(il,im)=zm(il,im)+(1.-sij(il,im,jm))*ment(il,im,jm)                 .and. csum(il, i) < m(il, i)) then
293           end do               nent(il, i)=0
294          end do               ment(il, i, i)=m(il, i)
295        end do               qent(il, i, i)=rr(il, 1)-ep(il, i)*clw(il, i)
296  !               uent(il, i, i)=u(il, nk(il))
297        do jm=1,nd               vent(il, i, i)=v(il, nk(il))
298          do im=1,nd               elij(il, i, i)=clw(il, i)
299            do il=1,ncum               !MAF sij(il, i, i)=1.0
300            if(zm(il,im).ne.0.) then               sij(il, i, i)=0.0
           ment(il,im,jm)=ment(il,im,jm)*m(il,im)/zm(il,im)  
301            endif            endif
302           end do         enddo ! il
303    
304        end DO
305    
306        ! MAF: renormalisation de MENT
307        do jm=1, nd
308           do im=1, nd
309              do il=1, ncum
310                 zm(il, im)=zm(il, im)+(1.-sij(il, im, jm))*ment(il, im, jm)
311              end do
312         end do         end do
313        end do      end do
314  !  
315        do jm=1,nd      do jm=1, nd
316         do im=1,nd         do im=1, nd
317          do 999 il=1,ncum            do il=1, ncum
318           qents(il,im,jm)=qent(il,im,jm)               if(zm(il, im) /= 0.) then
319           ments(il,im,jm)=ment(il,im,jm)                  ment(il, im, jm)=ment(il, im, jm)*m(il, im)/zm(il, im)
320  999     continue               endif
321              end do
322           end do
323        end do
324    
325        do jm=1, nd
326           do im=1, nd
327              do il=1, ncum
328                 qents(il, im, jm)=qent(il, im, jm)
329                 ments(il, im, jm)=ment(il, im, jm)
330              end do
331         enddo         enddo
332        enddo      enddo
333    
334      end SUBROUTINE cv3_mixing
335    
336        return  end module cv3_mixing_m
       end  
Legend:
Removed from v.82  
changed lines
  Added in v.145
  ViewVC Help
Powered by ViewVC 1.1.21