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