/[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/CV30_routines/cv30_mixing.f revision 198 by guez, Tue May 31 16:17:35 2016 UTC
# Line 1  Line 1 
1    module cv30_mixing_m
2    
3        SUBROUTINE cv3_mixing(nloc,ncum,nd,na,ntra,icb,nk,inb &    implicit none
4                            ,ph,t,rr,rs,u,v,tra,h,lv,qnk &  
5                            ,hp,tv,tvp,ep,clw,m,sig &  contains
6           ,ment,qent,uent,vent, nent, sij,elij,ments,qents,traent)  
7              use cvparam3    SUBROUTINE cv30_mixing(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 cv_thermo_m, only: cpd, cpv, rrv
18  ! inputs:      USE dimphy, ONLY: klev, klon
19        integer 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 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  ! ori        do 360 i=1,ncum*nlp            nent(i, j) = 0
56          do 361 j=1,nl         end do
57          do 360 i=1,ncum      end do
58            nent(i,j)=0  
59  ! in convect3, m is computed in cv3_closure      do j = 1, nl
60  ! ori          m(i,1)=0.0         do k = 1, nl
61   360    continue            do i = 1, ncum
62   361    continue               qent(i, k, j) = rr(i, j)
63                 uent(i, k, j) = u(i, j)
64  ! ori      do 400 k=1,nlp               vent(i, k, j) = v(i, j)
65  ! ori       do 390 j=1,nlp               elij(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)      ment(1:ncum, 1:klev, 1:klev) = 0.0
71              vent(i,k,j)=v(i,j)      sij(1:ncum, 1:klev, 1:klev) = 0.0
72              elij(i,k,j)=0.0  
73  !ym            ment(i,k,j)=0.0      zm(:, :) = 0.
74  !ym            sij(i,k,j)=0.0  
75   385      continue      ! CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
76   390    continue      ! RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
77   400  continue      ! FRACTION (sij)
78    
79  !ym      do i = minorig + 1, nl
80        ment(1:ncum,1:nd,1:nd)=0.0  
81        sij(1:ncum,1:nd,1:nd)=0.0         do j = minorig, nl
82              do il = 1, ncum
83  !      do k=1,ntra               if((i >= icb(il)).and.(i <= inb(il)).and. &
84  !       do j=1,nd  ! instead nlp                    (j >= (icb(il) - 1)).and.(j <= inb(il)))then
85  !        do i=1,nd ! instead nlp  
86  !         do il=1,ncum                  rti = rr(il, 1) - ep(il, i) * clw(il, i)
87  !            traent(il,i,j,k)=tra(il,j,k)                  bf2 = 1. + lv(il, j) * lv(il, j) * rs(il, j) / (rrv &
88  !         enddo                       * t(il, j) * t(il, j) * cpd)
89  !        enddo                  anum = h(il, j) - hp(il, i) + (cpv - cpd) * t(il, j) * (rti &
90  !       enddo                       - rr(il, j))
91  !      enddo                  denom = h(il, i) - hp(il, i) + (cpd - cpv) * (rr(il, i) &
92        zm(:,:)=0.                       - rti) * t(il, j)
93                    dei = denom
94  !=====================================================================                  if(abs(dei) < 0.01)dei = 0.01
95  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING                  sij(il, i, j) = anum / dei
96  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING                  sij(il, i, i) = 1.0
97  ! --- FRACTION (sij)                  altem = sij(il, i, j) * rr(il, i) + (1. - sij(il, i, j)) &
98  !=====================================================================                       * rti - rs(il, j)
99                    altem = altem / bf2
100        do 750 i=minorig+1, nl                  cwat = clw(il, j) * (1. - ep(il, j))
101                    stemp = sij(il, i, j)
102         do 710 j=minorig,nl                  if((stemp < 0.0.or.stemp > 1.0.or.altem > cwat) &
103          do 700 il=1,ncum                       .and.j > i)then
104           if( (i.ge.icb(il)).and.(i.le.inb(il)).and. &                     anum = anum - lv(il, j) * (rti - rs(il, j) - cwat * bf2)
105              (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then                     denom = denom + lv(il, j) * (rr(il, i) - rti)
106                       if(abs(denom) < 0.01)denom = 0.01
107            rti=rr(il,1)-ep(il,i)*clw(il,i)                     sij(il, i, j) = anum / denom
108            bf2=1.+lv(il,j)*lv(il,j)*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd)                     altem = sij(il, i, j) * rr(il, i) + (1. - sij(il, i, j)) &
109            anum=h(il,j)-hp(il,i)+(cpv-cpd)*t(il,j)*(rti-rr(il,j))                          * rti - rs(il, j)
110            denom=h(il,i)-hp(il,i)+(cpd-cpv)*(rr(il,i)-rti)*t(il,j)                     altem = altem - (bf2 - 1.) * cwat
111            dei=denom                  end if
112            if(abs(dei).lt.0.01)dei=0.01                  if(sij(il, i, j) > 0.0.and.sij(il, i, j) < 0.95)then
113            sij(il,i,j)=anum/dei                     qent(il, i, j) = sij(il, i, j) * rr(il, i) + (1. &
114            sij(il,i,i)=1.0                          - 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                          - sij(il, i, j)) * u(il, minorig)
117            cwat=clw(il,j)*(1.-ep(il,j))                     vent(il, i, j) = sij(il, i, j) * v(il, i) + (1. &
118            stemp=sij(il,i,j)                          - sij(il, i, j)) * v(il, minorig)
119            if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat) &                     elij(il, i, j) = altem
120                         .and.j.gt.i)then                     elij(il, i, j) = amax1(0.0, elij(il, i, j))
121             anum=anum-lv(il,j)*(rti-rs(il,j)-cwat*bf2)                     ment(il, i, j) = m(il, i) / (1. - sij(il, i, j))
122             denom=denom+lv(il,j)*(rr(il,i)-rti)                     nent(il, i) = nent(il, i) + 1
123             if(abs(denom).lt.0.01)denom=0.01                  end if
124             sij(il,i,j)=anum/denom                  sij(il, i, j) = amax1(0.0, sij(il, i, j))
125             altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j)                  sij(il, i, j) = amin1(1.0, sij(il, i, j))
126             altem=altem-(bf2-1.)*cwat               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, minorig) - ep(il, i) * clw(il, i)
137                 uent(il, i, i) = u(il, minorig)
138                 vent(il, i, i) = v(il, minorig)
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  
   
 !       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)  
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  !      do j=1,ntra               endif
274  !       do il=1,ncum            enddo
275  !        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)         end do
276  !     :     .and. csum(il,i).lt.m(il,i) ) then  
277  !         traent(il,i,i,j)=tra(il,nk(il),j)         do il = 1, ncum
278  !        endif            if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
279  !       enddo                 .and. csum(il, i) < m(il, i)) then
280  !      enddo               nent(il, i) = 0
281  789   continue               ment(il, i, i) = m(il, i)
282  !               qent(il, i, i) = rr(il, 1) - ep(il, i) * clw(il, i)
283  ! MAF: renormalisation de MENT               uent(il, i, i) = u(il, minorig)
284        do jm=1,nd               vent(il, i, i) = v(il, minorig)
285          do im=1,nd               elij(il, i, i) = clw(il, i)
286            do il=1,ncum               sij(il, i, i) = 0.0
           zm(il,im)=zm(il,im)+(1.-sij(il,im,jm))*ment(il,im,jm)  
          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)  
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.47  
changed lines
  Added in v.198
  ViewVC Help
Powered by ViewVC 1.1.21