/[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 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 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 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 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  
   
 !  
 !   ***   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)  
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  789   continue               endif
271  !            enddo
272  ! MAF: renormalisation de MENT         end do
273        do jm=1,nd  
274          do im=1,nd         do il=1, ncum
275            do il=1,ncum            if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
276            zm(il,im)=zm(il,im)+(1.-sij(il,im,jm))*ment(il,im,jm)                 .and. csum(il, i) < m(il, i)) then
277           end do               nent(il, i)=0
278          end do               ment(il, i, i)=m(il, i)
279        end do               qent(il, i, i)=rr(il, 1)-ep(il, i)*clw(il, i)
280  !               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            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)  
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.69  
changed lines
  Added in v.186
  ViewVC Help
Powered by ViewVC 1.1.21