source: LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/calc_Re.F90 @ 5185

  subroutine calc_Re(Q,Np,rho_a, &
             dtype,dmin,dmax,apm,bpm,rho_c,p1,p2,p3, &
             Re)
  use math_lib 
  implicit none

! Purpose:
!   Calculates Effective Radius (1/2 distribution 3rd moment / 2nd moment). 

!   For some distribution types, the total number concentration (per kg), Np
!   may be optionally specified.   Should be set to zero, otherwise.

!   Roj Marchand July 2010


! Inputs:

!   [Q]        hydrometeor mixing ratio (g/kg)  ! not needed for some distribution types
!   [Np]       Optional Total number concentration (per kg).  0 = use defaults (p1, p2, p3)

!   [rho_a]    ambient air density (kg m^-3)   

!   Distribution parameters as per quickbeam documentation.
!   [dtype]    distribution type
!   [dmin]     minimum size cutoff (um)
!   [dmax]     maximum size cutoff (um)
!   [apm]      a parameter for mass (kg m^[-bpm])
!   [bmp]      b params for mass 
!   [p1],[p2],[p3]  distribution parameters


! Outputs:
!   [Re]       Effective radius, 1/2 the 3rd moment/2nd moment (um)

! Created:
!   July 2010  Roj Marchand

! ----- INPUTS -----  
  
  real*8, intent(in) :: Q,Np,rho_a
 
  integer, intent(in):: dtype
  real*8, intent(in) :: dmin,dmax,rho_c,p1,p2,p3
    
  real*8, intent(inout) :: apm,bpm  
    
! ----- OUTPUTS -----

  real*8, intent(out) :: Re
  
! ----- INTERNAL -----
 
  integer :: local_dtype
  real*8  :: local_p3,local_Np

  real*8 :: pi, &
  N0,D0,vu,dm,ld, &         ! gamma, exponential variables
  rg,log_sigma_g
  
  real*8 :: tmp1,tmp2

  pi = acos(-1.0)

  ! // if density is constant, set equivalent values for apm and bpm
  if ((rho_c > 0) .AND. (apm < 0)) then
    apm = (pi/6)*rho_c
    bpm = 3.
  endif

  ! Exponential is same as modified gamma with vu =1
  ! if Np is specified then we will just treat as modified gamma
  if(dtype.eq.2 .AND. Np>0) then
    local_dtype=1;
    local_p3=1;
  else
    local_dtype=dtype;
    local_p3=p3;
  endif
  
  select case(local_dtype)
  
! ---------------------------------------------------------!
! // modified gamma                                        !
! ---------------------------------------------------------!
! :: Np = total number concentration (1/kg) = Nt / rho_a
! :: D0 = characteristic diameter (um)
! :: dm = mean diameter (um) - first moment over zeroth moment
! :: vu = distribution width parameter 

  case(1)  
  
    if( abs(local_p3+2) < 1E-8) then
  
    if(Np>1E-30) then
        ! Morrison scheme with Martin 1994 shape parameter (NOTE: vu = pc +1)
        ! fixed Roj. Dec. 2010 -- after comment by S. Mcfarlane
        vu = (1/(0.2714 + 0.00057145*Np*rho_a*1E-6))**2 ! units of Nt = Np*rhoa = #/cm^3
    else
        PRINT *, 'Error: Must specify a value for Np in each volume', &
             ' with Morrison/Martin Scheme.'
            stop    
    endif
    
    elseif (abs(local_p3+1) > 1E-8) then

      ! vu is fixed in hp structure  
      vu = local_p3 

    else

      ! vu isn't specified
      
      PRINT *, 'Error: Must specify a value for vu for Modified Gamma distribution'
      stop    
      
    endif
    

    if( Np.eq.0 .AND. p2+1 > 1E-8) then     ! use default value for MEAN diameter as first default
      
        dm = p2             ! by definition, should have units of microns
    D0 = gamma(vu)/gamma(vu+1)*dm
        
    else   ! use value of Np
        
        if(Np.eq.0) then
        
            if( abs(p1+1) > 1E-8 ) then  !   use default number concentration 
            
                local_Np = p1 ! total number concentration / pa --- units kg^-1
            else
            PRINT *, 'Error: Must specify Np or default value ', &
                 '(p1=Dm [um] or p2=Np [1/kg]) for ', &
                 'Modified Gamma distribution'
                stop
            endif
        else
            local_Np=Np;    
    endif
    
    D0 = 1E6 * ( Q*1E-3*gamma(vu)/(apm*local_Np*gamma(vu+bpm)) )**(1/bpm)  ! units = microns

    endif  
      
    Re = 0.5*D0*gamma(vu+3)/gamma(vu+2)
    
    
! ---------------------------------------------------------!
! // exponential                                           !
! ---------------------------------------------------------!
! :: N0 = intercept parameter (m^-4)
! :: ld = slope parameter (um)

  case(2)
  
    ! Np not specified (see if statement above) 
  
    if((abs(p1+1) > 1E-8) ) then   ! N0 has been specified, determine ld
    
        N0 = p1
    tmp1 = 1./(1.+bpm)
    ld = ((apm*gamma(1.+bpm)*N0)/(rho_a*Q*1E-3))**tmp1
    ld = ld/1E6                     ! set units to microns^-1
        
    elseif (abs(p2+1) > 1E-8) then  ! lambda=ld has been specified as default

        ld = p2     ! should have units of microns^-1 
    
    else
    
    PRINT *, 'Error: Must specify Np or default value ', &
         '(p1=No or p2=lambda) for Exponential distribution'
        stop
        
    endif

    Re = 1.5/ld ;
  
! ---------------------------------------------------------!
! // power law                                             !
! ---------------------------------------------------------!
! :: ahp = Ar parameter (m^-4 mm^-bhp)
! :: bhp = br parameter
! :: dmin_mm = lower bound (mm)
! :: dmax_mm = upper bound (mm)

  case(3)

    Re=0;  ! Not supporting LUT approach for power-law ...

    if(Np>0) then
    
        PRINT *, 'Variable Np not supported for ', &
         'Power Law distribution'
        stop
    endif
! ---------------------------------------------------------!
! // monodisperse                                          !
! ---------------------------------------------------------!
! :: D0 = particle diameter (um) == Re

  case(4)
  
        Re = p1
    
        if(Np>0) then
        PRINT *, 'Variable Np not supported for ', &
         'Monodispersed distribution'
        stop
    endif
    
! ---------------------------------------------------------!
! // lognormal                                             !
! ---------------------------------------------------------!
! :: N0 = total number concentration (m^-3)
! :: np = fixed number concentration (kg^-1)
! :: rg = mean radius (um)
! :: log_sigma_g = ln(geometric standard deviation)

  case(5)
  
        if( abs(local_p3+1) > 1E-8 ) then

            !set natural log width
            log_sigma_g = local_p3 
        else
            PRINT *, 'Error: Must specify a value for sigma_g ', &
             'when using a Log-Normal distribution'
            stop
        endif
     
    ! get rg ... 
    if( Np.eq.0 .AND. (abs(p2+1) > 1E-8) ) then ! use default value of rg
    
            rg = p2
              
        else
            if(Np>0) then
            
                local_Np=Np;
                
            elseif(abs(p2+1) < 1E-8) then
            
                local_Np=p1
            else
                PRINT *, 'Error: Must specify Np or default value ', &
                 '(p2=Rg or p1=Np) for Log-Normal distribution'
            endif
           
            log_sigma_g = p3
            tmp1 = (Q*1E-3)/(2.**bpm*apm*local_Np)
            tmp2 = exp(0.5*bpm**2.*(log_sigma_g))**2.      
            rg = ((tmp1/tmp2)**(1/bpm))*1E6
        endif
                
        Re = rg*exp(+2.5*(log_sigma_g**2)) 
        
  end select
  
  end subroutine calc_Re