!>
!! @brief Module to provide the idealized aerosol fields for the Easy Aerosol project
!!
!! @remarks
!!   This module contains routines that allow to set the shape of the
!!   idealized aerosol field via the namelist and that, based on this, calculate the
!!   aerosol optical properties in the format needed by the radiation routine.
!!
!! @author Aiko Voigt, MPI-M, Hamburg (2013-04-09):
!!         Aiko Voigt, LMD, Paris (revised in July 2013, August 2013)
!!         Aiko Voigt, LDEO, Columbia University, New York City (revised in Dec 2013)
!!
!! $ID: n/a$
!!
!! @par Origin
!!   The code was based on mo_aero_kinne.f90.
!
MODULE mo_aero_easyaerosols
  
  USE mo_kind,      ONLY: dp
  USE mo_control,   ONLY: ldebugs
  USE mo_debugs

  IMPLICIT NONE
  PRIVATE

  PUBLIC    :: init_aero_easyaerosols, set_aop_easyaerosols

! default parameters
!
! note that the naming of parameters does not for all parameters follow the naming
! in the Easy Aerosol description. However, this should cause no confusion as the parameters
! are explained in this module (see below). 
!
! attention: the values given below are valid for the zonal-mean
! case (l3plumes=.false.)
! if plumes are specified, make sure to use the following values instead:
!  REAL(dp)  :: lon1     = 20.0_dp
!  REAL(dp)  :: lon2     = 110.0_dp
!  REAL(dp)  :: lon3     = 280.0_dp
!  REAL(dp)  :: siglon1  = 30.0_dp
!  REAL(dp)  :: siglon2  = 24.0_dp
!  REAL(dp)  :: siglon3  = 25.0_dp
!  REAL(dp)  :: fact1    = 1.8_dp
!  REAL(dp)  :: fact2    = 2.0_dp
!  REAL(dp)  :: fact3    = 0.9_dp
!  REAL(dp)  :: aod0     = 0.0932_dp
!
! note that the values of lon1,lon2,lon3,siglon1,siglon2,siglon3,fact1,fact2,fact3 are
! irrelevant when l3plumes=.false., the values are those appropriate for l3plumes=.true.

  LOGICAL   :: l3plumes = .false.        !! false: zonally-uniform aerosol, true: 3 plumes
  REAL(dp)  :: lat0     = 35.0_dp        !! center latitude of meridional Gaussian
  REAL(dp)  :: siglat   = 25.0_dp        !! width of meridonal Gaussian
  REAL(dp)  :: lon1     = 20.0_dp        !! center longitude of Europe plume 
  REAL(dp)  :: lon2     = 110.0_dp       !! center longitude of SE Asia plume
  REAL(dp)  :: lon3     = 280.0_dp       !! center longitude of US plume
  REAL(dp)  :: siglon1  = 30.0_dp        !! zonal width of Europe plume
  REAL(dp)  :: siglon2  = 24.0_dp        !! zonal width of SE Asia plume
  REAL(dp)  :: siglon3  = 25.0_dp        !! zonal width of US plume
  REAL(dp)  :: fact1    = 1.8_dp         !! relative magnitude of Europe plume
  REAL(dp)  :: fact2    = 2.0_dp         !! relative magnitude of SE Asia plume
  REAL(dp)  :: fact3    = 0.9_dp         !! relative magnitude of US plume
  REAL(dp)  :: aod0     = 0.08_dp        !! total aod at 550 nm
  REAL(dp)  :: a        = 0.98657e-3_dp  !! factor for vertical profile of beta
  REAL(dp)  :: c        = 0.75e-3_dp     !! exponential decay constant 
  REAL(dp)  :: hexp     = 1250.0_dp      !! height up to which beta is constant
  REAL(dp)  :: ang      = 1.8_dp         !! Angstroem coefficient
  REAL(dp)  :: ssa0     = 0.926_dp       !! single scattering albedo
  REAL(dp)  :: asy0     = 0.65_dp        !! asymmetry factor
  
CONTAINS
  !----------------------------------------------------------------
  !>
  !! @brief Reads namelist easyaersolsctl
  !!
  !! @remarks This routine is called from setup_radiation when the
  !! aerosol flag iaero is set to 8 in radctl (iaero=8).
  !
  SUBROUTINE init_aero_easyaerosols

    USE mo_namelist,       ONLY: position_nml, positioned, missing, &
                                 open_nml
    USE mo_mpi,            ONLY: p_parallel_io, p_io, p_bcast
    USE mo_exception,      ONLY: finish, message
    USE mo_submodel,       ONLY: print_value

    IMPLICIT NONE

    INTEGER :: ierr, inml, iunit

    NAMELIST /easyaerosolsctl/ l3plumes,lat0, siglat, lon1, siglon1, lon2, siglon2, lon3, siglon3, &
                               fact1, fact2, fact3, aod0, a, c, hexp, ang, ssa0, asy0

    CALL message('','Initializing mo_aero_easyaerosols')

    IF (p_parallel_io) THEN
      inml = open_nml ('namelist.echam')                 
      iunit= position_nml('easyaerosolsctl', inml, status=ierr)
      SELECT CASE (ierr)
          CASE (positioned)
             READ (iunit, easyaerosolsctl)
             CALL message ('','EASYAEROSOLSCTL namelist read')
          CASE (missing)
             CALL finish ('', 'Namelist for easy aerosols is missing')
      END SELECT
    END IF

    ! make sure that lon1,lon2,lon3 are between 0 and 360 deg
    ! make sure a, c, hexp, and ang are positive
    IF (p_parallel_io) THEN
       IF ((lon1.lt.0.0_dp)   .or. (lon2.lt.0.0_dp)   .or. (lon3.lt.0.0_dp)    .or. &
           (lon1.ge.360.0_dp) .or. (lon2.ge.360.0_dp) .or. (lon3.ge.360.0_dp)) THEN
          CALL finish ('', ' init_easy_aerosols: lon1,lon2,lon3 must be >= 0 and < 360 deg!') 
       ENDIF
       IF (a.lt.0.0_dp)    CALL finish ('', ' init_easy_aerosols: a must be positive!') 
       IF (c.lt.0.0_dp)    CALL finish ('', ' init_easy_aerosols: c must be positive!') 
       IF (hexp.lt.0.0_dp) CALL finish ('', ' init_easy_aerosols: hexp must be positive!') 
       IF (ang.lt.0.0_dp)  CALL finish ('', ' init_easy_aerosols: ang must be positive!') 
    ENDIF    

    CALL p_bcast(l3plumes, p_io)
    CALL p_bcast(lat0, p_io)
    CALL p_bcast(siglat, p_io)
    CALL p_bcast(lon1, p_io)
    CALL p_bcast(siglon1, p_io)
    CALL p_bcast(lon2, p_io)
    CALL p_bcast(siglon2, p_io)
    CALL p_bcast(lon3, p_io)
    CALL p_bcast(siglon3, p_io)
    CALL p_bcast(fact1, p_io)
    CALL p_bcast(fact2, p_io)
    CALL p_bcast(fact3, p_io)
    CALL p_bcast(aod0, p_io)
    CALL p_bcast(a, p_io)
    CALL p_bcast(c, p_io)
    CALL p_bcast(hexp, p_io)
    CALL p_bcast(ang, p_io)
    CALL p_bcast(ssa0, p_io)
    CALL p_bcast(asy0, p_io)

    CALL message('','----------------------------')
    CALL print_value('l3plumes',l3plumes)
    CALL print_value('lat0',lat0)
    CALL print_value('siglat',siglat)
    CALL print_value('lon1',lon1)
    CALL print_value('siglon1',siglon1)
    CALL print_value('lon2',lon2)
    CALL print_value('siglon2',siglon2)
    CALL print_value('lon3',lon3)
    CALL print_value('siglon3',siglon3)
    CALL print_value('fact1',fact1)
    CALL print_value('fact2',fact2)
    CALL print_value('fact3',fact3)
    CALL print_value('aod0',aod0)
    CALL print_value('a',a)
    CALL print_value('c',c)
    CALL print_value('hexp',hexp)
    CALL print_value('ang',ang)
    CALL print_value('ssa0',ssa0)
    CALL print_value('asy0',asy0)
    CALL message('','----------------------------')

  END SUBROUTINE init_aero_easyaerosols
  !-----------------------------------------------------------------------
  !>
  !! @brief Sets the aerosol optical properties for Easy Aerosols
  !!
  !! @remarks This routine sets the aerosol optical properties for
  !! all spectral bands (solar and infrared) at each grid point
  !! and vertical level using the parameters
  !! specified through the namelist easyaerosolsctl. 
  !! The routine is called from rrtm_interace to which it passes the arrays
  !! paer_XXX_Xw_vr. The latter contain the aerosol optical properties needed 
  !! by the radiation code. Because we do not consider longwave aerosol effects,
  !! the longwave array paer_tau_lw_vr is set to zero. The actual calculation
  !! of the aerosol fields is done by calling the routine 
  !! make_aop_easyaerosols, so set_aop_easyaerosols acts as the interface between
  !! the model and the definition of the Easy Aerosol fields.
  !! Note: the arrays paer_XXX_Xw_vr are ordered opposite to the usual vertical 
  !! ordering of ECHAM6. They start from the surface and go to the top-of-atmosphere.
  !! while the standard ordering in ECHAM6 is from the top-of-atmosphere to the
  !! surface. The name extension "_vr" indicates this difference in the
  !! vertical ordering.
  !
  SUBROUTINE set_aop_easyaerosols ( &
          & kproma,                 kbdim,              klev,             &
          & krow,                   nb_lw,              nb_sw,            &
          & paer_tau_lw_vr,         paer_tau_sw_vr,     paer_piz_sw_vr,   &
          & paer_cg_sw_vr)
  
    USE mo_scan_buffer,         ONLY: alnpr                    ! for call to geopot 
    USE mo_physical_constants,  ONLY: vtmpc1                   ! for virtual temperature
    USE mo_memory_g1a,          ONLY: tm1, qm1, xlm1, xim1     ! for virtual temperature  
    USE mo_physical_constants,  ONLY: rgrav                    ! inverse of gravitational acceleration
    USE mo_srtm_config,         ONLY: wavenum1, wavenum2       ! indizes 16:29 contain wavenumbers for shortwave bands
                                                               ! the corresponding wavelength in nm is given by 1e7/wavenumber
    USE mo_geoloc,              ONLY: philon_2d, philat_2d     ! geographical longitude and latitude in degrees

    IMPLICIT NONE

    !INPUT PARAMETERS
    INTEGER,INTENT(in)  :: kproma, &! actual block length
                           kbdim,  &! maximum block length
                           krow,   &! block index
                           klev,   &! number of vertical levels
                           nb_lw,  &! number of longwave bands  (far infrared)
                           nb_sw    ! number of shortwave bands (visible + near infrared)

    !OUTPUT PARAMETERS
    REAL(dp),INTENT(out)  :: paer_tau_lw_vr(kbdim,klev,nb_lw),    & ! aerosol optical depth for longwave 
                             paer_tau_sw_vr(kbdim,klev,nb_sw),    & ! aerosol optical depth for shortwave
                             paer_piz_sw_vr(kbdim,klev,nb_sw),    & ! single scattering albedo 
                             paer_cg_sw_vr(kbdim,klev,nb_sw)        ! asymmetry factor 

    !LOCAL VARIABLES
    REAL(dp), DIMENSION(kbdim,klev)       :: ztv, zalpha
    REAL(dp), DIMENSION(kbdim,klev+1)     :: geoph
    REAL(dp), DIMENSION(kbdim)            :: zref   
    REAL(dp), DIMENSION(kbdim,klev,nb_sw) :: aod_vr,ssa_vr,asy_vr
    REAL(dp), DIMENSION(kbdim,klev)       :: heightt,              & ! height of level top above ground in m
                                             heightb,              & ! height of level bottom above ground in m
                                             dheight,              & ! layer thickness in m
                                             heightt_vr,           & ! reversed vertical order
                                             heightb_vr,           & ! reversed vertical order
                                             dheight_vr              ! reversed vertical order
    REAL(dp), DIMENSION(nb_sw)            :: lambda1,              & ! upper wavelength bound of spectral band (in nm)
                                             lambda2,              & ! lower wavelenght bound of spectral band (in nm)
                                             lambda                  ! central wavelenght of spectral band (in nm)
    REAL(dp), DIMENSION(kbdim)            :: lon,lat                 ! longitude and latitude of grid box in degrees 
    INTEGER                               :: jk, klevp1

    klevp1 = klev + 1

    ! compute geopotential height of half levels to get height of bottom and top of level
    ! virtual temperature
    ztv(1:kproma,:) = tm1(1:kproma,:,krow)*(1._dp+vtmpc1*qm1(1:kproma,:,krow)  &
                      -(xlm1(1:kproma,:,krow)+xim1(1:kproma,:,krow)))
    zref(1:kproma)  = 0._dp
    zalpha          = 0._dp
    ! geopotential height above ground (not above NN)
    CALL geopot(geoph(:,2:klevp1),ztv,alnpr(:,:,krow),zalpha,zref,kbdim,kproma) ! half level geopotential
    geoph(:,1) = geoph(:,2) + 1.0_dp   ! for the top layer: set upper boundary to 1 m above its lower boundary
                                       ! so that it's thickness will be 1 m, note: there is no aerosol in the
                                       ! top layer anyhow

    ! calculate layer thickness and altitude above NN, both in units of m
    ! they are needed for the extinction coefficient, make sure height is non-negative
    ! k=1 is the top of the atmosphere
    DO jk=1,klev
       heightt(1:kproma,jk) = MAX(geoph(1:kproma,jk)  *rgrav,0.0_dp)
       heightb(1:kproma,jk) = MAX(geoph(1:kproma,jk+1)*rgrav,0.0_dp)
    ENDDO
    ! layer thickness in m 
    dheight(1:kproma,1:klev)= heightt-heightb

    ! reverse the vertical ordering because vertical grids of ECHAM6 and radiation code are reversed
    ! heightt_vr, heightb_vr and dheight_vr run from surface level to uppermost level
    DO jk=1,klev
       heightt_vr(1:kproma,jk) = heightt(1:kproma,klev-jk+1)
       heightb_vr(1:kproma,jk) = heightb(1:kproma,klev-jk+1)
       dheight_vr(1:kproma,jk) = dheight(1:kproma,klev-jk+1)
    ENDDO

    ! latitude and longitude 
    lat(1:kproma) = philat_2d(1:kproma,krow)  
    lon(1:kproma) = philon_2d(1:kproma,krow)  
    ! make sure longitude goes from 0 to 360 
    WHERE (lon(1:kproma).LT.0) lon = lon+360.0_dp 

    ! define wavelength bounds of the nb_sw solar bands; lambda1 and lambda2 are in units of nm
    lambda1 = 1.e7_dp/wavenum1(16:29);
    lambda2 = 1.e7_dp/wavenum2(16:29);
    lambda  = 0.5_dp*(lambda1+lambda2)    ! central wave length of spectral band in nm

    ! set the aerosol optical properties using the idealized Easy Aerosol climatology 
    CALL make_aop_easyaerosols(aod_vr,ssa_vr,asy_vr,heightt_vr,heightb_vr,dheight_vr,lat,lon,lambda,nb_sw,kproma,kbdim,klev,krow)

    ! write results to output arrays
    paer_tau_sw_vr(1:kproma,1:klev,1:nb_sw) = aod_vr(1:kproma,1:klev,1:nb_sw)
    paer_piz_sw_vr(1:kproma,1:klev,1:nb_sw) = ssa_vr(1:kproma,1:klev,1:nb_sw)
    paer_cg_sw_vr(1:kproma,1:klev,1:nb_sw)  = asy_vr(1:kproma,1:klev,1:nb_sw)
    paer_tau_lw_vr(1:kproma,1:klev,1:nb_lw) = 0.0_dp  ! aerosol not active in longwave domain 

    ! write to debug stream 
    ! additional diagnostic in debug stream: optical properties at 550nm
    IF (ldebugs) THEN
        ddf1(1:kproma,1:klev,krow) = paer_tau_sw_vr(1:kproma,1:klev,10)
        ddf2(1:kproma,1:klev,krow) = paer_piz_sw_vr(1:kproma,1:klev,10)
        ddf3(1:kproma,1:klev,krow) = paer_cg_sw_vr(1:kproma,1:klev,10)
        ddf4(1:kproma,1:klev,krow) = paer_tau_lw_vr(1:kproma,1:klev,10)
    ENDIF

  END SUBROUTINE set_aop_easyaerosols
  !---------------------------------------------------------------------------
  !>
  !! @brief Calculates the aerosol fields
  !!
  !! @remarks This routine is called by set_aop_easyaerosols. It uses the fits
  !! based on the anthropogenic aerosols of the Kinne year 2000 aerosol climatology
  !! to calculate the aerosol fields.
  !! Note on ordering of vertical layers:
  !! The vertical ordering of aod,ssa, and asy is the same as the vertical ordering of the 
  !! input arrays height and dheight. This means that the vertical ordering of the output
  !! is steered by the vertical ordering of the input.
  !
  SUBROUTINE make_aop_easyaerosols(aod,ssa,asy,heightt,heightb,dheight,lat,lon,lambda,nb_sw,kproma,kbdim,klev,krow)

    IMPLICIT NONE

    ! INPUT VARIABLES
    INTEGER,INTENT(in)  :: kproma,   & ! actual block length of arrays,
                           kbdim,    & ! maximum block length
                           klev,     & ! number of vertical levels
                           nb_sw,    & ! number of shortwave bands (visible + near infrared)
                           krow
    REAL(dp),INTENT(in),DIMENSION(kbdim,klev) :: heightt,   & ! geometrical height of top of layer in m
                                                 heightb,   & ! geometrical height of bottom of layer in m
                                                 dheight      ! geometrical thickness of layer in m         
 
    REAL(dp),INTENT(in),DIMENSION(kbdim)      :: lon         ! central longitude of the considered horizontal grid box
    REAL(dp),INTENT(in),DIMENSION(kbdim)      :: lat         ! central latitude of the considered horizontal grid box
    REAL(dp),INTENT(in),DIMENSION(nb_sw)      :: lambda      ! central wavelength of the shortwave bands in nm

    REAL(dp),INTENT(out),DIMENSION(kbdim,klev,nb_sw):: &
                        aod,   & !vertically resolved aerosol optical depth
                        ssa,   & !vertically resolved single scattering albedo 
                        asy      !vertically resolved asymmetry factor

    ! LOCAL VARIABLES
    REAL(dp),DIMENSION(kbdim,klev)   :: beta                  ! fraction of total aerosol in given layer
    REAL(dp),DIMENSION(kbdim,nb_sw)  :: aod_tot               ! total aerosol optical depth at surface
    REAL(dp)                         :: dlon1,dlon2,dlon3     ! zonal distance from center of plumes in deg lon 
    REAL(dp)                         :: aux                   ! auxiliary field for zonal variation of aod_tot
                                                              ! in case we use plumes 
    INTEGER                          :: jj,jk,jw              ! loop counters

    ! beta is the fraction of total aod assigned to layer and is unitless  
    ! that said, beta is the layer-integrated normalized extinction coefficient
    ! because the vertical profile is exponential, we can integrate analytically 
    ! if know the vertical boundaries of the level 
    beta(1:kproma,1:klev)=0.0_dp
    DO jk=1,klev
    DO jj=1,kproma
       ! constant close to surface 
       IF (heightt(jj,jk).LE.hexp) THEN
          beta(jj,jk) = a*exp(-c*hexp)*dheight(jj,jk)
       ! exponential decay above
       ELSEIF (heightb(jj,jk).GE.hexp) THEN
          beta(jj,jk) = a/c*(exp(-c*heightb(jj,jk))-exp(-c*heightt(jj,jk)))
       ! if height at which profile changes from contant to exponential is within layer 
       ELSE
          beta(jj,jk) =  a*exp(-c*hexp)*(hexp-heightb(jj,jk)) + a/c*(exp(-c*hexp)-exp(-c*heightt(jj,jk)))
       ENDIF
    ENDDO
    ENDDO

    ! make sure that vertical integral of beta is one 
    DO jj=1,kproma
          beta(jj,:) = beta(jj,:) / SUM(beta(jj,:))
    ENDDO

    ! total aerosol optical depth, here approximated as a Gaussian bellcurve, possibly with a longitudinal structure
    ! aod_tot is total aerosol optical depth for the given spectral band and the given lat-lon position
    DO jj=1,kproma
    DO jw=1,nb_sw
       ! Easy Aerosol concentrated in three zonal plumes
       IF (l3plumes) THEN 
          dlon1 = MIN(ABS(lon(jj)-lon1),ABS(360.0_dp-ABS(lon(jj)-lon1)))
          dlon2 = MIN(ABS(lon(jj)-lon2),ABS(360.0_dp-ABS(lon(jj)-lon2)))
          dlon3 = MIN(ABS(lon(jj)-lon3),ABS(360.0_dp-ABS(lon(jj)-lon3)))
          aux   = fact1*exp(-(dlon1**2)/(2*(siglon1**2))) + fact2*exp(-(dlon2**2)/(2*(siglon2**2))) + &  
                  fact3*exp(-(dlon3**2)/(2*(siglon3**2))) 
       ELSEIF (.NOT.l3plumes) THEN ! zonally-uniform Easy Aerosol
          aux = 1
       END IF
       aod_tot(jj,jw) = aod0 * (lambda(jw)/550.0_dp)**(-1*ang) * aux * exp((-(lat(jj)-lat0)**2)/(2*(siglat**2)))
    ENDDO
    ENDDO

    ! vertically resolved aerosol optical depth: aod_tot * beta
    ! vertical sum of aod equals aod_tot because beta is normalized
    DO jk=1,klev
    DO jw=1,nb_sw
       aod(1:kproma,jk,jw) = aod_tot(1:kproma,jw)*beta(1:kproma,jk)
    ENDDO
    ENDDO

    ! constant single scattering albedo and asymmetry factor
    ssa(1:kproma,1:klev,1:nb_sw) = ssa0
    asy(1:kproma,1:klev,1:nb_sw) = asy0

    ! write to debug stream 
    ! additional diagnostic in debug stream: beta, height, dheight, aod_tot at three wavelenghts
    IF (ldebugs) THEN
       ddf5(1:kproma,1:klev,krow)=beta(1:kproma,1:klev)
       ddf6(1:kproma,1:klev,krow)=dheight(1:kproma,1:klev)
       zdf1(1:kproma,krow)       =aod_tot(1:kproma,10)    !! at  534 nm
       zdf2(1:kproma,krow)       =aod_tot(1:kproma,13)    !! at  232 nm
       zdf3(1:kproma,krow)       =aod_tot(1:kproma,1)     !! at 3462 nm
    END IF
  
  END SUBROUTINE make_aop_easyaerosols

END MODULE mo_aero_easyaerosols
