#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') 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:      ECEARTH_OPTICS
!
! !DESCRIPTION: Optics module to calculate optical depth from m7 output, based
!               on the AOP_Package of Michael Kahnert.
!\\
!\\
! !INTERFACE: 
!
MODULE ECEARTH_OPTICS
  !
  ! !USES:
  !
  use GO,           only : gol, goErr, goPr
  use global_types, only : d3_data
  use optics,       only : wavelendep   ! wavelength type
  use tm5_distgrid, only : dgrid, Get_DistGrid

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: ECEarth_Optics_Init, ECEarth_Optics_Done, ECEarth_Optics_Step
  !
  ! !PUBLIC TYPES:
  !
  type, public :: optics_data
     real, dimension(:,:,:,:), pointer :: Ext, a, g
  end type optics_data
  !
  ! !PUBLIC DATA MEMBERS:
  !
  integer,           public, parameter :: nregions_optics = 1  ! number of regions for which optical properties are calculated for EC-Earth
  type(optics_data), public, dimension(nregions_optics), target :: optics_dat  ! optical properties container
  !
  ! !PRIVATE DATA MEMBERS:
  !
  !integer,               parameter :: nwl = 6   ! number of wavelengths for EC-Earth v2
  integer,               parameter :: nwl = 14  ! number of wavelengths for McRad (EC-Earth v3 and later)
  type(wavelendep), dimension(nwl) :: ecearth_wvl
  !
  character(len=*), parameter  ::  mname = 'ECEarth_Optics'
  !
  ! !REVISION HISTORY: 
  !      Mar 2009 - Maarten Krol   - Implemented
  !      Apr 2009 - Twan van Noije - Adapted
  !    5 Sep 2013 - Ph. Le Sager   - adapted for TM5v4
  !   19 Sep 2013 - Ph. Le Sager   - adapted for lon-Lat MPI decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    ECEarth_Optics_Init
  !
  ! !DESCRIPTION: Read lookup table, initialize wavelengths and their
  !               parameters. Allocate work arrays. Called from initexit/start. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE ECEarth_Optics_Init( status )
    !
    ! !USES:
    !
    use dims,         only : lm
    use optics,       only : optics_init
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !    5 Sep 2013 - Ph. Le Sager - TM5v4 version
    !   19 Sep 2013 - Ph. Le Sager - adapted for lon-Lat MPI decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer       :: region, i1, i2, j1, j2, lmr

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

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

    write (gol,'("Initializing M7 aerosol optical properties for EC-Earth (ancillary fields) ...")'); call goPr 

    ! ------------------------------------------
    ! ecearth optics wavelengths
    ! for 6-band SW radiative code in IFS
    ! 0.185-0.25-0.44-0.69-1.19-2.38-4.00 micron
    ! ------------------------------------------     
    ! ecearth_wvl( 1)%wl = 0.2175
    ! ecearth_wvl( 2)%wl = 0.345
    ! ecearth_wvl( 3)%wl = 0.565
    ! ecearth_wvl( 4)%wl = 0.94
    ! ecearth_wvl( 5)%wl = 1.785
    ! ecearth_wvl( 6)%wl = 3.19

    ! -------------------------------------
    ! ecearth optics wavelengths
    ! for "McRad" radiative code in IFS,
    ! adapted to account for solar spectrum
    ! as done in IFS, e.g. for MACv2-SP.
    ! -------------------------------------
    !ecearth_wvl( 1)%wl = 0.2316
    !ecearth_wvl( 2)%wl = 0.3040
    !ecearth_wvl( 3)%wl = 0.3932
    !ecearth_wvl( 4)%wl = 0.5332
    !ecearth_wvl( 5)%wl = 0.7016
    !ecearth_wvl( 6)%wl = 1.0101
    !ecearth_wvl( 7)%wl = 1.2705
    !ecearth_wvl( 8)%wl = 1.4625
    !ecearth_wvl( 9)%wl = 1.7840
    !ecearth_wvl(10)%wl = 2.0460
    !ecearth_wvl(11)%wl = 2.3250
    !ecearth_wvl(12)%wl = 2.7885
    !ecearth_wvl(13)%wl = 3.4615
    !ecearth_wvl(14)%wl = 8.0205
    ecearth_wvl( 1)%wl = 0.257 
    ecearth_wvl( 2)%wl = 0.313
    ecearth_wvl( 3)%wl = 0.398
    ecearth_wvl( 4)%wl = 0.530
    ecearth_wvl( 5)%wl = 0.697
    ecearth_wvl( 6)%wl = 0.973
    ecearth_wvl( 7)%wl = 1.269
    ecearth_wvl( 8)%wl = 1.447
    ecearth_wvl( 9)%wl = 1.767
    ecearth_wvl(10)%wl = 2.040
    ecearth_wvl(11)%wl = 2.308
    ecearth_wvl(12)%wl = 2.752
    ecearth_wvl(13)%wl = 3.407
    ecearth_wvl(14)%wl = 5.254

    ecearth_wvl%split = .false. 
    ecearth_wvl%insitu = .false.

    ! ------------------------------------------
    ! initialise the optics
    ! ------------------------------------------
    ! read lookup table, initialize wavelength dependent parameters
    call optics_init( nwl, ecearth_wvl, status )
    IF_NOTOK_RETURN(status=1)

    ! ------------------------------------------    
    ! allocate optical properties arrays.
    ! ------------------------------------------    
    do region = 1, nregions_optics

       call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2,  J_STRT=j1, J_STOP=j2 )
       lmr = lm(region)

       allocate ( optics_dat(region)%Ext(i1:i2,j1:j2,lmr,nwl)); optics_dat(region)%Ext = 0.
       allocate ( optics_dat(region)%a  (i1:i2,j1:j2,lmr,nwl)); optics_dat(region)%a = 0.
       allocate ( optics_dat(region)%g  (i1:i2,j1:j2,lmr,nwl)); optics_dat(region)%g = 0.

    end do

    ! ok
    status = 0

  END SUBROUTINE ECEarth_Optics_Init
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    ECEARTH_OPTICS_STEP
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine ECEarth_Optics_Step( status )
    !
    ! !USES:
    !
    use optics,      only : optics_aop_get
    use dims,        only : lm
    use tm5_prism,   only : ifs_cpl_nlev, ifs_cpl_nlev_cutoff, refine_levels
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)  ::  status
    !
    ! !REVISION HISTORY: 
    !    5 Sep 2013 - Ph. Le Sager - TM5v4
    !   19 Sep 2013 - Ph. Le Sager - adapted for lon-Lat MPI decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer                             :: region, i1, i2, j1, j2, lmr
    integer                             :: lwl, lvec
    real, dimension(:,:,:), allocatable :: aop_out_ext
    real, dimension(:,:),   allocatable :: aop_out_a, aop_out_g
    logical, parameter :: new_units=.false.

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

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

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

    REG: do region = 1, nregions_optics

       call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2,  J_STRT=j1, J_STOP=j2 )

       if (ifs_cpl_nlev_cutoff /= ifs_cpl_nlev) then
          if (.NOT. refine_levels) then
             ! Limit optics calculations at IFS wavelengths to ifs_cpl_nlev_cutoff.
             ! The values at higher levels are not sent to IFS
             ! and can therefore be left zero.
             lmr = ifs_cpl_nlev_cutoff
          else
             write(gol,'("WARNING: EC-Earth optics called for the whole atmosphere,")'); call goPr 
             write(gol,'("WARNING: but not all levels are sent to IFS")'); call goPr
             lmr = lm(region)
          endif
       else
          lmr = lm(region)
       endif

       lvec = (i2-i1+1)*(j2-j1+1)*lmr

       ! allocate AOP fields
       allocate( aop_out_ext(lvec, nwl, 1 ) ) ! extinction
       allocate( aop_out_a  (lvec, nwl    ) ) ! single scattering albedo (extinction due to absorption if new_units)
       allocate( aop_out_g  (lvec, nwl    ) ) ! asymmetry factor (times extinction due to scattering if new_units) 

       ! Compute AOP (ie fill aop_out_* arrays)
       call optics_aop_get( lvec, region, nwl, ecearth_wvl, 1, new_units, &
                            aop_out_ext, aop_out_a, aop_out_g, status )
       IF_NOTOK_RETURN(status=1)

       ! ---------------------------------
       ! unpack results from aop computation   
       ! ---------------------------------
       do lwl = 1, nwl
          optics_dat(region)%ext(:,:,1:lmr,lwl) = reshape( aop_out_ext(:,lwl,1), (/(i2-i1+1),(j2-j1+1),lmr/) )
          optics_dat(region)%a  (:,:,1:lmr,lwl) = reshape( aop_out_a  (:,lwl),   (/(i2-i1+1),(j2-j1+1),lmr/) ) 
          optics_dat(region)%g  (:,:,1:lmr,lwl) = reshape( aop_out_g  (:,lwl),   (/(i2-i1+1),(j2-j1+1),lmr/) ) 
       end do

       ! free temporary arrays for results from calculate_aop
       deallocate( aop_out_ext, aop_out_a, aop_out_g )

    end do REG

    status = 0

  END SUBROUTINE ECEARTH_OPTICS_STEP
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    ECEARTH_OPTICS_DONE
  !
  ! !DESCRIPTION: Deallocate optical properties arrays. Called from initexit/exitus.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE ECEarth_Optics_Done( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !    5 Sep 2013 - Ph. Le Sager -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer :: region

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

    do region = 1, nregions_optics
       deallocate ( optics_dat(region)%Ext )
       deallocate ( optics_dat(region)%a   )
       deallocate ( optics_dat(region)%g   )
    enddo

    ! ok
    status = 0

  END SUBROUTINE ECEarth_Optics_Done
  !EOC

END MODULE ECEARTH_OPTICS