ecearth3/sources/tm5mp/proj/cb05/optics.F90

#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      OPTICS
!
! !DESCRIPTION: Optics module to calculate optical depth from m7 output,
!               based on the AOP_Package op Michael Kahnert.
!               
!               The optics code should be used in the following way (see
!               examples in photolysis.F90, ecearth_optics.F90, and
!               user_output_aerocom.F90): 
!               
!               1) prepare the optics calculation by calling 
!                  'OPTICS_INIT' --> lookuptables etc. 
!                  this routine reads the wavelengths, lookupable and calculates
!                  refr. indices at these wavelengths.
!               
!               2) allocate AOP arrays aop_out_ext/a/g with 
!                  (n_gridcells, n_wavelengths, n_split), with:
!                  'n_split' = 1 for (split == .false.) or 
!                  'n_split' = 11 for (split == .true. ), ie. 
!                    partial contributions by 
!                     Total, SO4, BC, OC, SS, DU, NO3, Water, Fine, Fine Dust, Fine SS
!                  additional fields are also provided for split==.true. in 
!                   aop_out_add, which has to have size 
!                   (n_gridcells, n_wavelengths, n_add) with nadd = 2 for 
!                    surface PM10 dry extinction and surface PM10 dry absorption
!               3) Compute AOP for current conditions by calling 'OPTICS_AOP_GET'
!               4) done: 'OPTICS_DONE'
!
!  IMPORTANT:   *) Skipped the ZOOM options! (temporary) 
!               *) OC is actually POM. Remember converting OC to POM 
!                  when sending it to optics_calculate_aop
!\\
!\\
! !INTERFACE: 
!
MODULE OPTICS
  !
  ! !USES:
  !
  use GO,              only : gol, goErr, goPr
  Use mo_aero_m7,      only : nmod
  Use global_data,     only : rcfile
  use GO,              only : TrcFile, Init, Done, ReadRc
  Use Dims,            only : nregions

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public  :: optics_init,   optics_done  ! init/done methods for a wavelendep object
  public  :: optics_aop_get              ! compute aerosol optical properties
  !
  ! !PRIVATE DATA MEMBERS:
  !
  type(TrcFile)      :: rcF  
  character(len=256) :: lookuptable, refractive_indices
  !
  ! !PUBLIC TYPES:
  !
  ! wavelength type (to be used by methods using the optics)
  type, public :: wavelendep 
     real  :: wl             ! user requested wavelength    unit = um (e.g. 0.550)
     real, dimension(7) :: n ! SO4, BC, OC, SOA, SS, DU, WATER
     real, dimension(7) :: k ! SO4, BC, OC, SOA, SS, DU, WATER
     logical :: split, insitu
  end type wavelendep
  !
  ! !PRIVATE TYPES:
  !
  ! AOP input type and field
  type aopi
     real, dimension(nmod) :: SO4, NO3, BC, OC, SOA, SS, DU, H2O, numdens, rg, rgd
  end type aopi

  type(aopi), dimension(:), allocatable :: aop_in 

  ! characteristics of the lookup-tables
  integer, parameter                      :: n_rii = 15
  integer, parameter                      :: n_rir = 40
  integer, parameter                      :: n_x   = 100
  real,dimension(n_rii)                   :: lkval        ! -log img part refr. index
  real,dimension(n_rii)                   :: kval         ! img part refr. index, 10^(-lkval)
  real,dimension(n_rir)                   :: n1r          ! real part refr. index

  Real, Dimension(n_x)                    :: xs
  Real, Dimension(n_x,n_rir,n_rii),Target :: cext_159, a_159, g_159
  real, dimension(n_x,n_rir,n_rii),Target :: cext_200, a_200, g_200

  character(len=*), parameter  ::  mname = 'Optics'

  ! lookup table information for interpolation to wavelength properties
  Integer, Parameter                :: opacdim       =   61 ! Wavelength dimension of OPAC data
  Integer, Parameter                :: echamhamdim   =   49 ! Wavelength dimension of ECHAM-HAM data
  Integer, Parameter                :: segelsteindim = 1261 ! Wavelength dimension of Segelstein data
  Real, Dimension(:,:), allocatable :: opac, echamham, segelstein

  !
  ! !REVISION HISTORY: 
  !    Implemented by Maarten Krol 03/2009
  !    Extended by Joost Aan de Brugh
  ! 
  !    6 Feb 2011 - Achim Strunk - complete revision 
  !                              - worked on DECOUPLING optics routines 
  !                                from (optics)_output, in order to 
  !                                (re)establish a flexible way of using it
  !                              - remaining parts have been moved to
  !                                (the new routine) user_output_optics
  !
  !   24 Jun 2011 - Achim Strunk - added NO3 explicitly in order to allow
  !                                for a slightly better split of the
  !                                optical properties of (SO4+NO3)
  !   20 Jun 2012 - P. Le Sager  - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OPTICS_INIT
  !
  ! !DESCRIPTION: Initialize a wavelendep object: Input/lookuptables are
  !               opened and read in, and wavelength dependent parameters are
  !               calculated.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE OPTICS_INIT( nwav, wdep, status )
    !
    ! !USES:
    !
    USE MDF, ONLY : MDF_OPEN, MDF_NETCDF, MDF_READ, MDF_Get_Var, MDF_Close, MDF_Inq_VarID
    !
    ! !INPUT PARAMETERS:
    !
    integer,                           intent(in)    :: nwav
    !
    ! !OUTPUT PARAMETERS:
    !
    type(wavelendep), dimension(nwav), intent(inout) :: wdep
    integer,                           intent(out)   :: status
    !
    ! !REMARKS:  read by each core (FIXME?)
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer                     :: fid, varid, i
    character(len=*), parameter :: rname = mname//'/Optics_Init'

    ! --- begin --------------------------------

    allocate( opac      (7, opacdim      ) )
    allocate( echamham  (5, echamhamdim  ) )
    allocate( segelstein(3, segelsteindim) )

    ! get filenames from rcfile
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'optics.lookuptable', lookuptable, status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'optics.refractive_indices', refractive_indices, status )
    IF_NOTOK_RETURN(status=1)

    ! read lookup table
    CALL MDF_Open( TRIM(lookuptable), MDF_NETCDF, MDF_READ, fid, status )
    IF_NOTOK_RETURN(status=1)
    
    CALL MDF_Inq_VarID( fid, 'mr', varid, status )
    IF_NOTOK_RETURN(status=1)
    CALL MDF_Get_Var( fid, varid, n1r, status )
    IF_NOTOK_RETURN(status=1)

    CALL MDF_Inq_VarID( fid, 'mi', varid, status )
    IF_NOTOK_RETURN(status=1)
    CALL MDF_Get_Var( fid, varid, kval, status )
    IF_NOTOK_RETURN(status=1)

    lkval = -log10(kval)

    CALL MDF_Inq_VarID( fid, 'x', varid, status )
    IF_NOTOK_RETURN(status=1)
    CALL MDF_Get_Var( fid, varid, xs, status )
    IF_NOTOK_RETURN(status=1)
    
    CALL MDF_Inq_VarID( fid, 'ext_159', varid, status )
    IF_NOTOK_RETURN(status=1)
    CALL MDF_Get_Var( fid, varid, cext_159, status )
    IF_NOTOK_RETURN(status=1)

    CALL MDF_Inq_VarID( fid, 'ext_200', varid, status )
    IF_NOTOK_RETURN(status=1)
    CALL MDF_Get_Var( fid, varid, cext_200, status )
    IF_NOTOK_RETURN(status=1)
    
    CALL MDF_Inq_VarID( fid, 'a_159', varid, status )
    IF_NOTOK_RETURN(status=1)
    CALL MDF_Get_Var( fid, varid, a_159, status )
    IF_NOTOK_RETURN(status=1)

    CALL MDF_Inq_VarID( fid, 'a_200', varid, status )
    IF_NOTOK_RETURN(status=1)
    CALL MDF_Get_Var( fid, varid, a_200, status )
    IF_NOTOK_RETURN(status=1)

    CALL MDF_Inq_VarID( fid, 'g_159', varid, status )
    IF_NOTOK_RETURN(status=1)
    CALL MDF_Get_Var( fid, varid, g_159, status )
    IF_NOTOK_RETURN(status=1)

    CALL MDF_Inq_VarID( fid, 'g_200', varid, status )
    IF_NOTOK_RETURN(status=1)
    CALL MDF_Get_Var( fid, varid, g_200, status )
    IF_NOTOK_RETURN(status=1)

    CALL MDF_Close( fid, status )
    IF_NOTOK_RETURN(status=1)

    ! Read refractive indices
    CALL MDF_Open( TRIM(refractive_indices), MDF_NETCDF, MDF_READ, fid, status )
    IF_NOTOK_RETURN(status=1)
    
    ! Get Opac data
    CALL MDF_Inq_VarID( fid, 'Opac', varid, status )
    IF_NOTOK_RETURN(status=1)
    CALL MDF_Get_Var( fid, varid, opac, status )
    IF_NOTOK_RETURN(status=1)

    ! Get ECHAM-HAM data
    CALL MDF_Inq_VarID( fid, 'ECHAM-HAM', varid, status )
    IF_NOTOK_RETURN(status=1)
    CALL MDF_Get_Var( fid, varid, echamham, status )
    IF_NOTOK_RETURN(status=1)

    ! Get Segelstein data
    CALL MDF_Inq_VarID( fid, 'Segelstein', varid, status )
    IF_NOTOK_RETURN(status=1)
    CALL MDF_Get_Var( fid, varid, segelstein, status )
    IF_NOTOK_RETURN(status=1)
    
    CALL MDF_Close( fid, status )
    IF_NOTOK_RETURN(status=1)

    call optics_wavelen_init( wdep, status )
    IF_ERROR_RETURN( status=1 )

    write(gol,*) 'Optical parameters:'; call goPr
    do i = 1, nwav
      write(gol,*) 'Wavelength :', wdep(i)%wl ; call goPr
      write(gol,*) ' SO4 (real/img) :', wdep(i)%n(1), wdep(i)%k(1) ; call goPr
      write(gol,*) ' BC  (real/img) :', wdep(i)%n(2), wdep(i)%k(2) ; call goPr
      write(gol,*) ' OC  (real/img) :', wdep(i)%n(3), wdep(i)%k(3) ; call goPr
      write(gol,*) ' SOA (real/img) :', wdep(i)%n(4), wdep(i)%k(4) ; call goPr
      write(gol,*) ' SS  (real/img) :', wdep(i)%n(5), wdep(i)%k(5) ; call goPr
      write(gol,*) ' DU  (real/img) :', wdep(i)%n(6), wdep(i)%k(6) ; call goPr
      write(gol,*) ' H2O (real/img) :', wdep(i)%n(7), wdep(i)%k(7) ; call goPr
    enddo
       
    ! Done
    ! -------------
    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)

    deallocate( opac, echamham, segelstein )

    status = 0

  END SUBROUTINE OPTICS_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OPTICS_DONE
  !
  ! !DESCRIPTION: dummy for now
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE OPTICS_DONE( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !    6 Feb 2011 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    ! ok
    status = 0
  END SUBROUTINE OPTICS_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OPTICS_WAVELEN_INIT
  !
  ! !DESCRIPTION: Initialise parameters which are depending on given 
  !               wavelengths.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE OPTICS_WAVELEN_INIT( wdep, status )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    ! wavelength properties (wavelength itself and real/img part of refractive index)
    type(wavelendep), intent(inout), dimension(:) :: wdep 
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,          intent(out)                 :: status
    !
    ! !REVISION HISTORY: 
    !    6 Feb 2011 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/optics_wavelen_init'
    
    integer :: nwl, i, idx
    real    :: wl, h
    real    :: nscale, kscale

    ! Real part n and imaginary part k of refractive index 
    ! for each wavelength:

    !>>> TvN
    ! The refractive index for 'sulfate' is taken from 
    ! the OPAC database (Hess et al., 1998; 
    ! Koepke et al., MPI report no. 243, 1997).
    ! The value used is that of 'sulfate solution',
    ! i.e. particles consisting of 75% of sulfuric acid (H2SO4),
    ! based on Fenn et al. (chapter 18 in Handbook of Geophysics
    ! and Space Environment, 1985).
    ! This is in line with Kinne et al. (JGR, 2003),
    ! who write that refractive indices for sulfate
    ! are usually based on 75% sulfuric acid solution.
    ! Actually, the OPAC refractive index agrees very well
    ! with the expression given by Boucher and Anderson (JGR, 1995)
    ! for H2SO4 at visible wavelengths.
    ! Thus, the OPAC data can be considered to apply to pure H2SO4,
    ! which is consistent of our application of mixing rules.
    
    ! For BC, the OPAC refractive index is scaled to one of
    ! the values proposed by Bond and Bergstrom (2006),
    ! valid at 550 nm (see their Table 5).
    ! Selected values for high, medium and low absorption are
    ! 1.95 + 0.79 i, 1.85 + 0.71 i, and 1.75 + 0.63 i, respectively.
    ! We select the low-absorption value,
    ! because it produces results in best agreement
    ! with AAOD from MODIS Collection 6 Deep Blue.
    ! Note that the medium-absorption value 1.85 + 0.71 i 
    ! was used in ECHAM simulations by Stier et al. (2007),
    ! and found to give reasonable results.
    ! Lowenthal et al. (Atmos. Environ., 2000)
    ! give a value of 1.90 + 0.6 i.
    ! The reference value used in the scaling
    ! is the OPAC value at 550 nm: 1.75 + 0.44 i.
    ! It would be better to include the scaling 
    ! in the input file.
    !nscale = 1.0 ! for using OPAC
    !kscale = 1.0 ! for using OPAC
    !nscale = 1.95/1.75
    !kscale = 0.79/0.44
    nscale = 1.85/1.75
    kscale = 0.71/0.44
    !nscale = 1.75/1.75
    !kscale = 0.63/0.44
    !<<< TvN

    nwl = size(wdep)
    do i = 1, nwl
       wl = wdep(i)%wl

       ! Interpolate Opac data
       findOpac: Do idx = 1,opacdim
          If(opac(1,idx) .gt. wl) exit findOpac
       End Do findOpac
       If(idx .gt. opacdim) then
          idx = opacdim
          h = 1.0
       Else If(idx .eq. 1) then
          idx = 2
          h = 0.0
       Else
          h = (wl-opac(1,idx-1))/(opac(1,idx)-opac(1,idx-1))
       End If

       ! SO4
       wdep(i)%n(1) = opac(2,idx-1)+h*(opac(2,idx)-opac(2,idx-1))
       wdep(i)%k(1) = opac(3,idx-1)+h*(opac(3,idx)-opac(3,idx-1))
       ! BC
       !>>> TvN
       !wdep(i)%n(2) = opac(4,idx-1)+h*(opac(4,idx)-opac(4,idx-1))
       !wdep(i)%k(2) = opac(5,idx-1)+h*(opac(5,idx)-opac(5,idx-1))
       wdep(i)%n(2) = nscale * ( opac(4,idx-1)+h*(opac(4,idx)-opac(4,idx-1)) )
       wdep(i)%k(2) = kscale * ( opac(5,idx-1)+h*(opac(5,idx)-opac(5,idx-1)) )
       !<<< TvN
       ! SS 
       wdep(i)%n(5) = opac(6,idx-1)+h*(opac(6,idx)-opac(6,idx-1))
       wdep(i)%k(5) = opac(7,idx-1)+h*(opac(7,idx)-opac(7,idx-1))

       ! Interpolate ECHAM-HAM data
       findEchamham: Do idx = 1,echamhamdim
          If(echamham(1,idx) .gt. wl) exit findEchamham
       End Do findEchamham
       If(idx .gt. echamhamdim) then
          idx = echamhamdim
          h = 1.0
       Else If(idx .eq. 1) then
          idx = 2
          h = 0.0
       Else
          h = (wl-echamham(1,idx-1))/(echamham(1,idx)-echamham(1,idx-1))
       End If

       ! OC 
       !>>> TvN
       ! The 'ECHAM-HAM' data currently used for POM
       ! are based on the data from Fenn et al. (1985) 
       ! for the 'water soluble' component,
       ! but at a reduced number of wavelengths up to 15 um.
       ! It would be better to use the original table
       ! in the input file.
       !<<< TvN
       wdep(i)%n(3) = echamham(2,idx-1)+h*(echamham(2,idx)-echamham(2,idx-1))
       wdep(i)%k(3) = echamham(3,idx-1)+h*(echamham(3,idx)-echamham(3,idx-1))


       ! SOA
       ! >>TB
       ! For now (March 2017) we use OC indices for SOA
       ! <<TB
       wdep(i)%n(4) = echamham(2,idx-1)+h*(echamham(2,idx)-echamham(2,idx-1))
       wdep(i)%k(4) = echamham(3,idx-1)+h*(echamham(3,idx)-echamham(3,idx-1))

       ! DU
       wdep(i)%n(6) = echamham(4,idx-1)+h*(echamham(4,idx)-echamham(4,idx-1))
       wdep(i)%k(6) = echamham(5,idx-1)+h*(echamham(5,idx)-echamham(5,idx-1))

       ! Interpolate Segelstein data
       findSegelstein: Do idx = 1,segelsteindim
          If(segelstein(1,idx) .gt. wl) exit findSegelstein
       End Do findSegelstein
       If(idx .gt. segelsteindim) then
          idx = segelsteindim
          h = 1.0
       Else If(idx .eq. 1) then
          idx = 2
          h = 0.0
       Else
          h = (wl-segelstein(1,idx-1))/(segelstein(1,idx)-segelstein(1,idx-1))
       End If

       ! H2O
       wdep(i)%n(7) = segelstein(2,idx-1)+h*(segelstein(2,idx)-segelstein(2,idx-1))
       wdep(i)%k(7) = segelstein(3,idx-1)+h*(segelstein(3,idx)-segelstein(3,idx-1))
    enddo

    status=0

  END SUBROUTINE OPTICS_WAVELEN_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OPTICS_GET
  !
  ! !DESCRIPTION: Main optical properties routine. Here the interpolated 
  !               values for extinction coefficient, single scattering
  !               albedo and assymetry parameter are retrieved and returned. 
  !\\
  !\\
  ! !INTERFACE:
  !
  PURE SUBROUTINE OPTICS_GET(m_eff, xg, Cext, a, g, cext_table, a_table, g_table, status, gol )
    !
    ! !INPUT PARAMETERS:
    !
    complex,                intent(in)    :: m_eff
    real,                   intent(in)    :: xg
    real, dimension(:,:,:), intent(in)    :: cext_table, a_table, g_table
    !
    ! !OUTPUT PARAMETERS:
    !
    real,             intent(out)   :: Cext, a, g
    integer,          intent(inout) :: status
    character(len=*), intent(inout) :: gol
    !
    ! !REVISION HISTORY: 
    !   15 Jul 2010 - P. Le Sager - Added check on i_n (out-of-bound)
    !    6 Feb 2011 - Achim Strunk -
    !
    ! !REMARKS:
    !   - By using "pure" subroutine, we cannot use any of the goPr, goErr,
    !      IF_NOTOK_RETURN, ...  routines. Traceback is limited:
    !      1) only one statement can be recorded in variable gol
    !      2) the calling routine should issue the "call goErr"
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real    :: n, k, n1, k1, dk1, dk2, lk
    real    :: modrad, modrad1, dr, dr1, slope, cext1, a1, g1
    integer :: i
    integer :: i_n, i_k, i_knew

    character(len=*), parameter ::  rname = mname//'/Optics_Get'

    ! --- begin --------------------------------

    status = 0

    n=real(m_eff)
    k=imag(m_eff)
    if(k.lt.kval(1))then
       k1=kval(1)
       i_knew=1
    elseif(k.gt.kval(15))then
       k1=kval(15)
       i_knew=15
    else
       get_k: do i=2,15
          if(k.lt.kval(i))then
             dk1=k-kval(i-1)
             dk2=kval(i)-k
             if(dk1.le.dk2)then
                k1=kval(i-1)
                i_knew=i-1
             else
                k1=kval(i)
                i_knew=i
             endif
             exit get_k
          endif
       enddo get_k
    endif
    lk=-log10(k1)
    !#n1=float(int(50*n+0.5))/50.
    !do i_n = 1, n_rir
    !   if (abs(n1 - n1r(i_n)) < 1e-4) exit
    !enddo

    !i_n = 1 + int((n-1.12)*50+0.5)
    !AJS: I guess n is a number on the (increasing) n1r axis; search nearest index:
    i_n = size(n1r)
    do i = 1, size(n1r)
       if ( n <= n1r(i) ) then
          i_n = i
          exit
       end if
    end do
    if ( i_n > 1 ) then
       if ( abs(n-n1r(i_n-1)) < abs(n-n1r(i_n)) ) i_n = i_n - 1
    else if ( i_n < n_rir ) then
       if ( abs(n-n1r(i_n+1)) < abs(n-n1r(i_n)) ) i_n = i_n + 1
    end if

    do i_k = 1, n_rii
       if (abs(lk - lkval(i_k)) < 1e-4) exit
    enddo

    !     ! PLS - test : ik can be determined without the preceding loop "do i_k = 1, n_rii"
    !     if (i_k.ne.i_knew) then
    !        status = 1
    !        write (gol,'(" PLSPLS  ik NE iknew = ",2(i2,2x))')i_k,i_knew 
    !     endif

    ! following the "get_k" loop above, the only way to get into this is to
    ! have a NaN for k in the first place ? 
    if (i_k > n_rii) then
       status = 1
       write(gol,*)'FATAL ERROR: i_k value outside LUT'
       !        write (gol,'("    lk = ",E16.4)')lk 
       !        write (gol,'("    k1 = ",E16.4)')k1 
       !        write (gol,'("     k = ",E16.4)')k  
       !        write (gol,'("     n = ",E16.4)')n  

       !        do i_n = 1, 15
       !           write (gol,'("    lkval(",i2,") = ",E16.4)')i_n,lkval(i_n)
       !           call goPr
       !           write (gol,'("     kval(",i2,") = ",E16.4)')i_n,kval(i_n)
       !           call goPr
       !        enddo

       return
    endif

    ! Added check (15-7-2010 - P. Le Sager)
    if (i_n > n_rir) then
       status = 1
       write(gol,*)'FATAL ERROR: i_n value out of range'
       return
    endif


    !>>> TvN
    ! Since xs(1) equals zero a problem may occur when xg is negative,
    ! because modrad.gt.xg would become true for i=1 in that case,
    ! while modrad1, Cext1, a1 and g1 are not yet set.
    ! It is not clear if negative xg values can ever occur,
    ! but if they do that should be prevented when calculating rg.
    ! 
    ! However, the problem may perhaps already occur
    ! when xg equals zero, because of the finite machine precision.
    !
    ! In any case, it is desired to initialize modrad1, Cext1, a1 and g1.
    ! Cext1, a1 and g1 should be set to their table entries for i=1,
    ! which are all zero:
    Cext1 = cext_table(1,i_n,i_k)
    a1 = a_table(1,i_n,i_k)
    g1 = g_table(1,i_n,i_k)
    ! modrad1 can be initialized to any value different from xs(1),
    ! to prevent division by zero.
    modrad1=xs(1)-9.99e-4
    ! This combination will force slope to zero
    ! and Cext, a and g to the table entries for i=1 (zero),
    ! in case modrad.gt.xg is true for i=1.
    !<<< TvN
    get_values: do i = 1, n_x
       modrad = xs(i)
       a = a_table(i,i_n, i_k)
       g = g_table(i,i_n, i_k)
       cext = cext_table(i,i_n, i_k)
       ! With small concentrations, like in the stratosphere, M7 does not trust its radii.
       ! See m7_averageproperties -> zinsvol. The M7-radius goes to zero
       ! Modrad may never be larger than xg the first time. Modrad1 is not set,
       ! will be zero (or something worse) and dr will be zero (or something worse).
       ! slope is demolished, makeing all values NaN. Therefore, it is modrad .gt. xg
       ! instead of modrad .ge. xg
       !
       ! PLS-TRANSLATION - It boils down to: "on the first iteration of this loop, if modrad=xg and
       !  you go into the if-statement below, you are in trouble, because modrad1 can be
       !  anything. To prevent that, replace GE with GT to avoid going into the if-statement at the
       !  first iteration."
       !
       if(modrad.gt.xg)then
          dr=modrad-modrad1
          dr1=xg-modrad1
          slope=(Cext-Cext1)/dr
          Cext=Cext1+slope*dr1
          slope=(a-a1)/dr
          a=a1+slope*dr1
          slope=(g-g1)/dr
          g=g1+slope*dr1
          exit get_values 
       endif
       modrad1=modrad
       Cext1=Cext
       a1=a
       g1=g
    enddo get_values

  END SUBROUTINE OPTICS_GET
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    GET_REFR_IDX
  !
  ! !DESCRIPTION: Compute refractive index of internally mixed aerosols by use
  !               of effective medium theory for the size-dependent aerosol
  !               mixtures assumed in M7.
  !\\
  !\\
  ! !INTERFACE:
  !
  PURE SUBROUTINE GET_REFR_IDX(wdep, SO4, BC, OC, SOA, SS, DU, water, mode, m_eff, status, gol)
    !
    ! !USES:
    !
    use binas,      only: rol
    use chem_param, only: ss_density, dust_density, carbon_density
    use chem_param, only: pom_density, so4_density, soa_density 
    !
    ! !INPUT PARAMETERS:
    !
    type(wavelendep), intent(in)    :: wdep   ! wavelength properties (wavelength, re/img part of refractive index)
    real,             intent(in)    :: SO4, BC, OC,SOA    ! mass mixing ratios or concentrations of sulphate, black carbon, organic carbon
    real,             intent(in)    :: SS,  DU, water ! sea salt, dust, and water
    integer,          intent(in)    :: mode   ! mode number (M7)
    !
    ! !OUTPUT PARAMETERS:
    !
    complex,          intent(out)   :: m_eff  ! effective refractive index of mixture
    integer,          intent(out)   :: status
    character(len=*), intent(inout) :: gol
    !
    ! !REVISION HISTORY: 
    !   12 Aug 2008 - Michael Kahnert, SMHI
    !    6 Feb 2011 - Achim Strunk -
    !
    ! !REMARKS:
    ! - see remarks of OPTICS_GET subroutine about PURE routine
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/get_refr_idx'

    ! refractive indices
    complex :: m_SO4, m_BC, m_OC, m_SOA, m_SS, m_DU, m_water

    ! volume fractions
    real    :: v_SO4, v_BC, v_OC, v_SOA, v_SS, v_DU, v_water, water_iv

    integer :: i_test
    real    :: vtot, v2
    complex :: m00, m0, m1, m2

    real    :: rpls, ipls

    ! --- begin --------------------------------

    ! Get the refractive indices from the lookup-tables and put them into complex numbers.
    m_so4    = cmplx( wdep%n(1),wdep%k(1) ) ! H2-SO4 + NH4NO3
    m_bc     = cmplx( wdep%n(2),wdep%k(2) ) ! BC
    m_oc     = cmplx( wdep%n(3),wdep%k(3) ) ! POM
    m_soa    = cmplx( wdep%n(4),wdep%k(4) ) ! SOA
    m_ss     = cmplx( wdep%n(5),wdep%k(5) ) ! SS
    m_du     = cmplx( wdep%n(6),wdep%k(6) ) ! DU
    m_water  = cmplx( wdep%n(7),wdep%k(7) ) ! Water

    status = 0
    !     no mixing for mode=6,7:
    if(mode.ge.6)then
       m_eff=m_DU
       return
    endif

    !     compute volume fractions:
    v_SO4=0.
    v_BC=0.
    v_OC=0.
    v_SOA=0.
    v_SS=0.
    v_DU=0.
    v_water=0.
    vtot=0.

    ! Added sanity check (15-7-2010 - P. Le Sager) : Avoid negative water
    ! mixing ratio!
    ! The bruggeman logically assumes that v_water is between 0 and 1, but
    ! this is never checked in the call chain : 
    ! ECEarth_Optics_Step -> calculate_aop -> get_refr_idx [here] ->
    ! Bruggeman
    ! We do it here, with a warning since it reflects a problem upstream:
    if(water.lt.0.0)then
       !write (gol,*)" WARNING - [Get_refr_idx] : negative relative humidity..." ; call goPr
       !write (gol,*)" WARNING - [Get_refr_idx] : .....set to 0" ; call goPr
       water_iv=0.0
    else
       water_iv=water
    endif

    if(mode.le.4)then
       v_SO4=SO4/so4_density
       v_water=water_iv/rol
       vtot=vtot+v_SO4+v_water
    endif
    if(mode.le.5)then
       !v_OC=OC/pom_density
       v_SOA=SOA/soa_density
       vtot=vtot+v_SOA
    end if
    if(mode.ge.2.and.mode.le.5)then
       v_BC=BC/carbon_density
       v_OC=OC/pom_density
       vtot=vtot+v_BC+v_OC
    endif
    if(mode.ge.3.and.mode.le.4)then
       v_SS=SS/ss_density
       vtot=vtot+v_SS
    endif
    if(mode.ge.3.and.mode.ne.5)then
       v_DU=DU/dust_density
       vtot=vtot+v_DU
    endif
    ! If vtot is zero, we will get 0.0/0.0's causing NaNs. In that case, the 
    ! refractive index does not matter and will be set to (1.0,1.0e-9). The 
    ! reason not to take (1.0,0.0) is that someone with humour might take 
    ! the logarithm of the imaginary part. Dust particles get their usual 
    ! refractive index, because they already returned m_DU. But that does 
    ! not matter, because there are zero aerosols in this case.
    If (vtot .le. 1e-20) Then
       m_eff = Cmplx(1.0,1.0e-9)
    Else
       v_SO4=v_SO4/vtot
       v_BC=v_BC/vtot
       v_OC=v_OC/vtot
       v_SOA=v_SOA/vtot
       v_SS=v_SS/vtot
       v_DU=v_DU/vtot
       v_water=v_water/vtot

       !-----------------------------------------------------------------------
       !     effective medium computations for each mode
       !-----------------------------------------------------------------------
       if(mode.eq.1)then
          !        Bruggeman mixing rule for SO4 OC and water:
          m1=m_SO4
          m2=m_SOA
          vtot=v_SO4+v_SOA
          v2=v_SOA/vtot
          call Bruggeman(m1,m2,v2,m0)
          m1=m0
          m2=m_water
          v2=v_water
          call Bruggeman(m1,m2,v2,m0)
       elseif(mode.eq.2)then
          !        iterative Bruggeman mixing rule for SO4, OC, and water:
          m1=m_SO4
          m2=m_OC
          vtot=v_SO4+v_OC
          v2=v_OC/vtot
          call Bruggeman(m1,m2,v2,m00)
          m1=m00
          m2=m_SOA
          vtot=vtot+v_SOA
          v2=v_SOA/vtot
          call Bruggeman(m1,m2,v2,m00)
          m1=m00
          m2=m_water
          vtot=vtot+v_water
          v2=v_water/vtot
          call Bruggeman(m1,m2,v2,m00)
          !        Maxwell-Garnett mixing rule for BC inclusions:
          m1=m00
          m2=m_BC
          v2=v_BC
          call Maxwell_Garnett(m1,m2,v2,m0)
       elseif(mode.eq.3.or.mode.eq.4)then
          !        iterative Bruggeman mixing rule for SO4, OC, SS, and water:
          m1=m_SO4
          m2=m_OC
          vtot=v_SO4+v_OC
          if ( vtot < TINY( vtot ) ) then
             v2=0.0
          else
             v2=v_OC/vtot
          end if
          call Bruggeman(m1,m2,v2,m00)
          m1=m00
          m2=m_SOA
          vtot=vtot+v_SOA
          if ( vtot < TINY( vtot ) ) then
             v2=0.0
          else
             v2=v_SOA/vtot
          end if
          call Bruggeman(m1,m2,v2,m00)
          m1=m00
          m2=m_SS
          vtot=vtot+v_SS
          if ( vtot < TINY( vtot ) ) then
             v2=0.0
          else
             v2=v_SS/vtot
          end if
          call Bruggeman(m1,m2,v2,m00)
          m1=m00
          m2=m_water
          vtot=vtot+v_water
          v2=v_water/vtot
          call Bruggeman(m1,m2,v2,m00)
          !        iterative Maxwell-Garnett mixing rule for BC and dust 
          !        inclusions:
          m1=m00
          m2=m_BC
          vtot=vtot+v_BC
          if ( vtot < TINY( vtot ) ) then
             v2=0.0
          else
             v2=v_BC/vtot
          end if
          call Maxwell_Garnett(m1,m2,v2,m00)
          m1=m00
          m2=m_DU
          v2=v_DU
          call Maxwell_Garnett(m1,m2,v2,m0)
       elseif(mode.eq.5)then

          m1=m_SOA
          m2=m_OC
          vtot=v_SOA+v_OC
          v2=v_OC/vtot
          call Bruggeman(m1,m2,v2,m00)
 !        Maxwell-Garnett mixing rule for BC inclusions:
          m1=m00
          m2=m_BC
          v2=v_BC
          call Maxwell_Garnett(m1,m2,v2,m0)          


       endif
       m_eff = m0
    End If

    ! Debug : trap for a NAN (13-7-2010 - P. Le Sager)
    rpls=real(m_eff)
    ipls=imag(m_eff)
    IF ((rpls.NE.rpls).or.(ipls.NE.ipls)) then
       status = 1
       write (gol,'(" GET_REFR_IDX-NAN:  ", 3(E16.4,2x),i4,2x,7(E16.4,2x))') rpls, ipls, vtot, mode,&
            &    SO4,BC,OC,SOA,SS,DU,water
    endif


  END SUBROUTINE GET_REFR_IDX
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    BRUGGEMAN
  !
  ! !DESCRIPTION: Compute effective refractive index of a mixture of 2 components
  !               by use of the Bruggeman mixing rule
  !\\
  !\\
  ! !INTERFACE:
  !
  PURE SUBROUTINE BRUGGEMAN(m1,m2,v2,m0)
    !
    ! !INPUT PARAMETERS:
    !
    complex, intent(in)  :: m1,m2
    real,    intent(in)  :: v2
    !
    ! !OUTPUT PARAMETERS:
    !
    complex, intent(out) :: m0
    !
    ! !REVISION HISTORY: 
    !   12 Aug 2008 - Michael Kahnert, SMHI
    !    6 Feb 2011 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    complex m1s,m2s,mt
    real fac1,fac2

    if(v2.eq.1.0)then
       m0=m2
       return
    elseif(v2.eq.0.0)then
       m0=m1
       return
    endif

    fac1=2.-3.*v2
    fac2=3.*v2-1.
    m1s=m1**2
    m2s=m2**2

    mt=m1s*fac1+m2s*fac2
    m0=1./16.*mt**2+0.5*m1s*m2s
    m0=csqrt(m0)
    m0=m0+0.25*mt
    m0=csqrt(m0)

  END SUBROUTINE BRUGGEMAN
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    MAXWELL_GARNETT
  !
  ! !DESCRIPTION: Compute effective refractive index for a medium consisting of
  !               a matrix with refractive index m1 and inclusions with refractive
  !               index m2 and volume fraction v2 by use of the Maxwell-Garnett 
  !               mixing rule.
  !\\
  !\\
  ! !INTERFACE:
  !
  PURE SUBROUTINE MAXWELL_GARNETT( m1, m2, v2, m0)
    !
    ! !INPUT PARAMETERS:
    !
    complex, intent(in)  :: m1, m2
    real,    intent(in)  :: v2
    !
    ! !OUTPUT PARAMETERS:
    !
    complex, intent(out) :: m0
    !
    ! !REVISION HISTORY: 
    !   12 Aug 2008 - Michael Kahnert, SMHI
    !    6 Feb 2011 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    complex :: m1s,m2s
    real :: fac1,fac2

    fac1=3.0-2.0*v2
    fac2=3.0-v2
    m1s=m1**2
    m2s=m2**2

    m0=m2s*(fac1*m1s+2.0*v2*m2s)/(v2*m1s+fac2*m2s)
    m0=csqrt(m0)

  END SUBROUTINE MAXWELL_GARNETT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OPTICS_CALCULATE_AOP
  !
  ! !DESCRIPTION: This routine writes on aop_out_* (module wide parameters) 
  !               the retrieved aerosol properties. The caller has to ensure
  !               that these fields have been allocated properly.
  !   IMPORTANT:  OC is actually POM. 
  !               Remember converting OC to POM when sending it to this method.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE OPTICS_CALCULATE_AOP( nwl, wdep, lvec, ecearth_units, &
                                   aop_out_ext, aop_out_a, aop_out_g, status, aop_out_add )
    !
    ! !USES:
    !
    use binas,        only : rol, twopi
    use chem_param,   only : ss_density, dust_density, carbon_density
    use chem_param,   only : pom_density, so4_density, soa_density
    Use mo_aero_m7,   only : nmod, cmedr2mmedr, sigma
    !
    ! !INPUT PARAMETERS:
    !
    integer,          intent(in)                 :: nwl
    type(wavelendep), intent(in), dimension(nwl) :: wdep   ! wavelength depending properties 
    integer,          intent(in)                 :: lvec
    logical,          intent(in)                 :: ecearth_units
    !
    ! !OUTPUT PARAMETERS:
    !
    Real, Dimension(:,:,:), intent(out)           :: aop_out_ext ! extinctions
    Real, Dimension(:,:),   intent(out)           :: aop_out_a   ! single scattering albedo
    Real, Dimension(:,:),   intent(out)           :: aop_out_g   ! assymetry factor
    integer,                intent(out)           :: status
    Real, Dimension(:,:,:), intent(out), optional :: aop_out_add ! additional parameters
    !
    ! !REVISION HISTORY: 
    !   12 Aug 2008 - Michael Kahnert, SMHI
    !    6 Feb 2011 - Achim Strunk -
    !   23 Jan 2013 - Narcisa Banda - included OMP (PLS: commented in TM5-MP)
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real,    dimension(:), allocatable :: NCsca, incext
    real                               :: Cexti, ai, gi, NCscai, xg
    complex                            :: m_eff
    real,    dimension(:), allocatable :: lnsigma

    integer :: i, imode, ivec!, lss, lee
    integer :: statusomp
    Logical :: coarse
    real    :: totvoldry, modfrac
    !    real                 :: t1, t2 , omp_get_wtime
    Real, Dimension(:,:,:), Pointer :: cext_table, a_table, g_table

    ! --- const ------------------------------

    character(len=*), parameter :: rname = mname//'/optics_calculate_AOP'

    ! --- begin --------------------------------


    allocate( NCsca ( lvec ) )
    !    allocate( NCscai( lvec ) )
    !    allocate( Cexti ( lvec ) )
    !    allocate( ai    ( lvec ) )
    !    allocate( gi    ( lvec ) )
    !    allocate( xg    ( lvec ) )
    !    allocate( m_eff ( lvec ) )
    allocate( incext( lvec ) )

    !t1= omp_get_wtime() 

    status = 0 

    !pls!   !$OMP PARALLEL &
    !pls!   !$OMP  default (none) &
    !pls!   !$OMP  shared  (nwl, lvec, NCsca, incext, &
    !pls!   !$OMP           wdep, aop_in, aop_out_Ext, aop_out_g, aop_out_a, aop_out_add, sigma, cmedr2mmedr,&
    !pls!   !$OMP           cext_200, a_200, g_200, cext_159, a_159, g_159, status, gol) &
    !pls!   !$OMP  private (lnsigma, i, imode, modfrac, coarse, totvoldry, lss, lee, ivec, & 
    !pls!   !$OMP           NCscai, Cexti, ai, gi, xg, m_eff, cext_table, a_table, g_table, statusomp )
    !pls!       call my_omp_domain (1, lvec, lss, lee)



    !     Sulphate based on OPAC (Hess et al., 1998):

    !=======================================================================
    !     Get refractive indices of each component at the given wavelength:
    !=======================================================================
    do i = 1, nwl   ! loop over wavelengths

       if( wdep(i)%split .or. wdep(i)%insitu ) then 
          allocate( lnsigma( nmod ) )
          lnsigma = log(sigma)
       end if

       aop_out_Ext(:,i,:)  = 0.0
       aop_out_g  (:,i)    = 0.0
       aop_out_a  (:,i)    = 0.0
       if( wdep(i)%insitu ) aop_out_add(:,i,:)  = 0.0 
       NCsca      (:)      = 0.0

       do imode = 1, 7        ! loop over M7 modes
          coarse = (imode .eq. 4 .or. imode .eq. 7)

          If (coarse) then
             cext_table => cext_200
             a_table => a_200
             g_table => g_200
          Else
             cext_table => cext_159
             a_table => a_159
             g_table => g_159
          End if


          !=======================================================================
          !     Compute effective refractive index of internally mixed aerosols
          !     for each grid cell and mode:
          !=======================================================================
          do ivec = 1, lvec 
             !write(1001,*) aop_in(ivec)%SO4(imode), aop_in(ivec)%NO3(imode)
             !write(1002,*) aop_in(ivec)%SOA(imode)
             !write(1003,*) aop_in(ivec)%BC(imode), aop_in(ivec)%OC(imode)
             !write(1004,*) aop_in(ivec)%SS(imode), aop_in(ivec)%DU(imode)


             call get_refr_idx( wdep(i), &
                  aop_in(ivec)%SO4(imode) + aop_in(ivec)%NO3(imode), & ! H2-SO4 + NH4NO3
                  aop_in(ivec)%BC (imode),  aop_in(ivec)%OC (imode), &
                  aop_in(ivec)%SOA (imode), &
                  aop_in(ivec)%SS (imode),  aop_in(ivec)%DU (imode), &
                  aop_in(ivec)%h2o(imode),  imode, m_eff, statusomp, gol)
             if (statusomp ==1) then
                ! Problem tracked to emilnox, should not happen anymore
                ! but keep the printout and hack for now

                !FIX to finish crescendo

                write (gol,'(" Problem with GET_REFR_IDX-NAN inputs:  ", 9(E16.4,2x))') &
                     aop_in(ivec)%SO4(imode), aop_in(ivec)%NO3(imode), aop_in(ivec)%h2o(imode), &
                     aop_in(ivec)%SOA(imode),aop_in(ivec)%BC(imode), aop_in(ivec)%OC(imode),&
                     aop_in(ivec)%SS(imode), aop_in(ivec)%DU(imode)
                call GoErr
                ! Assume edge case not caught. Set m_eff to default (1.0,1.0e-9), same as in ref_idx
                ! (subroutine problems with inputs SO4 and WATER?)
                status = 0
                m_eff= Cmplx(1.0,1.0e-9)
                ! If all edge cases were handled correctly, that should be:
                !status = 1
                !IF_ERROR_RETURN( status=1 )
             endif

             ! cmk added towpi for new netcdf lookup table
             xg = twopi*aop_in(ivec)%rg(imode) / wdep(i)%wl

             !=======================================================================
             !     get aerosol optical properties from data base for each mode
             !=======================================================================
             ! add extra safeguard for negative xg. It is unphysical, but have a safeguard anyway
             ! added for refr_idx_nan problem 

             if (xg .gt. 0.0) then
                call Optics_Get(m_eff, xg, Cexti, ai, gi, cext_table, a_table, g_table, statusomp, gol )
             else
                Cexti=0.0
                ai=1.0
                gi=1.0
             end if
             if (statusomp ==1) then
                call GoErr
                status = 1
                IF_ERROR_RETURN( status=1 )
             endif

             ! Multiply Cext with lambda^2 to get the cross section.
             Cexti = Cexti*(wdep(i)%wl**2)
             ! this here is extinction coefficient in this mode 
             incext(ivec) = aop_in(ivec)%numdens(imode) * Cexti
             ! sum up partial coefficients
             Aop_out_ext(ivec,i,1) = Aop_out_Ext(ivec,i,1) + incext(ivec)

             ! scattering portion
             NCscai = ai * incext(ivec)

             ! sum up weights for average (both albedo and asymmetry)
             NCsca    (ivec)   = NCsca    (ivec)   + NCscai 
             aop_out_g(ivec,i) = aop_out_g(ivec,i) + NCscai * gi

          end do
    
          ! Split extinction to separate contributions from constituents in modes.
          ! A volume mixing is assumed (in contrast to the explicit mixing in get_refr_ind). 
          if( wdep(i)%split .or. wdep(i)%insitu) then 
             do ivec = 1, lvec

                if (wdep(i)%split) then

                   ! The fine-mode contributions to the extinction
                   ! includes the contributions from particles
                   ! with (wet) diameters smaller than 1 micron.
                   ! For wet particles, only part of the accumulation mode
                   ! should be included, and the coarse mode should be excluded.
                   if (.not.coarse) then
                      ! Currently, the contribution of the accumulation mode
                      ! is approximated using weighting by volume scaling factor modfrac.
                      ! For extinction, area weighting would probably be more appropriate!!
                      if (imode .eq. 3 .or. imode .eq. 6 ) then
                         ! - convert number mean radius to volume mean radius (by cmedr2mmedr(imode))
                         ! - 1 micron diameter --> radius is 0.5 microns (rg is also in microns)
                         modfrac = 0.5 + 0.5 * ERF( log( 0.5 / (aop_in(ivec)%rg(imode) * cmedr2mmedr(imode)) ) / &
                              (sqrt(2.0) * lnsigma(imode)) )
                      else
                         ! Include full nucleation and Aitken mode contributions.
                         modfrac = 1.0
                      endif
                      aop_out_ext(ivec,i,10) = aop_out_ext(ivec,i,10) + modfrac * incext(ivec)
                   endif
                    ! total volume from so4/no3/bc/oc/ss/du in this mode (ATTENTION: DRY!!) 
                ! take no3 as so4
                totvoldry = aop_in(ivec)%so4(imode)/so4_density    + aop_in(ivec)%no3(imode)/so4_density + &
                            aop_in(ivec)%bc (imode)/carbon_density + aop_in(ivec)%oc (imode)/pom_density + &
                            aop_in(ivec)%soa (imode)/soa_density + &
                            aop_in(ivec)%ss (imode)/ss_density     + aop_in(ivec)%du (imode)/dust_density
                ! check whether there is some volume available 
                ! otherwise assign zeros to extinction increments
                if( totvoldry < tiny(totvoldry) ) then 
                   write(gol,'("WARNING: no volume in mode (",i3,"). assigning zero extinctions")') imode; call goPr
                   cycle
                end if
                ! H2-SO4 contribution
                aop_out_ext(ivec,i,2) = aop_out_ext(ivec,i,2) + incext(ivec) * (aop_in(ivec)%so4(imode)/so4_density   ) / totvoldry 
                ! BC contribution
                aop_out_ext(ivec,i,3) = aop_out_ext(ivec,i,3) + incext(ivec) * (aop_in(ivec)%bc (imode)/carbon_density) / totvoldry 
                ! POM contribution
                aop_out_ext(ivec,i,4) = aop_out_ext(ivec,i,4) + incext(ivec) * (aop_in(ivec)%oc (imode)/pom_density   ) / totvoldry 
                ! SOA contribution
                aop_out_ext(ivec,i,5) = aop_out_ext(ivec,i,5) + incext(ivec) * (aop_in(ivec)%soa (imode)/soa_density   ) / totvoldry 
                ! SS contribution
                aop_out_ext(ivec,i,6) = aop_out_ext(ivec,i,6) + incext(ivec) * (aop_in(ivec)%ss (imode)/ss_density    ) / totvoldry 
                ! DU contribution
                aop_out_ext(ivec,i,7) = aop_out_ext(ivec,i,7) + incext(ivec) * (aop_in(ivec)%du (imode)/dust_density  ) / totvoldry 
                ! NH4NO3 contribution
                aop_out_ext(ivec,i,8) = aop_out_ext(ivec,i,8) + incext(ivec) * (aop_in(ivec)%no3(imode)/so4_density   ) / totvoldry
                ! Fine-mode contributions for dust and sea salt 
                if (.not.coarse) then
                  aop_out_ext(ivec,i,11) = aop_out_ext(ivec,i,11) + incext(ivec) * (aop_in(ivec)%du (imode)/dust_density  ) / totvoldry * modfrac

                  aop_out_ext(ivec,i,12) = aop_out_ext(ivec,i,12) + incext(ivec) * (aop_in(ivec)%ss (imode)/ss_density  ) / totvoldry * modfrac
               endif
            endif

                   ! Water contribution: 
                   ! Get optical properties without water, the difference will be extinction due to 
                   ! water existence
                   ! - mis-use gi for this 
                   gi = 0.0 
                   call get_refr_idx( wdep(i), &
                        aop_in(ivec)%SO4(imode) + aop_in(ivec)%NO3(imode), &
                        aop_in(ivec)%BC (imode),  aop_in(ivec)%OC (imode), &
                        aop_in(ivec)%SOA (imode), &
                        aop_in(ivec)%SS (imode),  aop_in(ivec)%DU (imode), &
                        gi, imode, m_eff, statusomp, gol)
                   if (statusomp ==1) then
                       write (gol,'("Problem with GET_REFR_IDX-NAN inputs:  ", 9(E16.4,2x))') &
                            aop_in(ivec)%SO4(imode), aop_in(ivec)%NO3(imode), aop_in(ivec)%h2o(imode), &
                            aop_in(ivec)%SOA(imode),aop_in(ivec)%BC(imode), aop_in(ivec)%OC(imode),&
                            aop_in(ivec)%SS(imode), aop_in(ivec)%DU(imode)
                       call GoErr
                       ! Assume edge case not caught. Set m_eff to default (1.0,1.0e-9), same as in ref_idx
                       ! (subroutine problems with inputs SO4 and WATER?)
                       status = 0
                       m_eff= Cmplx(1.0,1.0e-9)
                       ! If all edge cases were handled correctly, that should be:
                       !status = 1
                       !IF_ERROR_RETURN( status=1 )
                       !status = 1
                   endif
                   
                   ! here we need the dry radius!!! 
                   !>>> TvN
                   ! 2*pi should be included (see comment above)
                   !xg = aop_in(ivec)%rgd(imode) / wdep(i)%wl
                   xg = twopi*aop_in(ivec)%rgd(imode) / wdep(i)%wl
                   !<<< TvN
                   ! TB
                   ! add extra safeguard for negative xg, unphysical but have a safeguard
                   ! added for refr_idx_nan problem 
                   if (xg .gt. 0.0) then
                      call Optics_Get(m_eff, xg, Cexti, ai, gi, cext_table, a_table, g_table, statusomp, gol )  
                   else
                      Cexti=0.0
                      ai=1.0
                      gi=1.0
                   end if
                   
                   if (statusomp ==1) then
                      call GoErr
                      status = 1
                      IF_ERROR_RETURN( status=1 )
                   endif
                   
                   Cexti = Cexti*(wdep(i)%wl**2)
                   
                   if (wdep(i)%split) then
                      ! add difference to water subarray
                      Aop_out_ext(ivec,i,9) = aop_out_ext(ivec,i,9) + (incext(ivec) - aop_in(ivec)%numdens(imode) * Cexti)
                   endif
                   
                   if (wdep(i)%insitu) then
                      ! Surface dry extinction and absorption: 
                      !>>> TvN
                      ! Remove cut off for the total values:
                      !modfrac = 0.5 + 0.5 * ERF( log( 5.0 / (aop_in(ivec)%rgd(imode) * cmedr2mmedr(imode)) ) / &
                      !     (sqrt(2.0) * lnsigma(imode)) )
                      ! extinction and absorption (extinction-scattering):
                      !Aop_out_add(ivec,i,1) = aop_out_add(ivec,i,1) + modfrac * aop_in(ivec)%numdens(imode) * Cexti
                      !Aop_out_add(ivec,i,2) = aop_out_add(ivec,i,2) + modfrac * aop_in(ivec)%numdens(imode) * Cexti * (1. - ai)
                      aop_out_add(ivec,i,1) = aop_out_add(ivec,i,1) + aop_in(ivec)%numdens(imode) * Cexti
                      aop_out_add(ivec,i,2) = aop_out_add(ivec,i,2) + aop_in(ivec)%numdens(imode) * Cexti * (1. - ai)
                      aop_out_add(ivec,i,3) = aop_out_add(ivec,i,3) + aop_in(ivec)%numdens(imode) * Cexti * ai * gi
                      ! Add fine-mode contributions
                      ! For dry aerosol, the fine-mode optical properties include 
                      ! the full accumulation mode, but not the coarse mode:
                      if (.not.coarse) then
                         aop_out_add(ivec,i,4) = aop_out_add(ivec,i,4) + aop_in(ivec)%numdens(imode) * Cexti
                         aop_out_add(ivec,i,5) = aop_out_add(ivec,i,5) + aop_in(ivec)%numdens(imode) * Cexti * (1. - ai)
                      endif
                      !<<< TvN
                   endif

             end do ! ivec
          end if

       enddo   ! modes

       if (ecearth_units) then
          ! return extinction due to absorption
          ! and asymmetry factor multiplied by extinction due to scattering
          ! and convert um**2/m**3 into 1/m (as for extinction below)
          Aop_out_a(:,i) = ( Aop_out_Ext(:,i,1) - NCsca(:) ) * 1.e-12
          Aop_out_g(:,i) = Aop_out_g(:,i) * 1.e-12
       else
          ! return single-scattering albedo and asymmetry factor
          ! take into account small extinction values
          where( Aop_out_Ext(:,i,1) > tiny(0.0) )
             Aop_out_a(:,i) = NCsca(:) / Aop_out_Ext(:,i,1)
          elsewhere
             ! No extinction -> Fill 1.0 for albedo because it looks clean.
             Aop_out_a(:,i) = 1.0
          endwhere

          ! take into account small extinction values
          where( NCsca(:) > tiny(0.0) )
             Aop_out_g(:,i) = Aop_out_g(:,i) / NCsca(:)
          elsewhere
             ! No scattering at all -> Fill 1.0 for asymmetry because it looks clean.
             Aop_out_g(:,i) = 1.0
          endwhere
       endif

       ! convert um**2/m**3 into 1/m
       Aop_out_Ext(:,i,:) = Aop_out_Ext(:,i,:) * 1e-12

       if( wdep(i)%insitu) then 
          ! take into account small extinction values
          where( aop_out_add(:,i,1) - aop_out_add(:,i,2) > tiny(0.0) )
            aop_out_add(:,i,3) = aop_out_add(:,i,3) / (aop_out_add(:,i,1)-aop_out_add(:,i,2))
          elsewhere
          ! No scattering at all -> Fill 1.0 for asymmetry because it looks clean.
            aop_out_add(:,i,3) = 1.0
          endwhere

          Aop_out_add(:,i,1) = aop_out_add(:,i,1) * 1e-12 
          Aop_out_add(:,i,2) = aop_out_add(:,i,2) * 1e-12 
          Aop_out_add(:,i,4) = aop_out_add(:,i,4) * 1e-12 
          Aop_out_add(:,i,5) = aop_out_add(:,i,5) * 1e-12
       endif

       if( wdep(i)%split .or. wdep(i)%insitu ) then
          deallocate( lnsigma )
       endif

       !=======================================================================
    enddo   ! loop over wavelengths

    

    Nullify(Cext_table)
    Nullify(a_table)
    Nullify(g_table)

    !$OMP END PARALLEL

    if (status==1) then 
       call goPr
       IF_NOTOK_RETURN(status=1)
    endif

    !do imode = 1,7
    !      print *, 'Radius,mode   :', imode, sum(aop_in(:)%rg(imode))/lvec
    !      print *, 'numden,mode   :', imode, sum(aop_in(:)%numdens(imode))/lvec
    !enddo

    deallocate( NCsca, incext )

    !do i = 1,nwl
    !   print *, 'AOD per grid box:', wdep(i)%wl, sum(Aop_out_Ext(:,i,1))/lvec
    !enddo

    ! ok
    status = 0

  END SUBROUTINE OPTICS_CALCULATE_AOP
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OPTICS_AOP_GET
  !
  ! !DESCRIPTION: Initialise the fields "aop_in" and then calculate the
  !               optical properties through a call to optics_calculate_aop. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE OPTICS_AOP_GET( lvec, region, nwav, wdep, ncontr, ecearth_units, &
                             aop_out_ext, aop_out_a, aop_out_g, status, aop_out_add )
    !
    ! !USES:
    !
    USE TM5_DISTGRID,   ONLY : dgrid, Get_DistGrid
    use Meteo,          only : Set
    Use Meteodata,      only : gph_dat, m_dat
    Use Dims,           only : im, jm, lm
    use global_data,    only : region_dat, mass_dat
    Use chem_param,     only : inus_n, iso4nus, isoanus, nmod
    use chem_param,     only : iais_n, iso4ais, ibcais, ipomais, isoaais
    Use chem_param,     only : iacs_n, iso4acs, ibcacs, ipomacs, isoaacs, issacs, iduacs, ino3_a
    use chem_param,     only : imsa
    use chem_param,     only : icos_n, iso4cos, ibccos, ipomcos, isoacos, isscos, iducos
    use chem_param,     only : iaii_n,          ibcaii, ipomaii, isoaaii
    use chem_param,     only : iaci_n, icoi_n,                           iduaci, iducoi
    use chem_param,     only : h2so4_factor, nh4no3_factor
    use chem_param,     only : so4_density, nh4no3_density, msa_density
    Use m7_data,        only : h2o_mode, rw_mode, rwd_mode
    !
    ! !INPUT PARAMETERS:
    !
    Integer,                           Intent(In) :: region, lvec
    Integer,                           Intent(In) :: nwav, ncontr
    type(wavelendep), dimension(nwav), Intent(In) :: wdep
    logical,                           intent(in) :: ecearth_units
    !
    ! !OUTPUT PARAMETERS:
    !
    real, dimension(lvec,nwav,ncontr), intent(out)  :: aop_out_ext ! extinctions
    real, dimension(lvec,nwav),        intent(out)  :: aop_out_a   ! single scattering albedo
    real, dimension(lvec,nwav),        intent(out)  :: aop_out_g   ! assymetry factor
    integer,                           intent(out)  :: status
    real, dimension(:,:,:),            intent(out), optional :: aop_out_add ! additional parameters 
    !
    ! !REVISION HISTORY: 
    !   29 Nov 2010 - Achim Strunk - v0
    !   24 Jun 2011 - Achim Strunk - Added NO3 to allow for explicit 
    !                                AOD split of (SO4+NO3)
    !   20 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter ::  rname = mname//'/optics_aop_get'

    real, dimension(:,:,:,:), pointer     :: rm
    real, dimension(:,:,:),   pointer     :: m
    real, dimension(:,:,:),   allocatable :: volume
    Integer                               :: i, j, l, iloop, i1, i2, j1, j2, lmr


    ! --- start ------------------------------

    ! ensure met fields being available
    call Set( gph_dat(region), status, used=.true. )

    ! air & tracers mass
    m  => m_dat(region)%data
    rm => mass_dat(region)%rm

    ! tile grid size
    CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

    ! Generalize to allow for reduced number of levels 
    ! when called from ECEarth_Optics_Step
    !lmr=lm(region)
    lmr = lvec / ( (i2-i1+1)*(j2-j1+1) )

    allocate( aop_in(lvec) )

    do i = 1, nmod
       aop_in(:)%so4    (i) = 0.0 ; aop_in(:)%bc   (i) = 0.0
       aop_in(:)%oc     (i) = 0.0 ; aop_in(:)%soa  (i) = 0.0
       aop_in(:)%ss     (i) = 0.0
       aop_in(:)%du     (i) = 0.0 ; aop_in(:)%h2o  (i) = 0.0
       aop_in(:)%numdens(i) = 0.0 ; aop_in(:)%rg   (i) = 0.0
       aop_in(:)%rgd    (i) = 0.0 ; aop_in(:)%no3  (i) = 0.0
    end do

!>>> TvN
! In M7 sulphate is assumed to be H2-SO4 with corresponding particle density so4_density
! The sulphate mass should therefore also include the small contribution from the H atoms
    ! NUS
    aop_in(:)%so4(1) = reshape( 1.e9 * h2so4_factor * rm(i1:i2,j1:j2,1:lmr,iso4nus) / &
                                                       m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%soa(1) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,isoanus) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )

    ! AIS 
    aop_in(:)%so4(2) = reshape( 1.e9 * h2so4_factor * rm(i1:i2,j1:j2,1:lmr,iso4ais) / &
                                                       m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%bc (2) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ibcais ) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%oc (2) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ipomais) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%soa(2) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,isoaais) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    ! ACS (additional: NO3)
! The contribution from methane sulfonate (MSA-) aerosol is added
! to that for sulfate.
! As the addition is done by volume,
! we need to account for the difference in densities
! (as done below for ammonium nitrate).
    aop_in(:)%so4(3) = reshape( 1.e9 * ( h2so4_factor * rm(i1:i2,j1:j2,1:lmr,iso4acs) + &
                          (so4_density / msa_density) * rm(i1:i2,j1:j2,1:lmr,imsa) ) /  &
                                                         m(i1:i2,j1:j2,1:lmr), (/lvec/) )
! Since nh4no3_density is the density of NH4NO3, the contribution from NH4 should be included.
! Moreover, assuming the same refractive index for NH4NO3 as for H2-SO4,
! the contributions from both components can be added by volume;
! thus we need to account for the difference in densities.
! Estimates of the refractive index of NH4NO3 are available from literature
! (e.g. Lowenthal et al., Atmos. Environ., 2000).
! For practical purposes, it can be set equal to the value used for sulfate,
! i.e. the value for a solution containing 75% H2SO4 (Fenn et al., 1985).
    aop_in(:)%no3(3) = reshape( 1.e9 * nh4no3_factor * (so4_density / nh4no3_density) * &
                                       rm(i1:i2,j1:j2,1:lmr,ino3_a ) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%bc (3) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ibcacs ) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%oc (3) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ipomacs) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%soa(3) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,isoaacs) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%ss (3) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,issacs ) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%du (3) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,iduacs ) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    ! COS
    aop_in(:)%so4(4) = reshape( 1.e9 * h2so4_factor * rm(i1:i2,j1:j2,1:lmr,iso4cos) / &
                                                       m(i1:i2,j1:j2,1:lmr), (/lvec/) )
!<<< TvN
    aop_in(:)%bc (4) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ibccos ) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%oc (4) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ipomcos) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%soa(4) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,isoacos) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%ss (4) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,isscos ) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%du (4) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,iducos ) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    ! AII
    aop_in(:)%bc (5) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ibcaii ) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%oc (5) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ipomaii) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%soa(5) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,isoaaii) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    ! ACI
    aop_in(:)%du (6) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,iduaci ) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    ! COI
    aop_in(:)%du (7) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,iducoi ) / &
                                        m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    ! Water in (hydrophillic) modes
    aop_in(:)%h2o(1) = reshape( 1.e9 * h2o_mode(region,1)%d3(i1:i2,j1:j2,1:lmr) / &
                                                           m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%h2o(2) = reshape( 1.e9 * h2o_mode(region,2)%d3(i1:i2,j1:j2,1:lmr) / &
                                                           m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%h2o(3) = reshape( 1.e9 * h2o_mode(region,3)%d3(i1:i2,j1:j2,1:lmr) / &
                                                           m(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%h2o(4) = reshape( 1.e9 * h2o_mode(region,4)%d3(i1:i2,j1:j2,1:lmr) / &
                                                           m(i1:i2,j1:j2,1:lmr), (/lvec/) )

    aop_in(:)%rg (1) = reshape( 1.e6 * rw_mode (region,1)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%rg (2) = reshape( 1.e6 * rw_mode (region,2)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%rg (3) = reshape( 1.e6 * rw_mode (region,3)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%rg (4) = reshape( 1.e6 * rw_mode (region,4)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%rg (5) = reshape( 1.e6 * rw_mode (region,5)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%rg (6) = reshape( 1.e6 * rw_mode (region,6)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%rg (7) = reshape( 1.e6 * rw_mode (region,7)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) )

    ! dry radius for soluble modes / rest equals the usual radii
    aop_in(:)%rgd(1) = reshape( 1.e6 * rwd_mode(region,1)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%rgd(2) = reshape( 1.e6 * rwd_mode(region,2)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%rgd(3) = reshape( 1.e6 * rwd_mode(region,3)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%rgd(4) = reshape( 1.e6 * rwd_mode(region,4)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) )
    aop_in(:)%rgd(5) = aop_in(:)%rg (5)
    aop_in(:)%rgd(6) = aop_in(:)%rg (6)
    aop_in(:)%rgd(7) = aop_in(:)%rg (7)

    ! volume for number density
    allocate( volume( i1:i2, j1:j2, 1:lmr ) )
    do j = j1, j2
          volume(:,j,:) = ( gph_dat(region)%data(i1:i2, j, 2:lmr+1) -   &
                            gph_dat(region)%data(i1:i2, j, 1:lmr  ) ) * &
                            region_dat(region)%dxyp(j) ! m3
    end do
       
    aop_in(:)%numdens(1) = reshape( rm(i1:i2,j1:j2,1:lmr,inus_n) / volume, (/lvec/) )
    aop_in(:)%numdens(2) = reshape( rm(i1:i2,j1:j2,1:lmr,iais_n) / volume, (/lvec/) )
    aop_in(:)%numdens(3) = reshape( rm(i1:i2,j1:j2,1:lmr,iacs_n) / volume, (/lvec/) )
    aop_in(:)%numdens(4) = reshape( rm(i1:i2,j1:j2,1:lmr,icos_n) / volume, (/lvec/) )
    aop_in(:)%numdens(5) = reshape( rm(i1:i2,j1:j2,1:lmr,iaii_n) / volume, (/lvec/) )
    aop_in(:)%numdens(6) = reshape( rm(i1:i2,j1:j2,1:lmr,iaci_n) / volume, (/lvec/) )
    aop_in(:)%numdens(7) = reshape( rm(i1:i2,j1:j2,1:lmr,icoi_n) / volume, (/lvec/) )
    deallocate( volume )

    ! check valid ranges in particle sizes (might be zero)
    where( aop_in%rg (1) .lt. 1.E-15 ) aop_in%rg (1) = 1.E-15
    where( aop_in%rgd(1) .lt. 1.E-15 ) aop_in%rgd(1) = 1.E-15
    where( aop_in%rg (2) .lt. 1.E-15 ) aop_in%rg (2) = 1.E-15
    where( aop_in%rgd(2) .lt. 1.E-15 ) aop_in%rgd(2) = 1.E-15
    where( aop_in%rg (3) .lt. 1.E-15 ) aop_in%rg (3) = 1.E-15
    where( aop_in%rgd(3) .lt. 1.E-15 ) aop_in%rgd(3) = 1.E-15
    where( aop_in%rg (4) .lt. 1.E-15 ) aop_in%rg (4) = 1.E-15
    where( aop_in%rgd(4) .lt. 1.E-15 ) aop_in%rgd(4) = 1.E-15
    where( aop_in%rg (5) .lt. 1.E-15 ) aop_in%rg (5) = 1.E-15
    where( aop_in%rgd(5) .lt. 1.E-15 ) aop_in%rgd(5) = 1.E-15
    where( aop_in%rg (6) .lt. 1.E-15 ) aop_in%rg (6) = 1.E-15
    where( aop_in%rgd(6) .lt. 1.E-15 ) aop_in%rgd(6) = 1.E-15
    where( aop_in%rg (7) .lt. 1.E-15 ) aop_in%rg (7) = 1.E-15
    where( aop_in%rgd(7) .lt. 1.E-15 ) aop_in%rgd(7) = 1.E-15

    nullify(rm)
    nullify(m)

    ! 
    aop_out_ext=0.0
    aop_out_a=0.0
    aop_out_g=0.0
    if (present(aop_out_add)) then
       call optics_calculate_aop( nwav, wdep, lvec, ecearth_units, &
                                  aop_out_ext, aop_out_a, aop_out_g, status, aop_out_add )
    else
       call optics_calculate_aop( nwav, wdep, lvec, ecearth_units, &
                                  aop_out_ext, aop_out_a, aop_out_g, status )
    endif
    IF_NOTOK_RETURN(status=1)

    ! OK
    Deallocate( aop_in )
    status = 0

  END SUBROUTINE OPTICS_AOP_GET
  !EOC

END MODULE OPTICS