subroutine stokes(p,t,rd,w,rho_aer)

!==================================================================
!     Purpose
!     -------
!     Compute the sedimentation velocity in a low pressure 
!     atmosphere.
!
!     Authors
!     ------- 
!     Francois Forget (1997)
!
!==================================================================

!      use radcommon_h, only : Rgas

      implicit none

#include "comcstfi.h"

!     input
!     -----
!     pressure (Pa), Temperature (K), particle radius (m), density
      real p, t, rd, rho_aer 

!     output
!     ------
!     sedimentation velocity (m/s, >0)
      real w

!     locally saved variables
!     -----------------------
      real a,b
      save a,b
  
      LOGICAL firstcall
      SAVE firstcall
      DATA firstcall/.true./

      if (firstcall) then

         !print*,'Routine not working: replace Rgas with r'
	 !stop
         !a = 0.707*Rgas/(4*pi*molrad**2 * avocado)
         a = 0.707*R/(4*pi*molrad**2 * avocado)
         b = (2./9.) * rho_aer * g / visc  
         firstcall=.false.
      end if

!     Sedimentation velocity =
!     Stokes' Law corrected for low pressures by the Cunningham
!     slip-flow correction according to Rossow (Icarus 36, 1-50, 1978)
      w = b * rd*rd * (1 + 1.333* (a*T/P)/rd )  

      return
    end subroutine stokes