!********************************************************
!  These subroutines are in a module so I can pass
!  allocated arrays back and forth
!********************************************************

MODULE retrieval_subs

CONTAINS

!********************************************************
!  SINGLE_SCATTER_LOOKUP SUBROUTINE
!********************************************************

subroutine single_scatter_lookup(mie_lwc_0p5,mie_re_0p5,mie_alpha_0p5_1mode,&
                                 mie_ext_0p5_1mode,mie_ssa_0p5_1mode,mie_g_0p5_1mode,&
                                 mie_lwc_2p1,mie_re_2p1,mie_alpha_2p1_1mode,&
                                 mie_ext_2p1_1mode,mie_ssa_2p1_1mode,mie_g_2p1_1mode,lwc,re,alpha)

use NETCDF

integer :: year,month,day,juldate
integer :: status,ncid,varid
integer :: lwc,re,alpha
character(len=100) :: filename

real, dimension(:),allocatable :: mie_lwc_0p5,mie_re_0p5,mie_alpha_0p5_1mode
real, dimension(:,:,:), allocatable :: mie_sca_0p5_1mode,mie_ext_0p5_1mode,&
  mie_g_0p5_1mode,mie_ssa_0p5_1mode

real, dimension(:), allocatable :: mie_lwc_2p1,mie_re_2p1,mie_alpha_2p1_1mode
real, dimension(:,:,:), allocatable :: mie_sca_2p1_1mode,mie_ext_2p1_1mode,&
  mie_g_2p1_1mode,mie_ssa_2p1_1mode

!print*,'in single scatter lookup'

filename='single_scatter_lookup_1mode_0p55.cdf'

!  Open the netcdf file
status=nf90_open(path=filename,mode=nf90_nowrite,ncid=ncid)

!  Read in dimensions
status=nf90_inq_dimid(ncid,"lwc",varid)
status=nf90_inquire_dimension(ncid,varid,len=lwc)
!print*,'lwc',lwc

status=nf90_inq_dimid(ncid,"re",varid)
status=nf90_inquire_dimension(ncid,varid,len=re)
!print*,'re',re

status=nf90_inq_dimid(ncid,"alpha",varid)
status=nf90_inquire_dimension(ncid,varid,len=alpha)
!print*,'alpha',alpha

!  Allocate the arrays to read in data from netcdf file
allocate ( mie_lwc_0p5(lwc))
allocate ( mie_re_0p5(re))
allocate ( mie_alpha_0p5_1mode(alpha))
allocate ( mie_sca_0p5_1mode(lwc,re,alpha))
allocate ( mie_ext_0p5_1mode(lwc,re,alpha))
allocate ( mie_g_0p5_1mode(lwc,re,alpha))
allocate ( mie_ssa_0p5_1mode(lwc,re,alpha))

!  Read in variables
status=nf90_inq_varid(ncid,"lwc_vector",varid)
status=nf90_get_var(ncid,varid,mie_lwc_0p5)

status=nf90_inq_varid(ncid,"re_vector",varid)
status=nf90_get_var(ncid,varid,mie_re_0p5)

status=nf90_inq_varid(ncid,"alpha_vector",varid)
status=nf90_get_var(ncid,varid,mie_alpha_0p5_1mode)

status=nf90_inq_varid(ncid,"sca_matrix",varid)
status=nf90_get_var(ncid,varid,mie_sca_0p5_1mode)

status=nf90_inq_varid(ncid,"ext_matrix",varid)
status=nf90_get_var(ncid,varid,mie_ext_0p5_1mode)

status=nf90_inq_varid(ncid,"g_matrix",varid)
status=nf90_get_var(ncid,varid,mie_g_0p5_1mode)

status=nf90_inq_varid(ncid,"ssa_matrix",varid)
status=nf90_get_var(ncid,varid,mie_ssa_0p5_1mode)

!  Close the netcdf file
status = nf90_close(ncid)

!print*,'extmatrix'
!print*,'ext',(mie_ext_0p5_1mode(1,1,i),i=1,alpha)
!print*,'ext',(mie_ext_0p5_1mode(i,1,1),i=1,lwc)
!print*,'ext',(mie_ext_0p5_1mode(1,i,1),i=1,re)
!print*,'lwc',(mie_lwc_0p5(i),i=1,lwc)
!print*,'re',(mie_re_0p5(i),i=1,re)

!***
!  Read in 2.1
!***

filename='single_scatter_lookup_1mode_2p1.cdf'

!  Open the netcdf file
status=nf90_open(path=filename,mode=nf90_nowrite,ncid=ncid)

!  Read in dimensions
status=nf90_inq_dimid(ncid,"lwc",varid)
status=nf90_inquire_dimension(ncid,varid,len=lwc)
!print*,'lwc',lwc

status=nf90_inq_dimid(ncid,"re",varid)
status=nf90_inquire_dimension(ncid,varid,len=re)
!print*,'re',re

status=nf90_inq_dimid(ncid,"alpha",varid)
status=nf90_inquire_dimension(ncid,varid,len=alpha)
!print*,'alpha',alpha

!  Allocate the arrays to read in data from netcdf file
allocate (mie_lwc_2p1(lwc))
allocate (mie_re_2p1(re))
allocate (mie_alpha_2p1_1mode(alpha))
allocate (mie_sca_2p1_1mode(lwc,re,alpha))
allocate (mie_ext_2p1_1mode(lwc,re,alpha))
allocate (mie_g_2p1_1mode(lwc,re,alpha))
allocate (mie_ssa_2p1_1mode(lwc,re,alpha))

!  Read in variables
status=nf90_inq_varid(ncid,"lwc_vector",varid)
status=nf90_get_var(ncid,varid,mie_lwc_2p1)

status=nf90_inq_varid(ncid,"re_vector",varid)
status=nf90_get_var(ncid,varid,mie_re_2p1)

status=nf90_inq_varid(ncid,"alpha_vector",varid)
status=nf90_get_var(ncid,varid,mie_alpha_2p1_1mode)

status=nf90_inq_varid(ncid,"sca_matrix",varid)
status=nf90_get_var(ncid,varid,mie_sca_2p1_1mode)

status=nf90_inq_varid(ncid,"ext_matrix",varid)
status=nf90_get_var(ncid,varid,mie_ext_2p1_1mode)

status=nf90_inq_varid(ncid,"g_matrix",varid)
status=nf90_get_var(ncid,varid,mie_g_2p1_1mode)

status=nf90_inq_varid(ncid,"ssa_matrix",varid)
status=nf90_get_var(ncid,varid,mie_ssa_2p1_1mode)

!  Close the netcdf file
status = nf90_close(ncid)

!print*,'re',(mie_re_2p1(i),i=1,re)

end subroutine single_scatter_lookup

!********************************************************
!  SUBROUTINE GET_LIQ_SINGLE_SCATTER
!********************************************************

subroutine get_liq_single_scatter(lwc,re,alpha,ext,ssa,g,mie_lwc,mie_re,mie_alpha,&
  mie_ext,mie_ssa,mie_g,n_lwc,n_re,n_alpha)

use interpolation

integer :: i,j,n_lwc,n_re,n_alpha
integer :: lwc_index,re_index,alpha_index
real :: alpha,ext,ssa,g,ext2,ssa2,g2

!double precision :: lwc,re
real :: lwc,re

real, dimension(:),allocatable :: mie_lwc,mie_re,mie_alpha
real, dimension(:,:,:), allocatable :: mie_ext,mie_ssa,mie_g

!***  Initialize these variables
ext=0.0
ssa=0.0
g=0.0

!***  Find the closest lwc value
j=1
do while (mie_lwc(j) .le. lwc .and. j .lt. n_lwc)
  !print*,mie_lwc(j),lwc,j,n_lwc
  j=j+1
enddo
lwc_index=j  !right index of the two values stradling lwc
!print*,'lwc_index',lwc_index
!print*,'lwc',lwc
!print*,'mie_lwc',(mie_lwc(i),i=1,n_lwc)

!***  Find the closest re value
j=1
do while (mie_re(j) .le. re .and. j .lt. n_re)
  !print*,mie_re(j),re,j,n_re
  j=j+1
enddo
re_index=j  !right index of the two values stradling re
!print*,'re_index',re_index
!print*,'re',re
!print*,'mie_re',(mie_re(i),i=1,n_re)

!***  Find the closest alpha value
j=1
do while (mie_alpha(j) .le. alpha .and. j .lt. n_alpha)
  !print*,mie_alpha(j),alpha,j,n_alpha
  j=j+1
enddo
alpha_index=j  !right index of the two values stradling alpha
!print*,'alpha_index',alpha_index
!print*,'alpha',alpha
!print*,'mie_alpha',(mie_alpha(i),i=1,n_alpha)

!***  Select out the closest values  right side
!ext=mie_ext(lwc_index, re_index, alpha_index)!*1.e5 !convert to 1/km from 1/cm
!ssa=mie_ssa(lwc_index, re_index, alpha_index)
!g=mie_g(lwc_index, re_index, alpha_index)

!print*,'***'
!print*,'closest',ext,ssa,g

!***  Select out the closest values  left side
!ext=mie_ext(lwc_index-1, re_index, alpha_index)!*1.e5 !convert to 1/km from 1/cm
!ssa=mie_ssa(lwc_index-1, re_index, alpha_index)
!g=mie_g(lwc_index-1, re_index, alpha_index)
!print*,'closest',ext,ssa,g

!***  Get the bilinear interpolated value
ext=interpolate(n_lwc,mie_lwc,n_re,mie_re,mie_ext(1:n_lwc,1:n_re,alpha_index),real(lwc),real(re))
ssa=interpolate(n_lwc,mie_lwc,n_re,mie_re,mie_ssa(1:n_lwc,1:n_re,alpha_index),real(lwc),real(re))
g=interpolate(n_lwc,mie_lwc,n_re,mie_re,mie_g(1:n_lwc,1:n_re,alpha_index),real(lwc),real(re))

!print*,ext
!print*,mie_ext(lwc_index-1,re_index-1,alpha_index),mie_ext(lwc_index,re_index-1,alpha_index)
!print*,mie_ext(lwc_index-1,re_index,alpha_index),mie_ext(lwc_index,re_index,alpha_index)
!print*,

!print*,ssa
!print*,mie_ssa(lwc_index-1,re_index-1,alpha_index),mie_ssa(lwc_index,re_index-1,alpha_index)
!print*,mie_ssa(lwc_index-1,re_index,alpha_index),mie_ssa(lwc_index,re_index,alpha_index)
!print*,

!print*,g
!print*,mie_g(lwc_index-1,re_index-1,alpha_index),mie_g(lwc_index,re_index-1,alpha_index)
!print*,mie_g(lwc_index-1,re_index,alpha_index),mie_g(lwc_index,re_index,alpha_index)
!print*,

!print*,ext2,ssa2,g2
!print*,ext,ssa,g

end subroutine get_liq_single_scatter

!***************************************************************

END MODULE retrieval_subs