/[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 47 by guez, Fri Jul 1 15:00:48 2011 UTC trunk/Sources/phylmd/CV3_routines/cv3_mixing.f revision 178 by guez, Fri Mar 11 18:47:26 2016 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 cvparam3    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 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 cv3_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      call zilch(asum, nloc*nd)
146        call zilch(csum, nloc*nd)
147        call zilch(csum, nloc*nd)
148    
149        do il=1, ncum
150         lwork(il) = .FALSE.         lwork(il) = .FALSE.
151        enddo      enddo
152    
153        DO 789 i=minorig+1,nl      DO i=minorig+1, nl
154    
155        num1=0         num1=0
156        do il=1,ncum         do il=1, ncum
157         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  
158         enddo         enddo
159  180   continue         if (num1 <= 0) cycle
160    
161        do 190 j=minorig,nl         do il=1, ncum
162         do il=1,ncum            if (i >= icb(il) .and. i <= inb(il)) then
163          if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &               lwork(il)=(nent(il, i) /= 0)
164           .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then               qp=rr(il, 1)-ep(il, i)*clw(il, i)
165           asum(il,i)=asum(il,i)+ment(il,i,j)               anum=h(il, i)-hp(il, i)-lv(il, i)*(qp-rs(il, i)) &
166           ment(il,i,j)=ment(il,i,j)*sig(il,j)                    +(cpv-cpd)*t(il, i)*(qp-rr(il, i))
167           bsum(il,i)=bsum(il,i)+ment(il,i,j)               denom=h(il, i)-hp(il, i)+lv(il, i)*(rr(il, i)-qp) &
168          endif                    +(cpd-cpv)*t(il, i)*(rr(il, i)-qp)
169         enddo               if(abs(denom) < 0.01)denom=0.01
170  190   continue               scrit(il)=anum/denom
171                 alt=qp-rs(il, i)+scrit(il)*(rr(il, i)-qp)
172                 if(scrit(il) <= 0.0.or.alt <= 0.0)scrit(il)=1.0
173                 smax(il)=0.0
174                 asij(il)=0.0
175              endif
176           end do
177    
178           do j=nl, minorig, -1
179    
180              num2=0
181              do il=1, ncum
182                 if (i >= icb(il) .and. i <= inb(il) .and. &
183                      j >= (icb(il)-1) .and. j <= inb(il) &
184                      .and. lwork(il)) num2=num2+1
185              enddo
186              if (num2 <= 0) cycle
187    
188              do il=1, ncum
189                 if (i >= icb(il) .and. i <= inb(il) .and. &
190                      j >= (icb(il)-1) .and. j <= inb(il) &
191                      .and. lwork(il)) then
192    
193                    if(sij(il, i, j) > 1.0e-16.and.sij(il, i, j) < 0.95)then
194                       wgh=1.0
195                       if(j > i)then
196                          sjmax=amax1(sij(il, i, j+1), smax(il))
197                          sjmax=amin1(sjmax, scrit(il))
198                          smax(il)=amax1(sij(il, i, j), smax(il))
199                          sjmin=amax1(sij(il, i, j-1), smax(il))
200                          sjmin=amin1(sjmin, scrit(il))
201                          if(sij(il, i, j) < (smax(il)-1.0e-16))wgh=0.0
202                          smid=amin1(sij(il, i, j), scrit(il))
203                       else
204                          sjmax=amax1(sij(il, i, j+1), scrit(il))
205                          smid=amax1(sij(il, i, j), scrit(il))
206                          sjmin=0.0
207                          if(j > 1)sjmin=sij(il, i, j-1)
208                          sjmin=amax1(sjmin, scrit(il))
209                       endif
210                       delp=abs(sjmax-smid)
211                       delm=abs(sjmin-smid)
212                       asij(il)=asij(il)+wgh*(delp+delm)
213                       ment(il, i, j)=ment(il, i, j)*(delp+delm)*wgh
214                    endif
215                 endif
216              end do
217    
218        do il=1,ncum         end do
219         if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then  
220          bsum(il,i)=amax1(bsum(il,i),1.0e-16)         do il=1, ncum
221          bsum(il,i)=1.0/bsum(il,i)            if (i >= icb(il).and.i <= inb(il).and.lwork(il)) then
222         endif               asij(il)=amax1(1.0e-16, asij(il))
223        enddo               asij(il)=1.0/asij(il)
224                 asum(il, i)=0.0
225        do 195 j=minorig,nl               bsum(il, i)=0.0
226         do il=1,ncum               csum(il, i)=0.0
227          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  
228         enddo         enddo
 195   continue  
229    
230        do 197 j=minorig,nl         do j=minorig, nl
231         do il=1,ncum            do il=1, ncum
232          if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &               if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
233           .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then                    .and. j >= (icb(il)-1) .and. j <= inb(il)) then
234           csum(il,i)=csum(il,i)+ment(il,i,j)                  ment(il, i, j)=ment(il, i, j)*asij(il)
235          endif               endif
236              enddo
237           end do
238    
239           do j=minorig, nl
240              do il=1, ncum
241                 if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
242                      .and. j >= (icb(il)-1) .and. j <= inb(il)) then
243                    asum(il, i)=asum(il, i)+ment(il, i, j)
244                    ment(il, i, j)=ment(il, i, j)*sig(il, j)
245                    bsum(il, i)=bsum(il, i)+ment(il, i, j)
246                 endif
247              enddo
248           end do
249    
250           do il=1, ncum
251              if (i >= icb(il).and.i <= inb(il).and.lwork(il)) then
252                 bsum(il, i)=amax1(bsum(il, i), 1.0e-16)
253                 bsum(il, i)=1.0/bsum(il, i)
254              endif
255         enddo         enddo
 197   continue  
256    
257        do il=1,ncum         do j=minorig, nl
258         if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &            do il=1, ncum
259             .and. csum(il,i).lt.m(il,i) ) then               if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
260          nent(il,i)=0                    .and. j >= (icb(il)-1) .and. j <= inb(il)) then
261          ment(il,i,i)=m(il,i)                  ment(il, i, j)=ment(il, i, j)*asum(il, i)*bsum(il, i)
262          qent(il,i,i)=rr(il,1)-ep(il,i)*clw(il,i)               endif
263          uent(il,i,i)=u(il,nk(il))            enddo
264          vent(il,i,i)=v(il,nk(il))         end do
265          elij(il,i,i)=clw(il,i)  
266  !MAF        sij(il,i,i)=1.0         do j=minorig, nl
267          sij(il,i,i)=0.0            do il=1, ncum
268         endif               if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
269        enddo ! il                    .and. j >= (icb(il)-1) .and. j <= inb(il)) then
270                    csum(il, i)=csum(il, i)+ment(il, i, j)
271  !      do j=1,ntra               endif
272  !       do il=1,ncum            enddo
273  !        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)         end do
274  !     :     .and. csum(il,i).lt.m(il,i) ) then  
275  !         traent(il,i,i,j)=tra(il,nk(il),j)         do il=1, ncum
276  !        endif            if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
277  !       enddo                 .and. csum(il, i) < m(il, i)) then
278  !      enddo               nent(il, i)=0
279  789   continue               ment(il, i, i)=m(il, i)
280  !               qent(il, i, i)=rr(il, 1)-ep(il, i)*clw(il, i)
281  ! MAF: renormalisation de MENT               uent(il, i, i)=u(il, nk(il))
282        do jm=1,nd               vent(il, i, i)=v(il, nk(il))
283          do im=1,nd               elij(il, i, i)=clw(il, i)
284            do il=1,ncum               !MAF sij(il, i, i)=1.0
285            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)  
286            endif            endif
287           end do         enddo ! il
288    
289        end DO
290    
291        ! MAF: renormalisation de MENT
292        do jm=1, nd
293           do im=1, nd
294              do il=1, ncum
295                 zm(il, im)=zm(il, im)+(1.-sij(il, im, jm))*ment(il, im, jm)
296              end do
297         end do         end do
298        end do      end do
299  !  
300        do jm=1,nd      do jm=1, nd
301         do im=1,nd         do im=1, nd
302          do 999 il=1,ncum            do il=1, ncum
303           qents(il,im,jm)=qent(il,im,jm)               if(zm(il, im) /= 0.) then
304           ments(il,im,jm)=ment(il,im,jm)                  ment(il, im, jm)=ment(il, im, jm)*m(il, im)/zm(il, im)
305  999     continue               endif
306              end do
307           end do
308        end do
309    
310        do jm=1, nd
311           do im=1, nd
312              do il=1, ncum
313                 qents(il, im, jm)=qent(il, im, jm)
314                 ments(il, im, jm)=ment(il, im, jm)
315              end do
316         enddo         enddo
317        enddo      enddo
318    
319      end SUBROUTINE cv3_mixing
320    
321        return  end module cv3_mixing_m
       end  
Legend:
Removed from v.47  
changed lines
  Added in v.178
  ViewVC Help
Powered by ViewVC 1.1.21