module set_crtm_aerosolmod
!$$$ module documentation block
!           .      .    .                                       .
! module:   set_crtm_aerosolmod
!  prgmmr: todling          org: gmao                date: 2011-06-01
!
! abstract: module providing interface to set-crtm-aerosol procedures
!
! program history log:
!   2011-06-01  todling
!   2011-09-20  hclin   - separate na and na_crtm for p25 handling
!
! subroutines included:
!   sub Set_CRTM_Aerosol_
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

implicit none

private

public Set_CRTM_Aerosol

interface Set_CRTM_Aerosol
  module procedure Set_CRTM_Aerosol_
end interface

contains

subroutine Set_CRTM_Aerosol_ ( km, na, na_crtm, aero_name, aero_conc, rh, aerosol)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    Set_CRTM_Aerosol_
!   prgmmr: hclin          org: ncar/mmm                date: 2011-09-20
!
! abstract: Set the CRTM Aerosol object given GOCART aerosol properties.
!           This version based on the WRFCHEM implementation of GOCART.
!
! program history log:
!   2011-09-20  hclin   - 
!
!   input argument list:
!     km        : number of CRTM levels
!     na        : number of aerosols
!     na_crtm   : number of aerosols seen by CRTM
!     aero_name : GOCART aerosol names
!     aero_conc : aerosol concentration (Kg/m2)
!     rh        : relative humdity [0,1]
!     aerosol   : CRTM Aerosol object
!
!   output argument list:
!     aero_conc : aerosol concentration (Kg/m2)
!     aerosol   : CRTM Aerosol object
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

! USES:

  use kinds, only: i_kind,r_kind
  use constants, only: tiny_r_kind
  use mpimod, only: mype
  use CRTM_Aerosol_Define, only: CRTM_Aerosol_type
  use mpeu_util, only: getindex
  use crtm_module, only: SULFATE_AEROSOL,BLACK_CARBON_AEROSOL,ORGANIC_CARBON_AEROSOL,&
      DUST_AEROSOL,SEASALT_SSAM_AEROSOL,SEASALT_SSCM1_AEROSOL,SEASALT_SSCM2_AEROSOL,SEASALT_SSCM3_AEROSOL
  use mpimod, only : mype ! ajl

  implicit none

! !ARGUMENTS:

  integer(i_kind) , intent(in)    :: km                ! number of levels
  integer(i_kind) , intent(in)    :: na                ! number of aerosols
  integer(i_kind) , intent(in)    :: na_crtm           ! number of aerosols seen by CRTM
  character(len=*), intent(in)    :: aero_name(na)     ! [na]    GOCART aerosol names
  real(r_kind),     intent(inout) :: aero_conc(km,na)  ! [km,na] aerosol concentration (Kg/m2)
  real(r_kind),     intent(in)    :: rh(km)            ! [km]    relative humdity [0,1]

  type(CRTM_Aerosol_type), intent(inout) :: aerosol(na_crtm)! [na]   CRTM Aerosol object

  integer(i_kind) :: i, k
  integer(i_kind) :: indx_p25, indx_dust1, indx_dust2, indx_dust3, indx_dust4, indx_dust5
  integer(i_kind) :: indx_bc1, indx_oc1

  indx_bc1=-1; indx_oc1=-1; indx_dust1=-1; indx_dust2=-1
  indx_dust3=-1; indx_dust4=-1; indx_dust5=-1; indx_p25=-1
!  write(900+mype,*)'aero_name',aero_name
!  call flush(900+mype)

  indx_p25   = getindex(aero_name,'p25')
  indx_dust1 = getindex(aero_name,'dust1')
  indx_dust2 = getindex(aero_name,'dust2')
  indx_dust3 = getindex(aero_name,'dust3')
  indx_dust4 = getindex(aero_name,'dust4')
  indx_dust5 = getindex(aero_name,'dust5')
  indx_bc1   = getindex(aero_name,'bc1')
  indx_oc1   = getindex(aero_name,'oc1')

  do i = 1, na

     if ( trim(aero_name(i)) == 'p25' ) cycle

     ! assign aerosol type
!     write(900+mype,*)'i',i,aero_name(i),indx_p25,'km',km
!     write(900+mype,*)'i',i,aero_name(i),'km',km
!     call flush(900+mype)
     select case ( trim(aero_name(i)) )
        case ('sulf')
           aerosol(i)%type  = SULFATE_AEROSOL
        case ('bc1','bc2')
           aerosol(i)%type  = BLACK_CARBON_AEROSOL
        case ('oc1','oc2')
           aerosol(i)%type  = ORGANIC_CARBON_AEROSOL
        case ('dust1','dust2','dust3','dust4','dust5')
           aerosol(i)%type  = DUST_AEROSOL
        case ('seas1')
           aerosol(i)%type  = SEASALT_SSAM_AEROSOL
        case ('seas2')
           aerosol(i)%type  = SEASALT_SSCM1_AEROSOL
        case ('seas3')
           aerosol(i)%type  = SEASALT_SSCM2_AEROSOL
        case ('seas4')
           aerosol(i)%type  = SEASALT_SSCM3_AEROSOL
     end select
     if(aerosol(i)%type<1)then
         write(6,*)'illegal aerosol ',aero_name(i),'i',i
         call flush(6)
     endif

     if ( indx_p25 > 0 ) then
        ! partition p25 to dust1 and dust2
        if ( i == indx_dust1 ) then
           aero_conc(:,i) = aero_conc(:,i)+ 0.78_r_kind*aero_conc(:,indx_p25)
        endif
        if ( i == indx_dust2 ) then
           aero_conc(:,i) = aero_conc(:,i)+ 0.22_r_kind*aero_conc(:,indx_p25)
        endif
     endif

     ! crtm aerosol structure
!     write(900+mype,*)'km',km,'i',i
     do k = 1, km
        aerosol(i)%concentration(k) = max(tiny_r_kind, aero_conc(k,i))
!        write(940+mype,*)'set aersol ugkg ',i,' k ',k,aerosol(i)%concentration(k)
        ! calculate effective radius
           
        aerosol(i)%effective_radius(k) &
           = GOCART_Aerosol_size(i, aerosol(i)%type, rh(k))
!         if(i.ge.11.and.k.eq.1)then
!           write(900+mype,*)'reff salt',i,aerosol(i)%effective_radius(k)
!           call flush(900+mype)
!         endif
        ! 5 dust bins
        aerosol(indx_dust1)%effective_radius(k) = 0.55_r_kind
        aerosol(indx_dust2)%effective_radius(k) = 1.4_r_kind
        aerosol(indx_dust3)%effective_radius(k) = 2.4_r_kind
        aerosol(indx_dust4)%effective_radius(k) = 4.5_r_kind
        aerosol(indx_dust5)%effective_radius(k) = 8.0_r_kind
     enddo

  enddo  ! na
!  write(900+mype,*)' end na'
!  call flush(900+mype)

  contains

  function GOCART_Aerosol_size( kk, itype,  & ! Input
                                       eh ) & ! Input in 0-1
                           result( R_eff  )   ! in micrometer
  use crtm_aerosolcoeff, only: AeroC
  use mpimod, only : mype
  implicit none
!
!   modified from a function provided by Quanhua Liu
!
  integer(i_kind) ,intent(in) :: kk, itype
  real(r_kind)    ,intent(in) :: eh

  real(r_kind), parameter :: reff_seasalt(4) = Reshape( (/0.3_r_kind, 1.0_r_kind, 3.25_r_kind, 7.5_r_kind/), (/4/) )
  integer(i_kind) :: j1,j2,k
  real(r_kind)    :: h1
  real(r_kind)    :: R_eff
!  write(900+mype,*)'itype',itype,'DUST',DUST_AEROSOL,'kk',kk,'indxbc',indx_bc1,indx_oc1
!  call flush(6)
! ajl for check
  R_eff=0.0
  j1=-999
!  if(mype.eq.0)then
!     write(6,*)'itype',itype,'eh',eh,'kk',kk
!     write(6,*)'dust',DUST_AEROSOL,'BLACK ',BLACK_CARBON_AEROSOL,'indx_ bc1',indx_bc1
!     write(6,*)'OGRAN',ORGANIC_CARBON_AEROSOL,'indx_ oc1',indx_oc1
!  endif
  

  if ( itype==DUST_AEROSOL ) then
     return
  else if ( itype==BLACK_CARBON_AEROSOL .and. kk==indx_bc1 ) then
     R_eff = AeroC%Reff(1,itype )
!     if(mype.eq.0)write(6,*)'BC reff',R_eff
     return
  else if ( itype==ORGANIC_CARBON_AEROSOL .and. kk==indx_oc1 ) then
     R_eff = AeroC%Reff(1,itype )
!     if(mype.eq.0)write(6,*)'OC reff',R_eff
     return
  endif
!  if(mype.eq.0)then
!    write(6,*)'rh1',AeroC%RH(1) ,'n_rh',AeroC%n_RH ,AeroC%RH(AeroC%n_RH)
!  endif

  j2 = 0
!  if ( eh < AeroC%RH(1) ) then
! ajl if eh =rh(1) j1 undefined 
  if ( eh <= AeroC%RH(1) ) then
     j1 = 1
  else if ( eh > AeroC%RH(AeroC%n_RH) ) then
     j1 = AeroC%n_RH
  else
     do k = 1, AeroC%n_RH-1
        if ( eh <= AeroC%RH(k+1) .and. eh > AeroC%RH(k) ) then
           j1 = k
           j2 = k+1
           h1 = (eh-AeroC%RH(k))/(AeroC%RH(k+1)-AeroC%RH(k))
!  if(mype.eq.0)then
!     write(6,*)'inbetween k',k,AeroC%RH(k:k+1)
!  endif
           exit
        endif
     enddo
  endif
  if(j1<0)then
!      write(900+mype,*)'j1',j1,'eh',eh,'aero1',AeroC%RH(1),'nrh',AeroC%RH(AeroC%n_RH)
!      call flush(900+mype)
       R_eff=0.
       return
  endif

  if ( j2 == 0 ) then
     R_eff = AeroC%Reff(j1,itype )
!     write(900+mype,*)'j2 zero j1',j1,' R_eff',R_eff
!     call flush(900+mype)
  else
     
     R_eff = (1.0_r_kind-h1)*AeroC%Reff(j1,itype ) + h1*AeroC%Reff(j2,itype )
!     write(900+mype,*)'AERO reff ',j1,AeroC%Reff(j1,itype)
!     call flush(900+mype)
!     write(900+mype,*)'AERO reff ',j2,AeroC%Reff(j2,itype)
     !call flush(900+mype)
  endif
!  write(900+mype,*)'R_eff',R_eff,'eh',eh,'rh1',AeroC%RH(1),'rhmax',AeroC%RH(AeroC%n_RH),'j1',j1,j2
!  call flush(900+mype)

  return
  end function GOCART_Aerosol_size

end subroutine Set_CRTM_Aerosol_

end module set_crtm_aerosolmod
