!>
!!
!! @brief Module mo_EVA: provides volcanic aerosol optical  
!properties as a function of lat,  height, time, and wavelength
!!
!! @remarks
!!
!! @author Matthew Toohey 
!!
!! $ID: n/a$
!!
!! @par Origin
!!
!! @par Copyright
!!

module mo_EVA

   USE netcdf

   implicit none     

   real, parameter :: pi = 3.1415927
   real, parameter :: deg2rad = pi/180.

   ! --------------------------------------------------------------------------------------
   ! Sulfate box model parameters

   !----------------------------------------------------------------------------------------
   ! Input files

!   CHARACTER(len=50)    :: eruption_list_filename = 'Eruption_list_GMD_1960_2015.nc'
   CHARACTER(len=50)    :: eruption_list_filename = 'eruptions.nc'
   CHARACTER(len=50)    :: parameter_set_filename = 'EVAv1_parameter_set_v1.0.nc'
   CHARACTER(len=50)    :: Lookuptable_filename   = 'eva_Mie_lookuptables.nc'

   ! --------------------------------------------------------------------------------------
   ! Sulfate box model parameters

   REAL, SAVE ::                            &
       background_SO2_in,                   &  ! background SO2 injection
       tau_loss,                            &  ! timescale of loss in extratropical plumes
       tau_prod,                            &  ! timescale of effective SO2->AOD production in ET
       tau_mix,                             &  ! timescale of mixing between EQ and ET plumes
       tau_res,                             &  ! timescale of residual mass circ from EQ-to-ET
       seasonal_amp,                        &  ! amplitude of seasonal cycle 
       tropical_edge                           ! edge of tropical pipe

   ! ----------------------------------------------------------------------------------------
   ! Spatial shape function parameters

   REAL, SAVE ::                            &
       v_width_ET,                          &  ! vertical width of the extratropical plumes
       v_width_EQ,                          &  ! vertical width of the equatorial plume
       v_offset_EQ,                         &  ! vertical offset of the equatorial plume
       h_center_ET,                         &  ! center latitude of the extratropical plumes
       h_width_ET,                          &  ! horizontal width of the extratropical plumes
       h_width_EQ                              ! horizontal width of the equatorial plum
   REAL, SAVE, allocatable ::               &
       cline(:),                            &  ! describes height of max ext wrt latitude, length nlat
       shape_EQ(:,:),                       &  ! 2D equatorial shape function
       shape_NH(:,:),                       &  ! 2D NH shape function
       shape_SH(:,:),                       &  ! 2D SH shape function
       v_shape_EQ(:,:),                     &  ! 2D equatorial vertical shape
       v_shape_ET(:,:),                     &  ! 2D extratropical vertical shape
       h_shape_EQ(:),                       &  ! equatorial horizontal shape
       h_shape_NH(:),                       &  ! NH horizontal shape
       h_shape_SH(:)                           ! SH horizontal shape
   ! ----------------------------------------------------------------------------------------
   ! Sulfate to AOD scaling parameters

   REAL, SAVE  ::                           &
       so4_to_aod_linear,                   &  ! Scaling factor for small eruptions
       so4_to_aod_nonlin,                   &  ! coefficient for 2/3 power realtionship for larger eruptions 
       so4_to_reff,                         &  ! Sacling factor for effective radius
       nonlin_threshold                        ! threshold for nonlinear aod to so4 scaling

   ! -----------------------------------------------------------------------------------------
   ! Mie look up tables (lut) parameters 

   REAL, ALLOCATABLE ::                     & 
       lut_extrat(:,:),                     &  ! extinction relative to ext at 550nm
       lut_ssa   (:,:),                     &  ! single scattering albedo
       lut_asy   (:,:),                     &  ! asymmetry factor
       lut_reff  (:),                       &  ! effective radius
       lut_lambda(:)                           ! wavelength

   INTEGER ::                               &
       nlut_reff,                           &  ! number of effective radii in lut
       nlut_wl                                 ! number of wavelengths in lut 

   ! NAMELISTS
   NAMELIST /EVA_INPUT/ eruption_list_filename, parameter_set_filename, Lookuptable_filename

   ! ------------------------------------------------------------------------------------------
   ! PROGRAM control parameters

   LOGICAL, SAVE     :: eva_initialized = .FALSE.

   ! -------------------------------------------------------------------------------------------

contains

   subroutine init_eva(lat, z, nlat, nz)
      real, intent(in) :: lat(nlat), z(nz)
      integer, intent(in) :: nlat, nz 
      write(*,*) 'Initializing EVA'

      call read_parameter_set
      call read_aero_volc_tables
      call def_vert_centerline(lat,nlat) 

      OPEN (UNIT=10, FILE='eva_namelist', STATUS='OLD')
      READ (10, NML=EVA_INPUT)
      CLOSE (10)

      write(*,*) 'Center line defined'
      call eva_shape(lat,z)
      write(*,*) 'Shape defined'
      eva_initialized = .TRUE.

   end subroutine init_eva

   subroutine eva_aop_profile(lat,z,lat0,SO4,wl,nw,nlat,nz,ext,ssa,asy)

      real, intent(in), dimension(nlat) :: lat
      real, intent(in), dimension(nz) :: z
      real, intent(in) :: lat0
      real, intent(in) :: SO4(3)
      real, intent(in) :: wl(nw)
      integer, intent(in) :: nlat, nz, nw
      real, intent(out), dimension(nz,nw) :: ext, ssa, asy

      real, dimension(nz) :: ext550, reff  

      IF (.NOT.eva_initialized) CALL init_eva(lat, z, nlat, nz)

      call get_ext550(lat, z, lat0, SO4, nlat, nz, ext550)
      call get_reff(SO4,lat,lat0,nlat,nz,reff)
      call convert_aero_props_interp(ext550,reff,wl,nw,nz,ext,ssa,asy) 

   end subroutine eva_aop_profile    

   subroutine eva_ext_reff(lat,z,lat0,SO4,nlat,nz,ext550,reff)

      real, intent(in), dimension(nlat) :: lat
      real, intent(in), dimension(nz) :: z
      real, intent(in) :: lat0
      real, intent(in) :: SO4(3)
      integer, intent(in) :: nlat, nz
      real, intent(out), dimension(nz) :: ext550, reff

      IF (.NOT.eva_initialized) CALL init_eva(lat, z, nlat, nz)

      call get_ext550(lat, z, lat0, SO4, nlat, nz, ext550)
      call get_reff(SO4,lat,lat0,nlat,nz,reff)

   end subroutine eva_ext_reff

   subroutine get_ext550(lat,z,lat0,SO4,nlat,nz,ext550)

     ! get extinction profile at a given latitude lat0

     real, intent(in), dimension(nlat) :: lat
     real, intent(in), dimension(nz) :: z
     real, intent(in) :: lat0
     real, intent(in) :: SO4(3)
     integer, intent(in) :: nlat, nz
     real, intent(out), dimension(nz) ::ext550 

     ! local variables
     integer :: lat0_ind, i
     real :: sulfate_aerosol_mass(3)
     real :: aod, SO4_lat
     real :: v_shape(nz), v_shape_norm(nz)

!     nonlin_threshold = (so4_to_aod_nonlin/so4_to_aod_linear)**(3.)
!     jan:
     nonlin_threshold = (so4_to_aod_nonlin/so4_to_aod_linear)**(3./(1.0+3.0*0.1))

      ! find index of lat0 in lat
      do i=1,nlat
         if ( lat(i) .eq. lat0 ) then
            lat0_ind=i
          end if
      end do

      SO4_lat = h_shape_EQ(lat0_ind)*SO4(2) &
                      + h_shape_SH(lat0_ind)*SO4(1) &
                      + h_shape_NH(lat0_ind)*SO4(3) 
     
      if ( SO4_lat < nonlin_threshold ) then
         aod=so4_to_aod_linear*SO4_lat
      else
!         aod=so4_to_aod_nonlin*(SO4_lat**(2.0/3.0))
!         jan
         aod=so4_to_aod_nonlin*(SO4_lat**(2.0/3.0-0.1))
      end if

      v_shape=shape_EQ(lat0_ind,:)*SO4(2) &
                      + shape_SH(lat0_ind,:)*SO4(1) &
                      + shape_NH(lat0_ind,:)*SO4(3)
      if ( sum(v_shape) > 0 ) then
         v_shape_norm=v_shape/sum(v_shape)
      else
         v_shape_norm=0
      end if

      ext550=aod*v_shape_norm

   end subroutine get_ext550 

   subroutine get_reff(SO4,lat,lat0,nlat,nz,reff)
      real, intent(in) :: SO4(3), lat(nlat), lat0
      integer, intent(in) :: nlat, nz
      real, intent(out) :: reff(nz)

      integer :: i, lat0_ind

      do i=1,nlat
         if ( lat(i) .eq. lat0 ) then
            lat0_ind=i
         end if
      end do

      reff = so4_to_reff*(shape_EQ(lat0_ind,:)*SO4(2) &
                      + shape_SH(lat0_ind,:)*SO4(1) &
                      + shape_NH(lat0_ind,:)*SO4(3))**(1./3.)
      reff=max(maxval(reff),0.2)

   end subroutine get_reff

   subroutine def_vert_centerline(lat,nlat)
      real, intent(in) :: lat(nlat)
      integer, intent(in) :: nlat
      integer :: i
      real :: cline_lat(37), cline_raw(37)
      INTEGER ::          &
         iret         , & !< netCDF reading return variable
         ncid         , & !< netCDF file ID
         VarID            !< pointer to generic dimension in netCDF file

      allocate(cline(nlat))

      iret = nf90_open(TRIM(parameter_set_filename), NF90_NOWRITE, ncid)
      IF (iret /= NF90_NOERR) STOP 'Error in opening parameter file'

      iret = nf90_inq_varid(ncid, "cline_lat", VarID)
      iret = nf90_get_var(ncid, VarID, cline_lat(:)  , start=(/1/) ,count=(/37/))
      IF (iret /= NF90_NOERR) STOP 'Error in reading cline_lat'

      iret = nf90_inq_varid(ncid, "cline", VarID)
      iret = nf90_get_var(ncid, VarID, cline_raw(:)  , start=(/1/) ,count=(/37/))
      IF (iret /= NF90_NOERR) STOP 'Error in reading cline'

      do i=1,nlat
         cline(i)=interp(cline_lat,cline_raw,lat(i),37)
      end do

   end subroutine def_vert_centerline

   subroutine EVA_shape(lat,z)

      ! input/output
      real, intent(in), dimension(:) :: lat, z

      ! local variables
      real, dimension(:), allocatable :: sinlat
      real :: h_sin_center_ET, h_sin_width_ET, h_sin_width_EQ
      integer :: nlat, nz, i

      real, dimension(:,:), allocatable :: v_shape_EQ, v_shape_ET

      ! allocate arrays
      nlat = size(lat)
      nz=size(z)
      allocate(sinlat(nlat))
      allocate(h_shape_EQ(nlat))
      allocate(h_shape_NH(nlat))
      allocate(h_shape_SH(nlat))
      allocate(v_shape_EQ(nlat,nz))
      allocate(v_shape_ET(nlat,nz))
      allocate(shape_EQ(nlat,nz))
      allocate(shape_NH(nlat,nz))
      allocate(shape_SH(nlat,nz))

      ! Horizontal shape calculations

      sinlat=sin( lat * pi / 180. )
      h_sin_center_ET = sin( h_center_ET * pi / 180. )
      h_sin_width_ET = sin( h_width_ET * pi / 180. )
      h_sin_width_EQ = sin( h_width_EQ * pi / 180. )

      h_shape_EQ=exp( -(sinlat-0)**2. / (2. * h_sin_width_EQ )**2. )
      h_shape_NH=exp( -(sinlat-h_sin_center_ET)**2. / (2. * h_sin_width_ET)**2. )
      h_shape_SH=exp( -(sinlat+h_sin_center_ET)**2. / (2. * h_sin_width_ET)**2. )

      ! normalize (needed because extratropical shapes won't extend to
      ! full width. Also, normalizing by area puts shapes into units of
      ! km**-2)

      h_shape_EQ = h_shape_EQ / area_gmean(h_shape_EQ,lat)
      h_shape_NH = h_shape_NH / area_gmean(h_shape_NH,lat)
      h_shape_SH = h_shape_SH / area_gmean(h_shape_SH,lat)

      nlat=size(cline)

      ! build 3D shape

      do i=1,nlat
         v_shape_ET(i,:) = exp( -(z-cline(i))**2. / (2. * v_width_ET)**2.) &
                                  /(v_width_ET*sqrt(2.*pi))
         v_shape_EQ(i,:) = exp(-(z-cline(i)-v_offset_EQ)**2. / (2. * v_width_EQ)**2. ) &
                                  /(v_width_EQ*sqrt(2.*pi))

         ! normalize (assuming equally spaced z)
         v_shape_ET(i,:)=v_shape_ET(i,:)/sum(v_shape_ET(i,:)*(z(2)-z(1)))
         v_shape_EQ(i,:)=v_shape_EQ(i,:)/sum(v_shape_EQ(i,:)*(z(2)-z(1)))

         ! include horizontal shape
         shape_EQ(i,:)=v_shape_EQ(i,:)*h_shape_EQ(i)
         shape_NH(i,:)=v_shape_ET(i,:)*h_shape_NH(i)
         shape_SH(i,:)=v_shape_ET(i,:)*h_shape_SH(i)

      end do

      deallocate(sinlat)

   end subroutine EVA_shape

   subroutine sulfate3_timeseries(year,month,n,SO4)

   ! input/output
      real, intent(out), dimension(3,n) :: SO4
      real, intent(in), dimension(n) :: year, month
      integer, intent(in) :: n
      ! local variables
      integer :: ntime, i, j, k

      real, dimension(3,n) :: SO2_in, SO2, SO4_new
      real, dimension(:), allocatable :: erup_ssi, erup_lat, erup_hemi, erup_dur
      real :: C, tau_loss_EQ_use, tau_loss_ET_use, SO4_tot
      real, dimension(n) :: hemi_corr_nh, hemi_corr_sh
      integer, dimension(:), allocatable :: erup_region, erup_year, erup_month, erup_day
      real :: SHmix, NHmix, SHtrans, NHtrans, seas_asy
      real :: S_tot_t0, SO2_tot_t0, SO4_tot_t0, SO4_ET2EQ_ratio_t0

   INTEGER ::          &
       iret         , & !< netCDF reading return variable
       ncid         , & !< netCDF file ID
       VarID        , & !< pointer to generic dimension in netCDF file
       nerup        

      call read_parameter_set

      ntime=size(year)
      
      SO2_in=0
      SO2_in(1,:)=background_SO2_in/12.0/2.0
      SO2_in(3,:)=background_SO2_in/12.0/2.0
      SO2=0
      SO4=0
      SO4_new=0
      hemi_corr_nh=1.0
      hemi_corr_sh=1.0

      ! initialize SO4 based on background SO2_in (no spin-up needed)
      ! (this is very close to being right, SO4 still dips every so slighlty in
      ! first months)
      S_tot_t0=background_SO2_in/12.0*tau_loss
      SO4_tot_t0 = S_tot_t0/(1.0+tau_prod/tau_loss)
      SO2_tot_t0 = S_tot_t0/(1.0+tau_loss/tau_prod)
      SO4_ET2EQ_ratio_t0 = ((tau_mix*tau_res/2.0) + tau_mix*tau_loss + tau_loss*tau_res) / (tau_loss*tau_res)
      SO4(2,1)=SO4_tot_t0 / ( 2.0*SO4_ET2EQ_ratio_t0 + 1.0)
      SO4(1,1)=SO4_tot_t0 / ( 2.0 + 1.0/SO4_ET2EQ_ratio_t0)
      SO4(3,1)=SO4_tot_t0 / ( 2.0 + 1.0/SO4_ET2EQ_ratio_t0)
      SO2(2,1)=SO2_tot_t0 /  ( 2.0*SO4_ET2EQ_ratio_t0 + 1.0)
      SO2(1,1)=SO2_tot_t0 / ( 2.0 + 1.0/SO4_ET2EQ_ratio_t0) 
      SO2(3,1)=SO2_tot_t0 / ( 2.0 + 1.0/SO4_ET2EQ_ratio_t0)


      ! Read eruption list
      write(*,*) 'Reading file: ', TRIM(eruption_list_filename)
      iret = nf90_open(TRIM(eruption_list_filename), NF90_NOWRITE, ncid)
      IF (iret /= NF90_NOERR) STOP 'Error in opening eruption list file'
      iret = nf90_inq_dimid(ncid, "nerup", VarID)
      iret = nf90_inquire_dimension(ncid, VarID, len = nerup)

      allocate(erup_year(nerup))
      allocate(erup_month(nerup))
      allocate(erup_day(nerup))
      allocate(erup_lat(nerup))
      allocate(erup_ssi(nerup))
      allocate(erup_region(nerup))
      allocate(erup_hemi(nerup))
      allocate(erup_dur(nerup))

      iret = nf90_inq_varid(ncid, "year", VarID)
      iret = nf90_get_var(ncid, VarID, erup_year(:)  , start=(/1/) ,count=(/nerup/))
      IF (iret /= NF90_NOERR) STOP 'Error in reading eruption year'

      iret = nf90_inq_varid(ncid, "month", VarID)
      iret = nf90_get_var(ncid, VarID, erup_month(:)  , start=(/1/) ,count=(/nerup/))
      IF (iret /= NF90_NOERR) STOP 'Error in reading eruption month'

      iret = nf90_inq_varid(ncid, "day", VarID)
      iret = nf90_get_var(ncid, VarID, erup_day(:)  , start=(/1/) ,count=(/nerup/))
      IF (iret /= NF90_NOERR) STOP 'Error in reading eruption day'

      iret = nf90_inq_varid(ncid, "lat", VarID)
      iret = nf90_get_var(ncid, VarID, erup_lat(:)  , start=(/1/) ,count=(/nerup/))
      IF (iret /= NF90_NOERR) STOP 'Error in reading eruption lat'

      iret = nf90_inq_varid(ncid, "ssi", VarID)
      iret = nf90_get_var(ncid, VarID, erup_ssi(:)  , start=(/1/) ,count=(/nerup/))
      IF (iret /= NF90_NOERR) STOP 'Error in reading eruption ssi'

      iret = nf90_inq_varid(ncid, "hemi", VarID)
      iret = nf90_get_var(ncid, VarID, erup_hemi(:)  , start=(/1/) ,count=(/nerup/))
      IF (iret /= NF90_NOERR) STOP 'Error in reading eruption hemi'

      iret = nf90_inq_varid(ncid, "duration", VarID)
      IF (iret /= NF90_NOERR) THEN 
         erup_dur=0
      ELSE
         iret = nf90_get_var(ncid, VarID, erup_dur(:)  , start=(/1/) ,count=(/nerup/))
         IF (iret /= NF90_NOERR) STOP 'Error in reading eruption duration'
      END IF


      iret = nf90_close  (ncid)

      !write(*,*) erup_ssi
 
      do i=1,nerup
         IF ( erup_lat(i) < -tropical_edge ) THEN
            erup_region(i) = 1
         ELSE IF ( erup_lat(i) > tropical_edge ) THEN
            erup_region(i) = 3
         ELSE
            erup_region(i) = 2
         END IF
      end do

      ! Insert S injections into timeseries
      ! for each eruption in list
      do i=1,nerup
         ! if date is >= Year_start & date <= Year_end
         if (erup_year(i) .ge. year(1) .and. erup_year(i) .le. year(ntime) ) then
            ! loop through length of year/month to find matching date
            do j=1,ntime
               if (year(j) .eq. erup_year(i) .and. month(j) .eq. erup_month(i)) then
                  ! set SO2_in for this date = elist_SO2
                  if (erup_dur(i) .gt. 0) then
                     SO2_in(erup_region(i),j:j+erup_dur(i))=SO2_in(erup_region(i),j:j+erup_dur(i)) + erup_ssi(i)/erup_dur(i)
                  else
                     SO2_in(erup_region(i),j)=SO2_in(erup_region(i),j) + erup_ssi(i)
                  end if
                  seas_asy=0.4*seasonal_amp*cos(2.0*pi*(month(j)+4.0)/12.0)+1 ! a rough approximation of what the 
                                                                 ! seasonal paramterized transport produces
                  if ( erup_hemi(i) > 0.0 ) then
                     if ( abs(seas_asy - erup_hemi(i)) > 0.1 ) then
                        if ( erup_hemi(i) < 1.0 ) then
                           do k=j,min(j+17,ntime)
                              hemi_corr_nh(k)=seas_asy/erup_hemi(i) 
                           end do
                        else
                           do k=j,min(j+17,ntime)
                              hemi_corr_sh(k) = erup_hemi(i)/seas_asy
                           end do
                        end if
                     end if
                  end if
               end if
            end do
         end if
      end do
   
      ! Run box-model calculations to produce SO4 timeseries
     
      !SO2(:,1)=SO2_in(:,1)

      ! SO2 is converted into SO4
      do i=2,ntime
         SO2(1,i)=SO2(1,i-1)*exp(-1/tau_prod)*exp(-1/tau_loss) +SO2_in(1,i)
         SO2(2,i)=SO2(2,i-1)*exp(-1/tau_prod)*exp(-1/tau_loss) +SO2_in(2,i)
         SO2(3,i)=SO2(3,i-1)*exp(-1/tau_prod)*exp(-1/tau_loss) +SO2_in(3,i)
         SO4_new(1,i)=SO2(1,i-1)*(1-exp(-1/tau_prod))
         SO4_new(2,i)=SO2(2,i-1)*(1-exp(-1/tau_prod))
         SO4_new(3,i)=SO2(3,i-1)*(1-exp(-1/tau_prod))
      end do

      ! SO4 grows and mixes

      do i=2,ntime
         ! compute mixing and transport timescales with seasonal
         ! variability. Seasonality modeled as cosine, with maxima in
         ! January (NH) and July (SH).
         SHmix = (1.0/(tau_mix*hemi_corr_sh(i)))*(seasonal_amp*cos(2*pi*mod(i-1+6,12)/12)+1)
         NHmix = (1.0/(tau_mix*hemi_corr_nh(i)))*(seasonal_amp*cos(2*pi*mod(i-1,12)/12)+1)
         SHtrans = (1.0/(hemi_corr_sh(i)*tau_res))*(seasonal_amp*cos(2*pi*mod(i-1+6,12)/12)+1)
         NHtrans = (1.0/(hemi_corr_nh(i)*tau_res))*(seasonal_amp*cos(2*pi*mod(i-1,12)/12)+1)

         ! accelerate loss for large sulfur loadings (special case for
         ! 1257 eruption: relationship based roughly on ECHAM5-HAM
         ! simulations of large eruptions. But eventually this should be
         ! replace, the loss rate would be better related to Reff, this
         ! requires merging the sulfate and Reff growth routines
         SO4_tot=sum(SO4(:,i-1))
         if (SO4_tot > 10) then
            tau_loss_EQ_use=((tau_loss-6.0)/0.3679)*exp(-SO4_tot/10.0)+6.0
            tau_loss_ET_use=((tau_loss-6.0)/0.3679)*exp(-SO4_tot/10.0)+6.0
         else
            tau_loss_EQ_use=tau_loss
            tau_loss_ET_use=tau_loss
         end if


         SO4(2,i) = (SO4(2,i-1)+SO4_new(2,i))*exp(-1/tau_loss_EQ_use) &
                          ! add production and subtract STE loss
                  - (SO4(2,i-1) - SO4(3,i-1))*SHmix &
                           ! subtract loss by mixing to SH
                  - (SO4(2,i-1)- SO4(1,i-1))*NHmix &
                           ! subtract loss by mixing to NH
                   - SO4(2,i-1)*NHtrans - SO4(2,i-1)*SHtrans
                           ! subtract loss by residual circ to both hemispheres

         SO4(1,i) = (SO4(1,i-1)+SO4_new(1,i))*exp(-1/tau_loss_ET_use) &
                           ! add production and subtract STE loss
                  + (SO4(2,i-1) - SO4(1,i-1))*SHmix &
                           ! add mixing from EQ
                   + SO4(2,i-1)*SHtrans
                           ! add residual circ from EQ

         SO4(3,i) = (SO4(3,i-1)+SO4_new(3,i))*exp(-1/tau_loss_ET_use) &
                           ! add production and subtract STE loss
                  + (SO4(2,i-1) - SO4(3,i-1))*NHmix &
                           ! add mixing from EQ
                   + SO4(2,i-1)*NHtrans
                           ! add residual circ from EQ
 
      end do
      
      !write(*,*) SO4(:,3)

      deallocate(erup_ssi)
      deallocate(erup_region)
      deallocate(erup_year)
      deallocate(erup_month)

   end subroutine sulfate3_timeseries

   function area_gmean(x,lat)        
      real, dimension(:) :: x, lat          
      real, dimension(:), allocatable :: w, rlat
      real :: area_gmean, lat_bounds(2)
      integer :: i, nlat
      real, parameter :: pi = 3.1415927
      nlat=size(lat)
      allocate(w(nlat))
      allocate(rlat(nlat))
      rlat = (pi/180.)*lat
      w=cos(rlat)
      w = w / sum( w )
      area_gmean = sum ( w * x )

      deallocate(w)
      deallocate(rlat) 

   end function area_gmean

   subroutine read_parameter_set
      INTEGER ::          &
         iret         , & !< netCDF reading return variable
         ncid         , & !< netCDF file ID
         VarID            !< pointer to generic dimension in netCDF file

    write(*,*) 'Reading parameter set: ', TRIM(parameter_set_filename)     
 
      iret = nf90_open(TRIM(parameter_set_filename), NF90_NOWRITE, ncid)
      IF (iret /= NF90_NOERR) STOP 'Error in opening parameter file'
      iret = nf90_inq_varid(ncid, "version", VarID)

      iret = nf90_inq_varid(ncid, "background_SO2_in", VarID)
      iret = nf90_get_var(ncid, VarID, background_SO2_in)
      IF (iret /= NF90_NOERR) STOP 'Error in reading background_SO2_in'

      iret = nf90_inq_varid(ncid, "tau_loss", VarID)
      iret = nf90_get_var(ncid, VarID, tau_loss)
      IF (iret /= NF90_NOERR) STOP 'Error in reading tau_loss'

      iret = nf90_inq_varid(ncid, "tau_prod", VarID)
      iret = nf90_get_var(ncid, VarID, tau_prod)
      IF (iret /= NF90_NOERR) STOP 'Error in reading tau_prod'

      iret = nf90_inq_varid(ncid, "tau_mix", VarID)
      iret = nf90_get_var(ncid, VarID, tau_mix)
      IF (iret /= NF90_NOERR) STOP 'Error in reading tau_mix'

      iret = nf90_inq_varid(ncid, "tau_res", VarID)
      iret = nf90_get_var(ncid, VarID, tau_res)
      IF (iret /= NF90_NOERR) STOP 'Error in reading tau_res'

      iret = nf90_inq_varid(ncid, "seasonal_amp", VarID)
      iret = nf90_get_var(ncid, VarID, seasonal_amp)
      IF (iret /= NF90_NOERR) STOP 'Error in reading seasonal_amp'

      iret = nf90_inq_varid(ncid, "tropical_edge", VarID)
      iret = nf90_get_var(ncid, VarID, tropical_edge)
      IF (iret /= NF90_NOERR) STOP 'Error in reading tropical_edge'

      iret = nf90_inq_varid(ncid, "v_width_ET", VarID)
      iret = nf90_get_var(ncid, VarID, v_width_ET)
      IF (iret /= NF90_NOERR) STOP 'Error in reading v_width_ET'

      iret = nf90_inq_varid(ncid, "v_width_EQ", VarID)
      iret = nf90_get_var(ncid, VarID, v_width_EQ)
      IF (iret /= NF90_NOERR) STOP 'Error in reading v_width_EQ'

      iret = nf90_inq_varid(ncid, "v_offset_EQ", VarID)
      iret = nf90_get_var(ncid, VarID, v_offset_EQ)
      IF (iret /= NF90_NOERR) STOP 'Error in reading v_offset_EQ'

      iret = nf90_inq_varid(ncid, "h_center_ET", VarID)
      iret = nf90_get_var(ncid, VarID, h_center_ET )
      IF (iret /= NF90_NOERR) STOP 'Error in reading h_center_ET'

      iret = nf90_inq_varid(ncid, "h_width_ET", VarID)
      iret = nf90_get_var(ncid, VarID, h_width_ET)
      IF (iret /= NF90_NOERR) STOP 'Error in reading h_width_ET'

      iret = nf90_inq_varid(ncid, "h_width_EQ", VarID)
      iret = nf90_get_var(ncid, VarID, h_width_EQ)
      IF (iret /= NF90_NOERR) STOP 'Error in reading h_width_EQ'

      iret = nf90_inq_varid(ncid, "so4_to_aod_linear", VarID)
      iret = nf90_get_var(ncid, VarID, so4_to_aod_linear)
      IF (iret /= NF90_NOERR) STOP 'Error in reading so4_to_aod_linear'

      iret = nf90_inq_varid(ncid, "so4_to_aod_nonlin", VarID)
      iret = nf90_get_var(ncid, VarID, so4_to_aod_nonlin)
      IF (iret /= NF90_NOERR) STOP 'Error in reading so4_to_aod_nonlin'

      iret = nf90_inq_varid(ncid, "so4_to_reff", VarID)
      iret = nf90_get_var(ncid, VarID, so4_to_reff)
      IF (iret /= NF90_NOERR) STOP 'Error in reading so4_to_reff'

      !write(*,*) tau_loss, tau_prod, so4_to_reff

   end subroutine read_parameter_set


   subroutine read_aero_volc_tables
      ! read in lookup tables from netcdf file
      INTEGER ::          &
         iret         , & !< netCDF reading return variable
         ncid         , & !< netCDF file ID
         VarID        , & !< pointer to generic dimension in netCDF file
         nreff        , & !< length of effective radius vector in lookuptable
         nwl          
      write(*,*) 'Reading look up table:', TRIM(Lookuptable_filename)
      ! read in the lookup tables
      !
      iret = nf90_open(TRIM(Lookuptable_filename), NF90_NOWRITE, ncid)
      IF (iret /= NF90_NOERR) STOP 'Error in opening lookup tables file'
      iret = nf90_inq_dimid(ncid, "reff", VarID)
      iret = nf90_inquire_dimension(ncid, VarID, len = nreff)
      iret = nf90_inq_dimid(ncid, "wl", VarID)
      iret = nf90_inquire_dimension(ncid, VarID, len = nwl)

      allocate(lut_lambda(nwl))
      allocate(lut_reff(nreff))
      allocate(lut_extrat(nreff,nwl))
      allocate(lut_ssa(nreff,nwl))
      allocate(lut_asy(nreff,nwl))

      iret = nf90_inq_varid(ncid, "reff", VarID)
      iret = nf90_get_var(ncid, VarID, lut_reff(:)  , start=(/1/)  ,count=(/nreff/))
      IF (iret /= NF90_NOERR) STOP 'Error in reading effective radii'
      iret = nf90_inq_varid(ncid, "wl", VarID)
      iret = nf90_get_var(ncid, VarID, lut_lambda(:)  , start=(/1/)  ,count=(/nwl/))
      IF (iret /= NF90_NOERR) STOP 'Error in reading wavelengths'
      iret = nf90_inq_varid(ncid, "extrat", VarID)
      iret = nf90_get_var(ncid, VarID, lut_extrat(:,:), start=(/1,1/),count=(/nreff,nwl/))
      IF (iret /= NF90_NOERR) STOP 'Error in reading extrat'
      iret = nf90_inq_varid(ncid, "ssa", VarID)
      iret = nf90_get_var(ncid, VarID, lut_ssa(:,:), start=(/1,1/),count=(/nreff,nwl/))
      IF (iret /= NF90_NOERR) STOP 'Error in reading ssa'
      iret = nf90_inq_varid(ncid, "asy", VarID)
      iret = nf90_get_var(ncid, VarID, lut_asy(:,:), start=(/1,1/),count=(/nreff,nwl/))
      IF (iret /= NF90_NOERR) STOP 'Error in reading asy'

      iret = nf90_close  (ncid)

      nlut_reff=nreff
      nlut_wl=nwl
      !write(*,*) lut_lambda
   end subroutine read_aero_volc_tables

   subroutine convert_aero_props_interp(ext550, reff, wl, nw, nz, ext, ssa, asy)
      ! Conversion of ext550 and Reff to ext(wl), ssa(wl) and asy(wl)
      ! based on interpolation in wavelength to look up table values
      ! To do: we will replace the look up tables for this option with
      !  something more useful for other models, i.e., higher res

      real, intent(in) :: ext550(nz), reff(nz), wl(nw)
      integer, intent(in) :: nw, nz
      real, intent(out) :: ext(nz,nw), ssa(nz,nw), asy(nz,nw)
      integer :: i, j

      ! For each height and wavelength, find ext, ssa and asym by bilinear interpolation
      !write(*,*) lut_extrat
      do i=1,nz
         do j=1,nw
            ext(i,j) = ext550(i) * interp2(lut_reff, lut_lambda, lut_extrat, reff(i), wl(j), nlut_reff, nlut_wl) 
            ssa(i,j) = interp2(lut_reff, lut_lambda, lut_ssa, reff(i), wl(j), nlut_reff, nlut_wl) 
            asy(i,j) = interp2(lut_reff, lut_lambda, lut_asy, reff(i), wl(j), nlut_reff, nlut_wl) 
         end do
      end do

   end subroutine convert_aero_props_interp 

   function interp(xin, yin, xout, n)
      ! simple routine to linearly interpolate
      ! xin should be monotonically increasing
      ! xout should be a single value
      ! For now, for xout values outside the range of xin, 
      !  the function returns the "nearest neighbour" 

      real :: xin(n), yin(n)
      integer :: n
      real :: interp, xout, slope
      integer :: i

      i=1
      if (xout < xin(1)) then
         interp=yin(1)
         return
      end if

      if (xout > xin(n)) then
         interp=yin(n)
         return
      end if

      do while (xin(i) < xout)
         i=i+1
      end do

      slope = (yin(i)-yin(i-1))/(xin(i)-xin(i-1))
      interp = yin(i-1) + slope * ( xout - xin(i-1) )

   end function interp


   function interp2(xin, yin, zin, xout, yout, nx, ny)
   ! simple routine to bilinearly interpolate
   ! xin, yin should be monotonically increasing
   ! xout, yout should be single values
   ! For now, for xout,yout values outside the range of xin, yin u
   !  the function returns 0 

      real, intent(in) :: xin(nx), yin(ny), zin(nx,ny), xout, yout
      real :: interp2, z_x_y1, z_x_y2
      integer :: i, j, nx, ny

      if (xout < xin(1)) then
         interp2=0
         return
      end if

      if (xout > xin(nx)) then
         interp2=0
         return
      end if
      if (yout < yin(1)) then
         interp2=0
         return
      end if

      if (yout > yin(ny)) then
         interp2=0 
         return
      end if 
        
      j=1 
      do while (yin(j) < yout)
         j=j+1
      end do
    
      ! first interpolate in x at y1 and y2
      z_x_y1 = interp(xin, zin(:,j-1), xout, nx)
      z_x_y2 = interp(xin, zin(:,j), xout, nx)

      ! then interpolate in y    
      interp2 = interp(yin(j-1:j), (/ z_x_y1, z_x_y2 /), yout, 2)

   end function interp2

 
end module mo_EVA
