/[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

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