         subroutine brume(ngrid,tab1,x1,
     &                    xnz,xnrad,taused,ihor,
     &                    pmu0,pfract,
     &                    precip)


*1  c nombre de particules de la grille de rayon r a l'altitude z 
*2  dt, pas de temps, en heure.  

*--------------------------------------------------------------*
*                                                              *
*            ENTRE 0 ET 1000 KILOMETRES                        *     
*                                                              * 
*    la dimension fractale est en tableau, attention au        *
*    raccordement entre le regime moleculaire et le regime     *
*    fluide                                                    *
*                                                              *
*    Modele microphysique:    Cabane et al.,1992 /             *
*    Modele version fractale: Cabane et al.,1993 /             *
*                                                              *
*--------------------------------------------------------------*
* VERSION DU 2 JUIN 1993  --- AUT 1994 --- 11/04/96
*
* changer: altitude de production z0=/taux de production ctot= 
*        : la charge/micron, ne
*        : df(h),rf... 
* raccordement aknc 
*
* declaration des blocs communs
*------------------------------

         use dimphy
         IMPLICIT NONE
#include "dimensions.h"
#include "microtab.h"
#include "varmuphy.h"


         real pmu0,pfract

         common/ctps/li,lf,dt
         common/con/c
         common/coag/k
         common/effets/xsaison

 
* declaration des variables communes
* ----------------------------------

         integer xnz,xnrad,ngrid
         integer li,lf,ihor
         real dt
         real c0(nz,nrad),c(nz,nrad,2)
         real k(nz,nrad,nrad),knu
         real xsaison
         real taused(nz,nrad)
         real precip(ngrid,5)

        

*  variables internes
*  ------------------

         integer h,ti,itime,i,j
         real tab1(nz,nrad)
         real x1 
         real somme,v1,dice3
         
         real vitesse

         save itime
         data itime/0/

* initialisation
* --------------

* effet saisonnier
* ----------------


         xsaison=0.
         xsaison=pmu0*4.*pfract
                              !=Pi si fract=1/2 (equinoxe) et 
                              !    si mu0(ihor)=1 sous le soleil
                              !    exactement.

c        xsaison=0.
c        if (ihor.le.9.or.ihor.ge.41) xsaison=8.  ! rapport des surfaces
c        xsaison=1.

* controles
* ---------

         if (nrad.ne.xnrad)  stop 'nrad.ne.xnrad'
         if (nz.ne.xnz)    stop 'nz.ne.xnz'

         do i=1,nz
           do j=1,nrad 
             c(i,j,1)=tab1(i,j)
             c(i,j,2)=0.0 
           enddo 
         enddo 

         dt=x1 

         if (itime.eq.0) then 
           ITIME=1
           call init
           call calcoag
         endif

         do i=1,nz,1
           do j=1,nrad
             v1=vitesse(i,j,0)
c             ho que c'est moche ! -> taused = RT/(Mn2*g)*1/vaer = H/v
             taused(i,j)=(8.314*t(i)/28.e-3/1.35)/v1
           enddo
         enddo

         call coagul

         call production(ihor)
 
         li=3-li
         lf=3-lf

         call sedif(dice3)
c        En theorie, dice3 est NEGATIF (en sedimentant on ne fait que perdre des aerosols)
c        Les precipitations sont comptees positivement. (ET ON NE PREND QUE DES VALEURS POSITIVES)
         precip(ihor,5)=AMAX1(-dice3/rhol,0.)    ! m3/m2=m 

         li=3-li
         lf=3-lf

         do i=1,nz
           do j=1,nrad
             tab1(i,j)=c(i,j,li)     ! li=1
           enddo
         enddo 
         
         return  

         end 
*________________________________________________________________________ 

         subroutine coagul

*********************************************************
* ce programme calcule la nouvelle concentration dans   *
* le a ieme intervalle de rayon, a l'altitude h, a      *
* l'instant t+dt                                        *
*********************************************************
         use dimphy
         IMPLICIT NONE
#include "dimensions.h"
#include "microtab.h"
#include "varmuphy.h"


* declaration des blocs communs
*------------------------------

         common/ctps/li,lf,dt
         common/con/c

* declaration des variables
* --------------------------

         integer li,lf
         real dt
         real c(nz,nrad,2)

* declaration des variables propres au ss-programme
* -------------------------------------------------

         integer h,a
         real pr,pe

*  traitement
*  ----------

         do h=nztop,nz
           do a=1,nrad
             call pertpro(h,a,pe,pr)
c            if((1+dt*pe).lt.0.) stop 'denom.eq.0'
             c(h,a,lf)=(c(h,a,li)+pr*dt)/(1+dt*pe)
           enddo
         enddo 

         if (nztop.ne.1) then
           do h=1,nztop-1
             do a=1,nrad
               c(h,a,lf)=c(h,a,li)
             enddo
           enddo 
         endif

         return
         end
      
      
*__________________________________________________________________________

         subroutine  calcoag

***************************************************************
*                                                             *
*  Ce programme calcule les coefficients de collection  d'une *
* particule de rayon x avec une particule de rayon b a une    *
* altitude donnee h                                           *
*************************************************************** 

* declaration des blocs communs
*------------------------------
         use dimphy
         IMPLICIT NONE
#include "dimensions.h"
#include "microtab.h"
#include "varmuphy.h"

         common/ctps/li,lf,dt
         common/con/c
         common/coag/k

* declaration des variables
* --------------------------

         integer li,lf,i
         real dt
         real c(nz,nrad,2)
         real knu,nud,k(nz,nrad,nrad)

* declaration des variables propres au ss-programme
* -------------------------------------------------

         integer h,b,x
         real nua,lambb,lambx,knb,knx,alphab,alphax,d,e,f,kcg
         real db,dx,rm,dm,deltab,deltax,del,g,beta,gx,gb
         real rfx,rfb,rpr
         real*8 ne,qe,epso
         real*8 corelec,yy

         real kco,vx,vb,vitesse,sto,ee,a,dd,bb,p0,t0,l0,ccol
         real st(37),ef(37)

* initialisation
* --------------
c        print*,'**** calcoag'

*  -nombres de STOCKES

         data(st(i),i=1,37)/1.35,1.5,1.65,1.85,2.05,2.25,2.5,2.8,3.1,
     s    3.35,3.6,3.95,4.3,4.7,5.05,5.45,5.9,6.4,7.,7.6,8.3,9.05,9.9,
     s       10.9,11.1,13.5,15.3,17.25,20.5,24.5,30.4,39.3,48,57,86.,
     s       187.,600./

*  -coef. d'efficacite de collection

         ef(1)=3.75
         ef(2)=8.75
         do i=3,37
           ef(i)=ef(i-1)+2.5
         enddo

         do i=1,37
           ef(i)=ef(i)*1e-2
         enddo

         qe=1.6e-19
         ne=-30.e+6         ! "vieille valeur!"
         ne=-15.e+6          ! pour fitter DISR !

         epso=1e-9/(36*pi)

         d=1.257
         e=0.4
         f=-1.1

*  iteration sur z
 
         do 1 h=1,nz
           nua=nud(h,1)      

*  iteration sur les rayons

           do 1 b=1,nrad

             knb=knu(h,b,1)
             vb=vitesse(h,b,1)

             do 1 x=1,b

               knx=knu(h,x,1)
               vx=vitesse(h,x,1)

**  COAGULATION  ****************************************************
** --------------****************************************************
* calcul du terme correcteur 'slip-flow'

               alphab=d+e*exp(f/knb)
               alphax=d+e*exp(f/knx)

* calcul du coefficient de diffusion

               rfb=(r_e(b)**(3./df(b)))*((rf(b))**(1.-3./df(b)))
               rfx=(r_e(x)**(3./df(x)))*((rf(x))**(1.-3./df(x)))
               db=kbz*t(h)*(1+alphab*knb)/(6*pi*nua*rfb)
               dx=kbz*t(h)*(1+alphax*knx)/(6*pi*nua*rfx)

* calcul du coefficient de coagulation

               rpr=rfb+rfx
               kcg=4*pi*rpr*(db+dx)

* calcul de la vitesse thermique

               gx=sqrt(6*kbz*t(h)/(rhol*pi**2*r_e(x)**3))
               gb=sqrt(6*kbz*t(h)/(rhol*pi**2*r_e(b)**3))

* calcul du libre parcours apparent des aerosols

               lambb=8*db/(pi*gb)
               lambx=8*dx/(pi*gx)

*calcul du terme correcteur beta

               rm=rpr/2.
               dm=(dx+db)/2.
               g=sqrt(gx**2+gb**2)
               deltab=(((2*rfb+lambb)**3-(4*rfb**2+lambb**2)**1.5)
     s         /(6*rfb*lambb)-2*rfb)*sqrt(2.)
               deltax=(((2*rfx+lambx)**3-(4*rfx**2+lambx**2)**1.5)
     s         /(6*rfx*lambx)-2*rfx)*sqrt(2.)
               del=sqrt(deltab**2+deltax**2)
               beta=1/((rm/(rm+del/2))+(4*dm/(g*rm)))

* calcul du coefficient de coagulation corrige

               kcg=kcg*beta
c       print*,' kcg:', rfb,rfx,knb,knx,db,dx
c       print*,' beta:', gx,gb,lambb,lambx,rm,dm
c       print*,' beta:', g,deltab,deltax,del,beta
c       print*,' '

**  COALESCENCE  **************************************************
**  -------------**************************************************

               kco=0.

               if (b.eq.x) goto 9

* calcul du nombre de Stockes de la petite particule

               sto=2*rhol*rfx**2*abs(vx-vb)/(9*nua*rfb)

* calcul du coef. de Cunningham-Millikan

               a=1.246
               bb=0.42
               dd=0.87
               l0=0.653e-7
               p0=101325.
               t0=288.

               ee=1+(l0*t(h)*p0*
     &         (a+bb*exp(-dd*rfx*t0*p(h)/(l0*t(h)*p0))))
     s         /(rfx*t0*p(h))

* calcul du nombre de Stockes corrige

               sto=sto*ee

               if (sto .le. 1.2) goto 9

               if (sto .ge. 600.) then 
                 ccol=1.
                 goto 8
               endif

*  recherche du coefficient de collection

               do 3 i=1,37
                 if (sto .gt. st(i)) then
                   goto 3
                 endif
                 if (sto .eq. st(i)) then
                   ccol=ef(i+1)
                 else
                   ccol=ef(i)
                 endif
                 goto 8
3              continue

*  calcul du coefficient de coalescence

8              kco=pi*(rfb+rfx)**2*ccol*abs(vb-vx)

9              continue

**  CORRECTION ELECTRICITE *******************************
**  ------------------------******************************

               yy=1.d0*ne**2*r_e(x)*r_e(b)*qe**2
     &         /(1.d0*kbz*t(h)*(r_e(b)+r_e(x))*4*pi*epso)


               corelec=0.
               if (yy.lt.50.) corelec=yy/(exp(yy)-1.)
               if (yy.le.1.e-3)  corelec=1.

               k(h,b,x)=(kcg+kco)*corelec
               k(h,x,b)=k(h,b,x)
 

1        continue
         return
         end

*______________________________________________________________________

         real function lambda(j,indic)
*
*------------------------------------------------------------------*
*  fonction calculant le libre parcours moyen des molecules        *
*  atmospheriques( rayon =ra) se trouvant dans la couche no j.     *
*  pour indic=0  ...... la particule se trouve a la frontiere entre*
*                        les couches j et j-1                      *
*  pour indic=1  ...... la particule se trouve au milieu de la     *
*                         la couche j                              *
*------------------------------------------------------------------*
*
* declaration des blocs communs
*------------------------------
         IMPLICIT NONE
#include "dimensions.h"
#include "microtab.h"
#include "varmuphy.h"

* declaration des variables communes
* ----------------------------------

         integer i,j

*  declaration des variables internes
*  ----------------------------------

         integer indic
         real pp,ra

         ra=1.75e-10

* traitement
* ----------

         if (indic.eq.0) then
           pp=pb(j)
         else
           if (indic.ne.1) then
             print*,'erreur argument fonction lambda'
             return
           endif
           pp=p(j)
         endif

         lambda=kbz*t(j)/(4*sqrt(2.)*pi*(ra**2)*pp)
         end

*******************************************************************************

         real function knu(j,k,indic)
*
*--------------------------------------------------------------*
*  fonction calculant le nombre de knudsen d'une particule     *
*  d'aerosol de rayon r_e(k) se trouvant dans la couche no j   *
*  indic ......  idem function lambda                          *
*--------------------------------------------------------------*
*
* declaration des blocs communs
*------------------------------
         IMPLICIT NONE
#include "dimensions.h"
#include "microtab.h"
#include "varmuphy.h"

*  declaration des variables internes
*  ----------------------------------

         integer indic,j,k
         real lambda,rfk

* traitement
* ----------

         if (indic.ne.0 .and.indic.ne.1) then
           print*,'erreur argument fonction knu'
           return
         endif

         rfk=(r_e(k)**(3./df(k)))*((rf(k))**(1.-3./df(k)))
         knu=lambda(j,indic)/rfk
         end

*****************************************************************************

         real function nud(j,indic)
*
*--------------------------------------------------------------*
*  fonction calculant la viscosite dynamique (en USI) de l'air *
*  d'apres la formule de Sutherlant a l'altitude j             *
*  indic  ......... idem fonction lambda                       *
*--------------------------------------------------------------*
*
         IMPLICIT NONE
#include "dimensions.h"
#include "microtab.h"
#include "varmuphy.h"


         integer indic,j
         real nud0,c,tt
*
         nud0=1.74e-5
         c=109.
 
         if(indic.ne.0.and.indic.ne.1) then
           print*,'erreur argument fonction nud'
           return
         endif
 
         if(indic.eq.0) tt=tb(j)
         if (indic.eq.1) tt=t(j)
 
         nud=nud0*sqrt(tt/293)*(1+c/293)/(1+c/tt)
         end

****************************************************************************

         real function vitesse(j,k,indic)
*
*-----------------------------------------------------------------*
*   fonction calculant la vitesse de chute d'une particule de rayon*
*   k se trouvant a l'altitude j  suivant la valeur du nombre de   *
*    Knudsen                                                       *
*   indic ....... idem function lambda                             *
*-----------------------------------------------------------------*
*
*  declaration des blocs communs
*------------------------------
         IMPLICIT NONE
#include  "dimensions.h"
#include  "microtab.h"
#include  "varmuphy.h"

*   declaration des variables internes
*   ----------------------------------

         integer indic,j,k
         real w,g,m,a0,zz,nud,knud,tt,rhoh
         real rbis, rfk,vb,zbx
         real akncx,knu

*   traitement
*   ----------

         if (indic.ne.0.and.indic.ne.1) then
           print*,'erreur argument fonction vitesse'
           return
         endif

         if(indic.eq.0) then
           zz=z(j)+dz(j)/2.
           tt=tb(j)
           rhoh=rhob(j)
         endif
         if(indic.eq.1) then
           zz=z(j)
           tt=t(j)
           rhoh=rho(j)
         endif

         g=g0*(rtit/(rtit+zz))**2
         a0=0.74
         m=(ach4(j)*mch4+aar(j)*mar+an2(j)*mn2)/nav
         knud=knu(j,k,indic)


         rfk=(r_e(k)**(3./df(k)))*((rf(k))**(1.-3./df(k)))
c        rfk=r_e(7)


          w=2./9.*rfk**(df(k)-1.)*rf(k)**(3.-df(k))*g*rhol/nud(j,indic)

         w=w*(1+1.2517*knud+0.4*knud*exp(-1.1/knud))


c        if (p(j).lt.500..and.k.eq.nrad) then
c           w=0.
c        endif

         vitesse=w

         end
***********************************************************************

         real function kd(h)
*
*--------------------------------------------------------------------*
*   cette fonction calcule le coefficient du terme de 'eddy diffusion'*
*   a l altitude j						     *
*--------------------------------------------------------------------*
*
         IMPLICIT NONE
#include  "dimensions.h"
#include  "microtab.h"
#include  "varmuphy.h"

         real zbx

         integer h

         zbx=z(h)+dz(h)/2.
         if(zbx.le.42000.) then
           kd=1.64e+12*(pb(h)/(kbz*tb(h)))**(-1./2.)
           kd=4.
         else
           kd=0.0*kd
           kd=1.64e+12*(pb(h)/(kbz*tb(h)))**(-1./2.)
         endif

         if(zbx.le.50000.) then
           kd=1.64e+12*(pb(h)/(kbz*tb(h)))**(-1./2.)
         endif

         kd=0.0*kd

         return
         end


*____________________________________________________________________________

         subroutine init
*
*--------------------------------------------------------------------*
*   cette routine effectue  :                                         *
*                1) interpolation a partir des donnees initiales des  *
*                    valeurs de p,t,rho,ach4,aar,an2  sur la grille   *
*                2) initialisation des constantes (common/phys/)      *
*                3) initialisation des variables temporelles (common  *
*                     /temps/)                                        *
*                4) definition des grilles en rayon et verticale      *
*                5)  initialisation de c(z,r,t) avec les donnees du   *
*                      fichier unit=1                                 *
*                                                                     *
*   les donnees sont des valeurs caracterisques de l'atmosphere de    *
*     TITAN  ( voir Lelouch and co )                                  *
*--------------------------------------------------------------------*

*  declaration des blocs communs
*------------------------------
         use dimphy
         IMPLICIT NONE
#include  "dimensions.h"
#include  "microtab.h"
#include  "varmuphy.h"


         common/ctps/li,lf,dt
         common/con/c

*  declaration des variables communes
*  ----------------------------------

         real c(nz,nrad,2)
         integer li,lf
         integer i,ii
         real dt

*  declaration des variables internes
*  ----------------------------------

         integer nzd
         parameter (nzd=254)
         integer limsup,liminf,j1,j2
         real zd(nzd),ach4d(nzd),rap
         real m


*  initialisation des variables temporelles
*  ----------------------------------------

         li=1
         lf=2


*  interpolation de xch4,xar et xn2 sur la grille
*  ----------------------------------------------

*  donnees initiales (Lellouch et al,87) 
*  ------------------------------------- 

c        print*,'****** init'
         do 1 i=1,168
           zd(i)=(1000.-5*(i-1))*1000.
1        continue
         do 2 i=1,78
           zd(168+i)=(160.-2*(i-1))*1000.
2        continue
         do 3 i=1,4
           zd(246+i)=(5.-(i-1))*1000.
3        continue
         do 4 i=1,4
           zd(250+i)=(1.5-(i-1)*0.5)*1000.
4        continue

         data (ach4d(i),i=1,168)/168*1.5e-2/
         data (ach4d(i),i=169,254)/63*1.5e-2,1.6e-2,1.8e-2,1.8e-2,
     &   1.9e-2,2.e-2,2.1e-2,2.3e-2,2.5e-2,2.8e-2,3.1e-2,3.6e-2,
     &   4.1e-2,4.7e-2,5.7e-2,6.7e-2,7.5e-2,7*8.e-2/

         liminf=0
         limsup=0

*  interpolation des taux de melange de ch4,ar,n2  
*-----------------------------------------------   

!        do 20 j1=1,nz
!           do 21 j2=1,nzd
!               if( zd(j2).le.z(j1)) goto 22
!21	    continue
!22         liminf=j2

          do 20 j1=1,nz
             do 21 j2=1,nzd
                 if( zd(j2).le.z(j1)) goto 22
  21         continue
  22         if (j2.ge.254) j2=254
             liminf = j2

          if (zd(liminf).eq.z(j1) )then
            ach4(j1)=ach4d(liminf)
            goto 20
          endif
          if (j2.ne.1) then
            limsup=j2-1
          else
            limsup=j2
          endif

          if (limsup.eq.liminf) then
            ach4(j1)=ach4(limsup)
          else
            ach4(j1)=ach4d(liminf)-(ach4d(limsup)-ach4d(liminf))/
     s       (zd(limsup)-zd(liminf))*(zd(liminf)-z(j1))
          endif
20      continue

*   rap= aar/an2  cst sur l'altitude

         rap=0.191
         do 23 i=1,nz
           an2(i)=(1.-ach4(i))/(1.+rap)
           aar(i)=rap*an2(i)
23       continue

         do 24 i=1,nz
           m=ach4(i)*mch4+an2(i)*mn2+aar(i)*mar
           rho(i)=p(i)*m/(rgp*t(i))
24       continue

         do 34 i=1,nz
           m=ach4(i)*mch4+an2(i)*mn2+aar(i)*mar
           rhob(i)=pb(i)*m/(rgp*tb(i))
c          print*,pb(i),m,rgp,tb(i),rhob(i),rho(i)
34       continue

*  fin d'interpolation des taux de melange
*----------------------------------------  

c        print*,'**** fin init'
540      continue
         return

500       print*,'erreur lecture initialisation de c...erreur=',ii
          stop

        end

*____________________________________________________________________________

         subroutine pertpro(h,a,l_,pr_)

*****************************************************************************
*                                                                           *
*  ce programme permet le calcul du terme de production (pr) et de perte (l)*
*  pour le phenomene de coagulation                                         *
*  dans le a ieme intervalle de rayon a une altitude h                      *
*****************************************************************************

*  declaration des blocs communs
*------------------------------
         use dimphy
         IMPLICIT NONE
#include  "dimensions.h"
#include  "microtab.h"
#include  "varmuphy.h"


         common/ctps/li,lf,dt
         common/con/c
         common/coag/k

*  declaration des variables
*  --------------------------

         integer li,lf
         real dt
         real c(nz,nrad,2),k(nz,nrad,nrad)

*  declaration des variables propres au ss-programme
*  -------------------------------------------------

         integer h,b,a,x,i
         real*8 pr,ss,s,l
         real pr_,l_,vol,del

*  traitement
*  -----------

*   production
*+++++++++++++
         s=0.d0
         ss=0.d0
         pr=0.

         if (a .eq. 1) goto 2  
         b=a-1

         if (c(h,b,lf) .eq. 0 .and. c(h,b,li) .eq. 0) goto 2
  
         do 1 i=1,b

         if(c(h,i,li) .eq. 0 .and. c(h,i,lf) .eq. 0) goto 1
     
         if (i .ne. b)del=1.
         if (i .eq. b) del=.5

         s=(v_e(i)*1.d0)*del*(k(h,b,i)*1.d0)*(c(h,i,li)*1.d0)+s
         ss=(v_e(i)*1.d0)*del*(k(h,b,i)*1.d0)*(c(h,i,lf)*1.d0)+ss

c        if (a.eq.2) print*,'SS>',v_e(i),k(h,b,i),c(h,b,lf)
c        if (a.eq.2) print*,'SS>',del*v_e(i)*k(h,b,i)*c(h,b,lf)


1        continue

*   calcul du terme de production

          pr=(c(h,b,lf)*s/(vrat_e-1.)+c(h,b,li)*ss)/v_e(a)
c         if (a.eq.2) print*,'PR>',s,ss,c(h,b,lf),v_e(a)

2         continue


*   perte
*-  - - - -

          l=0


*     condition limite : pas de perte dans le dernier intervalle

         if (a .eq. nrad) goto 9

         do 10 x=1,nrad

         if (c(h,x,li) .eq. 0) goto 10

         if (a .lt. x) vol=1.
         if (a .eq. x) vol=.5*vrat_e/(vrat_e-1)
         if (a .gt. x) vol=v_e(x)/(v_e(a)*(vrat_e-1))

         l=l+k(h,a,x)*c(h,x,li)*vol*1.d0

      
10       continue
9        continue

#ifdef  CRAY
         l_=l
         pr_=pr
#else
         l_=sngl(l)
         pr_=sngl(pr)
#endif
c        l_=sngl(l)
c        pr_=sngl(pr)
c          if (a.eq.2) print*,'pr_,l_',h,a,pr_,l_
c          if (a.eq.2) print*,'-----------------------'
c          if (a.eq.2) STOP



         return

         end

*_____________________________________________________________________________

         subroutine production(ihor)
*
*--------------------------------------------------------------------*
*   routine calculant le terme de production des molecules organiques *
*   composant les aerosols . rini= rayon des aerosols initiaux        *
*--------------------------------------------------------------------*
*
         use dimphy
         IMPLICIT NONE
#include  "dimensions.h"
#include  "microtab.h"
#include  "varmuphy.h"
#include  "clesphys.h"

c#include  "aerprod.h"


         integer ndz,i,ihor,k,k1
         real c(nz,nrad,2)
         integer li,lf
         real dt
         real zprod,zy,c0,ctot,prod,rini,rfron
         real xsaison,p0
         common/ctps/li,lf,dt
         common /con/c
         common/effets/xsaison

! Pressure level of aerosol production
          p0=p_prodaer

          do i=1,nz-1
           if (pb(i).lt.p0.and.pb(i+1).gt.p0) zprod=(z(i)+z(i+1))/2.
          enddo

         ctot=3.5e-13*tx ! ATTENTION, ??COHERENT AVEC INITPAR??
         ctot=ctot*xsaison  ! 

c        z0=385.e+3
         zy=20.e+3
         rini=1.3e-9
         rini=r_e(1)   
         ndz=50
*
         do 10 i=1,nrad 
           if(rini.lt.r_e(i)) goto 100
10       continue
100      continue
         if (i.eq.1) then
           rini=r_e(1)
         else
           rfron=(r_e(i)+r_e(i-1))/2
           if (rini .lt.rfron) then 
             rini=r_e(i-1)
             i=i-1
           else
             rini=r_e(i)
           endif
         endif
*
         c0=ctot/(sqrt(2.*pi)*zy)
         c0=c0*3./(4.*pi*rhol*rini**3)
*
         do 20 k=nztop,nz
         prod=0.
         do 201 k1=1,ndz
           prod=prod+c0*exp(-0.5*(((z(k)+dz(k)/2.-k1*dz(k)/(2.*ndz)
     s     -zprod)/zy)**2))*dt/ndz
201      continue

          if (prod .le. 1) prod=0.
            c(k,i,lf)=c(k,i,lf)+prod
20       continue

!        do 30 k=nztop,nz
!        c(k,i,lf)=c(k,i,lf)+prodaer(ihor,nz-k+1,1)
!     & *3./(4.*pi*rhol*rini**3)
!30      continue

         return
         end

*-------------------------------------------------------------------*



        subroutine sedif(dice3)
*
*------------------------------------------------------------------*
*   cette routine calcule l'evolution de la fonction de distribution*
*   c(z,r,t) pour les phenomenes de sedimentation et de diffusion   *
*------------------------------------------------------------------*
*
*
*  declaration des blocs communs
*------------------------------
         use dimphy
         IMPLICIT NONE
#include  "dimensions.h"
#include  "microtab.h"
#include  "varmuphy.h"

         common/ctps/li,lf,dt
         common/con/c

*  declaration des variables communes
*  ----------------------------------

         integer li,lf
         integer i,j,k,nb
         real dt
         real c(nz,nrad,2),dice3,bilan4,bilan14

*  declaration des variables internes
*  ----------------------------------

         real w,w1,dzbX,dc
         double precision sigma,theta,hc,l,rap,cmp,wp 
         double precision fs(nz+1),ft(nz+1)
         real as(nz),bs(nz),cs(nz),ds(nz)
         double precision asi(nztop:nz),bsi(nztop:nz),
     &                    csi(nztop:nz)
         double precision dsi(nztop:nz),xsol(nztop:nz)
         real vitesse,kd

         external dtridgl
*  resolution
*------------

            
            bilan4=0.
            do  k=1,nrad
            do  j=nztop,nz
              bilan4=bilan4+c(j,k,li)*dzb(j)*
     &        4./3.*pi*rf(k)**3.*vrat_e**(k-imono)
            enddo
            enddo


            do 10 k=1,nrad  
            do 20 j=nztop,nz

            if (j.eq.1) goto 20

*  calcul de la vitesse corrigee

              dzbX=(dz(j)+dz(j-1))/2.
              w=-1*vitesse(j,k,0)
              if (kd(j).ne.0.) then
                theta=0.5*(w*dzbX/kd(j)+log(rho(j-1)/rho(j)))
                if (theta.ne.0) then
                  sigma=1./dtanh(theta)-1./theta
                else
                  sigma=1.
                endif
              else
                sigma=1.
              endif
              if(c(j,k,li).eq.0.) then
                rap=10.
              else
                rap=c(j-1,k,li)/c(j,k,li)
                if( rap.gt.10.) rap=10.
                if( rap.lt.0.1) rap=0.1
              endif
              if (rap.gt.0.9 .and. rap.lt.1.1) then
                w1=w
              else
                if(w.ne.0) then
                  hc=dzbX/dlog(rap)
                  l=dzbX/(w*dt)*(dexp(-w*dt/hc)-1.)/(1.-rap)
                  wp=w*1.d0
                  cmp=dlog(-wp)+abs(sigma)*dlog(l)
                  if (cmp.gt.38) then
                    goto 20
                  endif
                  w1=-dexp(cmp)

                else
                  w1=0.
                endif
              endif

c               if(w1.ge.0. .or. w1.le.0.) then 
c                continue
c               else
c                print*,j,k,w1,hc,dzbx,rap,dlog(rap),' w1'
c               endif

*   calcul des flux aux interfaces


             if (kd(j).ne.0.) then
               if (theta.ne.0.) then
                 ft(j)=(w1+log(rho(j-1)/rho(j))*kd(j)/dzbX)/(dexp(2.*
     s           theta)-1.)
                 fs(j)=ft(j)*dexp(2.*theta)
               else
                 ft(j)=kd(j)/dzbX
                 fs(j)=kd(j)/dzbX
               endif
            else
              if (w1.lt.0.)then
                ft(j)=-w1
                fs(j)=0.
              else
                ft(j)=0.
                fs(j)=w1
              endif
            endif

20          continue

*  conditions aux limites pour les flux aux interfaces

            fs(1)=0.
            ft(1)=0.
            fs(nz+1)=0.
            ft(nz+1)=-w1

*  calcul des coefficients de l'equation discrete

            do 21 j=nztop,nz
              as(j)=-dz(j)/dt
              bs(j)=-ft(j)
              cs(j)=ft(j+1)+fs(j)-dz(j)/dt
              ds(j)=-fs(j+1)

              if ( cs(j).gt.0) goto 100
21         continue

*  cas explicite (mu=0) : calcul de la fonction c(z,r,t+1)

           do 22 j=nztop,nz-1

           if (j.eq.nztop) then
             dc=(cs(nztop)*c(nztop,k,li)+ds(nztop)
     &                    *c(nztop+1,k,li))/as(nztop)
             c(nztop,k,lf)=dc
             goto 22
           endif

           dc=(bs(j)*c(j-1,k,li)+cs(j)*c(j,k,li)+ds(j)*c(j+1,k,li))
     s         /as(j)
           c(j,k,lf)=dc


22          continue

            dc=(bs(nz)*c(nz-1,k,li)+cs(nz)*c(nz,k,li))/as(nz)
            c(nz,k,lf)=dc

           if (nztop.ne.1) then
           do 32 j=1,nztop-1
             c(j,k,lf)=c(j,k,li)
32         continue
           endif

           goto 10

100        continue

*  cas implicite (mu=1) : calcul de la fonction c(z,r,t+1)

          do 101 j=nztop,nz
            asi(j)=ft(j)
            bsi(j)=-(ft(j+1)+fs(j)+dz(j)/dt)
            csi(j)=fs(j+1)
            dsi(j)=-dz(j)/dt*c(j,k,li)
            xsol(j)=0.
101       continue

*  inversion de la matrice tridiagonale 

         nb=nz-nztop+1

         call dtridgl(nb,asi,bsi,csi,dsi,xsol) 


         do 102 j=nztop,nz
          c(j,k,lf)=xsol(j)
102      continue

         if (nztop.ne.1) then
           do 110 j=1,nztop-1
           c(j,k,lf)=c(j,k,li)
110        continue
         endif



10       continue

         bilan14=0.
         do  k=1,nrad
         do  j=nztop,nz
           bilan14=bilan14+c(j,k,lf)*dzb(j)*
     &     4./3.*pi*rf(k)**3.*vrat_e**(k-imono)
         enddo
         enddo


         dice3=(bilan14-bilan4)*rhol


         return

          end

