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