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

Diff of /trunk/phylmd/CV30_routines/cv30_mixing.f90

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

Legend:
Removed from v.47  
changed lines
  Added in v.332

  ViewVC Help
Powered by ViewVC 1.1.21