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

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

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

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