!
! Contact: thomas.diehl@nasa.gov
!
! Fortran subroutine stub for reading the volcanic SO2 emission data into a 
! CTM or climate model.
! This routine assumes that the variables 
! -> nevents, jdn, elevation, cloud_column_height, so2 <-
! from the volcanic emission files have been read into the variables 
! -> nve, vjdn, velev, vcch, vso2 <-
! elsewhere. It also assumes that the grid indeces for each volcano have been 
! computed from the variables lon and lat in the emission files elsewhere and 
! stored in the variables vi and vj.
! It is also assumed that vso2 has been converted from kt SO2/event/day to
! kg SO2/event/day.
!
! imx: number of grid boxes in east-west direction
! jmx: number of grid boxes in north-south direction
! lmx: number of vertical grid boxes
! delz(imx,jmx,lmx): height of grid boxes [m]
! surface_altitude(imx,jmx): surface altitude [m]
! eso2_volc(imx,jmx,lmx): amount of so2 injected into a grid box [kg SO2/box/s]

! Module variables:
! LOGICAL :: levolc ! (volcanic emissions are computed if set to true)
! INTEGER :: idx1, next_idx1
! LOGICAL :: first = .true.
! INTEGER :: nve
! allocate these variables with dimension nve:
! INTEGER, ALLOCATABLE :: vi(:), vj(:), vjdn(:), vcch(:), velev(:)
! REAL,    ALLOCATABLE :: vso2(:), vlon(:), vlat(:)
! dt: dt(1)%jday contains the Julian Day Number at the beginning of the current
! model time step.
! jday_m1: Julian Day Number at beginning of previous time step

SUBROUTINE src_vso2_stub ( imx, jmx, lmx, &
                           delz, surface_altitude, &
                           eso2_volc )
 
  USE mo_control,      ONLY: levolc 
  USE mo_volc,         ONLY: nve, vjdn, vcch, velev, vi, vj, vso2, &
                             idx1, first, next_idx1
  USE mo_time_control, ONLY: dt, jday_m1
  USE mo_exception,    ONLY: finish

  IMPLICIT NONE

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

  INTEGER :: k, l, lv1, lv2, nvolc
  INTEGER :: idx1_tmp(1)
  REAL    :: hght, slab, slab1
  REAL    :: zh(0:lmx), dz(lmx)
  LOGICAL :: found
  CHARACTER(LEN=100) :: msg

  ! executable statements

! ****************************************************************************
! *  Calculate volcanic emission of SO2.                                     *
! ****************************************************************************

  eso2_volc(:,:,:) = 0.0

  if_volc_emis: IF (levolc) THEN
        
     IF (first) THEN
        nvolc = COUNT(vjdn == dt(1)%jday)
        IF (nvolc == 0) THEN
           msg = 'Current date not found in volcanic database.'
           CALL finish('source_su', msg)
        END IF
        idx1_tmp = MAXLOC(vjdn, MASK = vjdn == dt(1)%jday)
        idx1 = idx1_tmp(1)
     END IF

     IF (dt(1)%jday /= jday_m1) THEN
        idx1 = next_idx1
        jday_m1 = dt(1)%jday
     END IF
     
     first = .FALSE.

     found = .FALSE.

     nevents_loop: DO k = idx1, nevents

        ! all volcanic events of one day are stored sequentially; therefore, 
        ! once we encounter the first element which belongs to the next day 
        ! the loop can be terminated 
        IF (dt(1)%jday /= vjdn(k)) THEN
           next_idx1 = k
           EXIT
        END IF
        found = .TRUE.

        IF (velev(k) == vcch(k)) THEN

           ! the volcano inventory is constructed with vcch=velev for 
           ! non-eruptive degassing (no plume is assumed); 
           ! therefore, the emission is placed only in the level which 
           ! contains the crater elevation

           hght = REAL(velev(k))
           zh(0) = 0.0 + surface_altitude(vi(k),vj(k)) ! convert to m a.s.l.

           IF (hght < zh(0)) THEN 
              eso2_volc(vi(k),vj(k),1) = &
                   eso2_volc(vi(k),vj(k),1) + vso2(k)/(24.0*3600.0)
           ELSE
              DO l = 1,lmx
                 dz(l) = delz(vi(k),vj(k),lmx-l+1) ! delz: TOA->sfc
                 zh(l) = zh(l-1) + dz(l)
                 IF (zh(l-1) <= hght .AND. zh(l) > hght) THEN 
                    eso2_volc(vi(k),vj(k),l) = eso2_volc(vi(k),vj(k),l) + &
                         vso2(k)/(24.0*3600.0)
                    EXIT
                 END IF
              END DO
           END IF

        ELSE

           ! --------------------------------------------------
           !   Define a slab at the top 1/3 of the volcano plume.
           ! --------------------------------------------------
           hght  = REAL(vcch(k))
           slab1 = hght - (hght - REAL(velev(k)))/3.0   ! slab bottom height
           slab  = hght - slab1                         ! slab thickness
           zh(0) = 0.0 + surface_altitude(vi(k),vj(k))  ! convert to m a.s.l.
       
           IF (hght <= zh(0)) THEN 
              lv1 = 1
              lv2 = 1
              dz(1) = slab
           ELSE
              level_loop: DO l = 1,lmx
                 dz(l) = delz(vi(k),vj(k),lmx-l+1)      ! delz: TOA->sfc
                 zh(l) = zh(l-1) + dz(l)                ! top edge height
                 IF (l == lmx .AND. hght > zh(l)) THEN  ! max model erup.height
                    hght = zh(l)
                    slab1 = slab1 - (hght - zh(l))
                 END IF
                 
                 ! --------------------------------------------------
                 !   If Zh(l) <= bottom of the slab, go to next level.
                 ! --------------------------------------------------
                 IF (zh(l) <= slab1) CYCLE
                 
                 ! --------------------------------------------------
                 !   If the slab is only in current level:
                 ! --------------------------------------------------
                 IF (zh(l-1) <= slab1 .AND. zh(l) >= hght) THEN
                    lv1 = l
                    lv2 = l
                    dz(l) = slab
                    EXIT
                    
                    ! --------------------------------------------------
                    !   The slab extends more than one level.  Find the 
                    !   lowest (lv1) and the highest (lv2) levels: 
                    ! --------------------------------------------------
                 ELSE IF (zh(l-1) < slab1 .AND. zh(l) > slab1)  THEN
                    lv1 = l
                    dz(l) = zh(l) - slab1
                    
                 ELSE IF (zh(l-1) <= hght  .AND. zh(l) > hght)  THEN
                    lv2 = l
                    dz(l) = hght - zh(l-1)
                    EXIT     
                 END IF
                 
              END DO level_loop

           END IF

           ! --------------------------------------------------
           !   Calculate SO2 emission in the levels between  
           !   lv1 and lv2.  Convert vso2 from kgSO2/box/day
           !   to kgSO2/box/sec.  eso2_volc is distributed evenly
           !   with altitude among the slab.
           !   Summation is necessary since more than one volcano
           !   could be located within the same grid box.
           ! --------------------------------------------------
           DO l = lv1,lv2
              eso2_volc(vi(k),vj(k),l) = eso2_volc(vi(k),vj(k),l) &
                   + vso2(k)/(24.0*3600.0) * dz(l)/slab
           END DO

        END IF
           
     END DO nevents_loop
           
     IF (.NOT. found) THEN
        msg = 'Current date not found in volcanic database.'
        CALL finish('source_su', msg)
     END IF

  END IF if_volc_emis
        
END SUBROUTINE src_vso2_stub
