! ------------------------------------------------------------------------
!
! Contact: thomas.diehl@nasa.gov
!
! ------------------------------------------------------------------------
!
! Fortran subroutine stub for converting fuel burned to aircraft emissions
! of SO2, SO4, BC and OC.
! This routine makes the following assumptions:
! 1. The variables 
!    -> ilev, lev, fuel <-
!    from the aircraft emission files have been read into the variables 
!    -> alt_lev_int, alt_lev_mid, fuel_alt <-
! 2. The dimension 
!    -> lev <-
!    from the aircraft emission files has been read into the variable 
!    -> nlev <-
! 3. ilev and lev have been converted from km to m.
!
! imx: number of model grid boxes in east-west direction
! jmx: number of model grid boxes in north-south direction
! lmx: number of model vertical grid boxes
! nmx: number of model species, with
!      n=NBC1 hydrophobic BC, n=NBC2 hydrophilic BC, 
!      n=NOC1 hydrophobic OC, n=NOC2 hydrophilic OC
!      n=NSO2 SO2,            n=NSO4 SO4
! delz(imx,jmx,lmx): height of model grid boxes [m]
! surface_altitude(imx,jmx): surface altitude above sea level [m]
! emi_ac(imx,jmx,lmx,nmx): amount of species injected into a grid 
!                          box in [mol/box/s] or [kg/box/s]

! Module variables:
! LOGICAL :: laircraft ! aircraft emissions are computed if set to true
! REAL, PARAMETER :: cmw = 12.0 ! molecular weight of carbon
!
! dt: dt(1)%yr contains the year at the beginning of the current
!     model time step.

SUBROUTINE aircraft_emissions_stub ( imx, jmx, lmx, nmx, &
                                     delz, surface_altitude, &
                                     emi_ac )

  USE mo_control,      ONLY: laircraft
  USE mo_aircraft,     ONLY: nlev, alt_lev_int, alt_lev_mid, fuel_alt
  USE mo_tracer,       ONLY: NBC1, NBC2, NOC1, NOC2, NSO2, NSO4
  USE mo_constants,    ONLY: cmw
  USE mo_time_control, ONLY: dt
  USE mo_exception,    ONLY: finish

  IMPLICIT NONE

  INTEGER, INTENT(IN)    :: imx, jmx, lmx, nmx 
  REAL,    INTENT(IN)    :: delz(imx,jmx,lmx) 
  REAL,    INTENT(IN)    :: surface_altitude(imx,jmx) 
  REAL,    INTENT(OUT)   :: emi_ac(imx,jmx,lmx,nmx)

  INTEGER :: i, j, k1, k2
  REAL, PARAMETER :: fac = 0.75 ! BC fraction of total carbon from aircraft 
                                ! emissions
  REAL :: ei(18), ei_hght(18), ei_int
  DATA ei      /   0.08,  0.096,   0.087,   0.075,   0.068,   0.062, &
                   0.058, 0.043,   0.042,   0.041,   0.031,   0.025, &
                   0.019, 0.048,   0.068,   0.069,   0.078,   0.082 /
  DATA ei_hght /   0.0,   500.0,  1500.0,  2500.0,  3500.0,  4500.0, &
                5500.0,  6500.0,  7500.0,  8500.0,  9500.0, 10500.0, &
               11500.0, 12500.0, 13500.0, 14500.0, 15500.0, 16500.0 /
  REAL,    SAVE :: ei_mid(17)
  LOGICAL, SAVE :: lfirst=.TRUE.
  LOGICAL :: found
  REAL    :: z_bottom(lmx), z_top(lmx), ovlap
  REAL    :: delz_alt_lev(nlev)
  REAL    :: fuel(imx,jmx,lmx)
  REAL    :: ei_so2
  CHARACTER(LEN=100) :: msg

  ! executable statements

! ****************************************************************************
! *  Calculate aircraft emissions of SO2, SO4, BC, and OC                    *
! ****************************************************************************

  if_aircraft_emis: IF (laircraft) THEN

     IF (lfirst) THEN
        DO k2 = 1,SIZE(ei_mid)
           ei_mid(k2) = (ei(k2) + ei(k2+1))*0.5
        END DO
       lfirst = .FALSE.
     END IF

     DO k2 = 1,SIZE(alt_lev_mid)
        delz_alt_lev(k2) = alt_lev_int(k2+1) - alt_lev_int(k2)
     END DO

     ! no hydrophobic BC/OC
     emi_ac(:,:,:,NBC1) = 0.0
     emi_ac(:,:,:,NOC1) = 0.0

     ! regrid fuel_alt (altitude levels) to fuel (model levels)
     ! unit: kg/box/d
     ! Aircraft emissions are an extensive quantity (they are given in kg,
     ! not as a concentration), so we have to add the values in the fine
     ! resolution boxes to get the value in the coarse grid box (instead of
     ! averaging) 
     DO j = 1,jmx
        DO i = 1,imx
           
           z_bottom(1) = 0.0 + surface_altitude(i,j) ! convert to m a.s.l.
           lfound = .FALSE.
           DO k1 = 1,lmx

              z_top(k1) = z_bottom(k1) + delz(i,j,k1)
              IF (k1 < lmx) z_bottom(k1+1) = z_top(k1)

              IF (z_bottom(k1) >= alt_lev_int(SIZE(alt_lev_int))) THEN
                 fuel(i,j,k1) = 0.0
                 lfound = .TRUE.
              END IF

              IF (.NOT. lfound) THEN
                 ! compute overlap of current box with each of the boxes from 
                 ! the altitude grid
                 fuel(i,j,k1) = 0.0
                 DO k2 = 1,SIZE(alt_lev_mid)
                    ! condition for ovlap > 0.0
                    IF (alt_lev_int(k2)   < z_top(k1) .AND. &
                        alt_lev_int(k2+1) > z_bottom(k1)) THEN
                       IF (z_top(k1) >= alt_lev_int(k2+1)) THEN
                          ovlap = MIN(alt_lev_int(k2+1) - z_bottom(k1), &
                                      delz(i,j,k1), &
                                      alt_lev_int(k2+1) - alt_lev_int(k2) )
                       ELSE ! ztop(k1) < alt_lev_int(k2+1)
                          ovlap = MIN(z_top(k1) - alt_lev_int(k2), &
                                      delz(i,j,k1), &
                                      alt_lev_int(k2+1) - alt_lev_int(k2) )
                       END IF
                       fuel(i,j,k1) = fuel(i,j,k1) + &
                                      ovlap/delz_alt_lev(k2)*fuel_alt(i,j,k2)
                    END IF
                 END DO
                 
              END IF
           END DO
        END DO
     END DO

     ! BC aircraft emission is calculated by assuming a height-dependent 
     ! emission index EI in g/(kg fuel burned)  
     ! Further assumptions:
     ! - there is no diurnal variation of the emission rate
     ! - approx 75% of total carbon is emitted as black carbon
     ! - all BC and OC is emitted as hydrophilic
     ! Finally, convert the units from g/box/d to desired units.
     ! (here we have [mol/box/s] for BC and OC and [kg/box/s] for SO2 and SO4)

     DO j = 1,jmx
        DO i = 1,imx
           
           z_bottom(1) = 0.0 + surface_altitude(i,j) ! convert to m a.s.l.

           lfound = .FALSE.
           DO k1 = 1,lmx

              z_top(k1) = z_bottom(k1) + delz(i,j,k1)
              IF (k1 < lmx) z_bottom(k1+1) = z_top(k1)

              IF (z_bottom(k1) >= ei_hght(SIZE(ei_hght))) THEN
                 ei_int = ei(SIZE(ei))
                 lfound = .TRUE.
              END IF

              IF (.NOT. lfound) THEN
                 ! compute overlap of current box with each of the boxes from 
                 ! the Doepelheuer grid
                 ei_int = 0.0
                 DO k2 = 1,SIZE(ei_mid)
                    ! condition for ovlap > 0.0
                    IF (ei_hght(k2)   < z_top(k1) .AND. &
                        ei_hght(k2+1) > z_bottom(k1)) THEN
                       IF (z_top(k1) >= ei_hght(k2+1)) THEN
                          ovlap = MIN(ei_hght(k2+1) - z_bottom(k1), &
                                      delz(i,j,k1), &
                                      ei_hght(k2+1) - ei_hght(k2) )
                       ELSE ! ztop(k1) < ei_hght(k2+1)
                          ovlap = MIN(z_top(k1) - ei_hght(k2), &
                                      delz(i,j,k1), &
                                      ei_hght(k2+1) - ei_hght(k2) )
                       END IF
!                       ovlap = MIN ( z_top(k1)      - ei_hght(k2), &
!                                     ei_hght(k2+1) - z_bottom(k1), &
!                                     delz(i,j,k1), &
!                                     ei_hght(k2+1) - ei_hght(k2) )
                       ei_int = ei_int + ovlap/delz(i,j,k1)*ei_mid(k2)
                    END IF
                 END DO

                 IF (ei_int == 0.0) THEN
                    msg = 'Error occurred while computing ei_int.'
                    CALL finish('aircraft_emissions_stub', msg)
                 END IF
                 
                 IF (z_top(k1) > ei_hght(SIZE(ei_hght))) THEN
                    ovlap = ABS(z_top(k1) - ei_hght(SIZE(ei_hght))) 
                    ei_int = ei_int + ovlap/delz(i,j,k1)*ei(SIZE(ei))
                 END IF
                 
              END IF

              ! in this example: final units are [mol(species)/box/s]
              emi_ac(i,j,k1,NBC2) = fuel(i,j,k1)*ei_int /cmw/(24.0*3600.0) 
              emi_ac(i,j,k1,NOC2) = (1.0-fac)/fac * &
                                    fuel(i,j,k1)*ei_int /cmw/(24.0*3600.0)
           END DO

        END DO
     END DO

     ! SO2 emission is calculated by assuming an emission index EI 
     ! [g(species)/kg(fuel burned] of 0.8 up to 1999 and of 0.4 for 2015. The
     ! EI is linearly interpolated between these years.
     ! It is also assumed that there is no diurnal variation of the emission 
     ! rate.

     IF (dt(1)%yr <= 1999) THEN
        ei_so2 = 0.8
     ELSE IF (dt(1)%yr == 2015) THEN
        ei_so2 = 0.4
     ELSE
        ei_so2 = 0.8 + (0.4 - 0.8)*(dt(1)%yr-1999)/(2015-1999)
     END IF

     ! in this example: final units are [kg(species)/box/s]
     emi_ac(:,:,:,NSO2) = fuel(:,:,:) * ei_so2 * 0.995 * 1.0E-3 / (24.0*3600.0)
     emi_ac(:,:,:,NSO4) = fuel(:,:,:) * ei_so2 * 0.015 * 1.0E-3 / (24.0*3600.0)

  ELSE

     emi_ac(:,:,:,NBC1) = 0.0
     emi_ac(:,:,:,NOC1) = 0.0

     emi_ac(:,:,:,NBC2) = 0.0
     emi_ac(:,:,:,NOC2) = 0.0

     emi_ac(:,:,:,NSO2) = 0.0
     emi_ac(:,:,:,NSO4) = 0.0

  END IF if_aircraft_emis

END SUBROUTINE aircraft_emissions_stub
