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

Diff of /trunk/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 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 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 dimphy, ONLY: klev, klon
18  ! inputs:      use suphec_m, only: rcpd, rcpv, rv
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, 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  ! ori        do 360 i=1,ncum*nlp      do j = 1, nl
56          do 361 j=1,nl         do i = 1, ncum
57          do 360 i=1,ncum            nent(i, j) = 0
58            nent(i,j)=0         end do
59  ! in convect3, m is computed in cv3_closure      end do
60  ! ori          m(i,1)=0.0  
61   360    continue      do j = 1, nl
62   361    continue         do k = 1, nl
63              do i = 1, ncum
64  ! ori      do 400 k=1,nlp               qent(i, k, j) = rr(i, j)
65  ! ori       do 390 j=1,nlp               uent(i, k, j) = u(i, j)
66        do 400 j=1,nl               vent(i, k, j) = v(i, j)
67         do 390 k=1,nl               elij(i, k, j) = 0.0
68            do 385 i=1,ncum            end do
69              qent(i,k,j)=rr(i,j)         end do
70              uent(i,k,j)=u(i,j)      end do
71              vent(i,k,j)=v(i,j)  
72              elij(i,k,j)=0.0      ment(1:ncum, 1:klev, 1:klev) = 0.0
73  !ym            ment(i,k,j)=0.0      sij(1:ncum, 1:klev, 1:klev) = 0.0
74  !ym            sij(i,k,j)=0.0  
75   385      continue      zm(:, :) = 0.
76   390    continue  
77   400  continue      ! CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
78        ! RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
79  !ym      ! FRACTION (sij)
80        ment(1:ncum,1:nd,1:nd)=0.0  
81        sij(1:ncum,1:nd,1:nd)=0.0      do i = minorig + 1, nl
82    
83  !      do k=1,ntra         do j = minorig, nl
84  !       do j=1,nd  ! instead nlp            do il = 1, ncum
85  !        do i=1,nd ! instead nlp               if((i >= icb(il)).and.(i <= inb(il)).and. &
86  !         do il=1,ncum                    (j >= (icb(il) - 1)).and.(j <= inb(il)))then
87  !            traent(il,i,j,k)=tra(il,j,k)  
88  !         enddo                  rti = rr(il, 1) - ep(il, i) * clw(il, i)
89  !        enddo                  bf2 = 1. + lv(il, j) * lv(il, j) * rs(il, j) / (rv &
90  !       enddo                       * t(il, j) * t(il, j) * rcpd)
91  !      enddo                  anum = h(il, j) - hp(il, i) + (rcpv - rcpd) * t(il, j) * (rti &
92        zm(:,:)=0.                       - rr(il, j))
93                    denom = h(il, i) - hp(il, i) + (rcpd - rcpv) * (rr(il, i) &
94  !=====================================================================                       - rti) * t(il, j)
95  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING                  dei = denom
96  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING                  if(abs(dei) < 0.01)dei = 0.01
97  ! --- FRACTION (sij)                  sij(il, i, j) = anum / dei
98  !=====================================================================                  sij(il, i, i) = 1.0
99                    altem = sij(il, i, j) * rr(il, i) + (1. - sij(il, i, j)) &
100        do 750 i=minorig+1, nl                       * rti - rs(il, j)
101                    altem = altem / bf2
102         do 710 j=minorig,nl                  cwat = clw(il, j) * (1. - ep(il, j))
103          do 700 il=1,ncum                  stemp = sij(il, i, j)
104           if( (i.ge.icb(il)).and.(i.le.inb(il)).and. &                  if((stemp < 0.0.or.stemp > 1.0.or.altem > cwat) &
105              (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then                       .and.j > i)then
106                       anum = anum - lv(il, j) * (rti - rs(il, j) - cwat * bf2)
107            rti=rr(il,1)-ep(il,i)*clw(il,i)                     denom = denom + lv(il, j) * (rr(il, i) - rti)
108            bf2=1.+lv(il,j)*lv(il,j)*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd)                     if(abs(denom) < 0.01)denom = 0.01
109            anum=h(il,j)-hp(il,i)+(cpv-cpd)*t(il,j)*(rti-rr(il,j))                     sij(il, i, j) = anum / denom
110            denom=h(il,i)-hp(il,i)+(cpd-cpv)*(rr(il,i)-rti)*t(il,j)                     altem = sij(il, i, j) * rr(il, i) + (1. - sij(il, i, j)) &
111            dei=denom                          * rti - rs(il, j)
112            if(abs(dei).lt.0.01)dei=0.01                     altem = altem - (bf2 - 1.) * cwat
113            sij(il,i,j)=anum/dei                  end if
114            sij(il,i,i)=1.0                  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                          - sij(il, i, j)) * rti
117            cwat=clw(il,j)*(1.-ep(il,j))                     uent(il, i, j) = sij(il, i, j) * u(il, i) + (1. &
118            stemp=sij(il,i,j)                          - sij(il, i, j)) * u(il, minorig)
119            if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat) &                     vent(il, i, j) = sij(il, i, j) * v(il, i) + (1. &
120                         .and.j.gt.i)then                          - sij(il, i, j)) * v(il, minorig)
121             anum=anum-lv(il,j)*(rti-rs(il,j)-cwat*bf2)                     elij(il, i, j) = altem
122             denom=denom+lv(il,j)*(rr(il,i)-rti)                     elij(il, i, j) = amax1(0.0, elij(il, i, j))
123             if(abs(denom).lt.0.01)denom=0.01                     ment(il, i, j) = m(il, i) / (1. - sij(il, i, j))
124             sij(il,i,j)=anum/denom                     nent(il, i) = nent(il, i) + 1
125             altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j)                  end if
126             altem=altem-(bf2-1.)*cwat                  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  
   
 !       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)  
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  !      do j=1,ntra               endif
276  !       do il=1,ncum            enddo
277  !        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)         end do
278  !     :     .and. csum(il,i).lt.m(il,i) ) then  
279  !         traent(il,i,i,j)=tra(il,nk(il),j)         do il = 1, ncum
280  !        endif            if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
281  !       enddo                 .and. csum(il, i) < m(il, i)) then
282  !      enddo               nent(il, i) = 0
283  789   continue               ment(il, i, i) = m(il, i)
284  !               qent(il, i, i) = rr(il, 1) - ep(il, i) * clw(il, i)
285  ! MAF: renormalisation de MENT               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
           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)  
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.47  
changed lines
  Added in v.201
  ViewVC Help
Powered by ViewVC 1.1.21