!!---------------------------------------------------------------
!! Author: Yaxun Tang, tang@sep.stanford.edy
!!---------------------------------------------------------------
module base_outgrd_mod
  use base_struct_mod
  implicit none
contains
  !!.. igrid has to be initialized before calling this subroutine
  subroutine base_outgrd_gen3d(ogrid, igrid, pgrid)
    type(out_grid) :: ogrid
    type(img_grid) :: igrid
    type(pro_grid) :: pgrid

    real           :: aix_max, aiy_max, aiz_max
    real           :: aix_min, aiy_min, aiz_min
    real           :: apx_max, apy_max, apz_max
    real           :: apx_min, apy_min, apz_min
    real           :: aix_length, aiy_length, aiz_length
    real           :: apx_length, apy_length, apz_length

    ogrid%ax  = igrid%ax
    ogrid%ay  = igrid%ay
    ogrid%az  = igrid%az
    ogrid%ahx = igrid%ahx
    ogrid%ahy = igrid%ahy
    ogrid%ahz = igrid%ahz
    ogrid%aht = igrid%aht

    aix_min = igrid%ax%o
    apx_min = pgrid%ax%o
    aix_max = igrid%ax%o + (igrid%ax%n-1)*igrid%ax%d
    apx_max = pgrid%ax%o + (pgrid%ax%n-1)*pgrid%ax%d

    aiy_min = igrid%ay%o
    apy_min = pgrid%ay%o
    aiy_max = igrid%ay%o + (igrid%ay%n-1)*igrid%ay%d
    apy_max = pgrid%ay%o + (pgrid%ay%n-1)*pgrid%ay%d

    aiz_min = igrid%az%o
    apz_min = pgrid%az%o
    aiz_max = igrid%az%o + (igrid%az%n-1)*igrid%az%d
    apz_max = pgrid%az%o + (pgrid%az%n-1)*pgrid%az%d


    if (pgrid%ax%n > 1) then
       if(apx_min < aix_min) then
          if(apx_max < aix_min) then
             ogrid%ax%n = 0
          else if(apx_max >= aix_min .and. apx_max < aix_max) then
             ogrid%ax%o = igrid%ax%o
             ogrid%ax%n = 1.5 + (apx_max-ogrid%ax%o)/ogrid%ax%d
          else
             ogrid%ax%o = igrid%ax%o
             ogrid%ax%n = 1.5 + (aix_max-ogrid%ax%o)/ogrid%ax%d
          end if
       else if(apx_min >= aix_min .and. apx_min <= aix_max) then
          if(apx_max < aix_max) then
             ogrid%ax%o = pgrid%ax%o
             ogrid%ax%n = 1.5 + (apx_max-ogrid%ax%o)/ogrid%ax%d
          else
             ogrid%ax%o = pgrid%ax%o
             ogrid%ax%n = 1.5 + (aix_max-ogrid%ax%o)/ogrid%ax%d
          end if
       else
          ogrid%ax%n = 0
       end if
    end if

    if (pgrid%ay%n > 1) then
       if(apy_min < aiy_min) then
          if(apy_max < aiy_min) then
             ogrid%ay%n = 0
          else if(apy_max >= aiy_min .and. apy_max < aiy_max) then
             ogrid%ay%o = igrid%ay%o
             ogrid%ay%n = 1.5 + (apy_max-ogrid%ay%o)/ogrid%ay%d
          else
             ogrid%ay%o = igrid%ay%o
             ogrid%ay%n = 1.5 + (aiy_max-ogrid%ay%o)/ogrid%ay%d
          end if
       else if(apy_min >= aiy_min .and. apy_min <= aiy_max) then
          if(apy_max < aiy_max) then
             ogrid%ay%o = pgrid%ay%o
             ogrid%ay%n = 1.5 + (apy_max-ogrid%ay%o)/ogrid%ay%d
          else
             ogrid%ay%o = pgrid%ay%o
             ogrid%ay%n = 1.5 + (aiy_max-ogrid%ay%o)/ogrid%ay%d
          end if
       else
          ogrid%ay%n = 0
       end if
    end if


    if (pgrid%az%n > 1) then
       if(apz_min < aiz_min) then
          if(apz_max < aiz_min) then
             ogrid%az%n = 0
          else if(apz_max >= aiz_min .and. apz_max < aiz_max) then
             ogrid%az%o = igrid%az%o
             ogrid%az%n = 1.5 + (apz_max-ogrid%az%o)/ogrid%az%d
          else
             ogrid%az%o = igrid%az%o
             ogrid%az%n = 1.5 + (aiz_max-ogrid%az%o)/ogrid%az%d
          end if
       else if(apz_min >= aiz_min .and. apz_min <= aiz_max) then
          if(apz_max < aiz_max) then
             ogrid%az%o = pgrid%az%o
             ogrid%az%n = 1.5 + (apz_max-ogrid%az%o)/ogrid%az%d
          else
             ogrid%az%o = pgrid%az%o
             ogrid%az%n = 1.5 + (aiz_max-ogrid%az%o)/ogrid%az%d
          end if
       else
          ogrid%az%n = 0
       end if
    end if

    !!.. computing the first and the last x, y index of the output image on the final image grid
    ogrid%index_outimg_xmin = 1.5 + (ogrid%ax%o - igrid%ax%o)/igrid%ax%d
    ogrid%index_outimg_ymin = 1.5 + (ogrid%ay%o - igrid%ay%o)/igrid%ay%d
    ogrid%index_outimg_zmin = 1.5 + (ogrid%az%o - igrid%az%o)/igrid%az%d
    ogrid%index_outimg_xmax = ogrid%index_outimg_xmin + ogrid%ax%n - 1
    ogrid%index_outimg_ymax = ogrid%index_outimg_ymin + ogrid%ay%n - 1
    ogrid%index_outimg_zmax = ogrid%index_outimg_zmin + ogrid%az%n - 1

    !!.. computing the first and the last x, y index of the output image on the propagation grid
    ogrid%index_outpro_xmin = 1.5 + (ogrid%ax%o - pgrid%ax%o)/pgrid%ax%d
    ogrid%index_outpro_ymin = 1.5 + (ogrid%ay%o - pgrid%ay%o)/pgrid%ay%d
    ogrid%index_outpro_zmin = 1.5 + (ogrid%az%o - pgrid%az%o)/pgrid%az%d
    ogrid%index_outpro_xmax = ogrid%index_outpro_xmin + ogrid%ax%n - 1
    ogrid%index_outpro_ymax = ogrid%index_outpro_ymin + ogrid%ay%n - 1
    ogrid%index_outpro_zmax = ogrid%index_outpro_zmin + ogrid%az%n - 1

    !!.. compute the maximum extrapolation depth steps
    pgrid%nzstep = ogrid%index_outpro_zmax
    pgrid%nzstep = min(pgrid%nzstep, pgrid%az%n)
    return
  end subroutine
end module
