/[lmdze]/trunk/phylmd/CV30_routines/cv30_mixing.f90
ViewVC logotype

Diff of /trunk/phylmd/CV30_routines/cv30_mixing.f90

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

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