module wei_greens_mod
  use sep
  use wei_type_mod
  use wei_ssf_mod
  use wei_vel_mod
  use wei_taper_mod

  implicit none
  complex, dimension(:,:,:)  , pointer, private :: wfld_d
  complex, dimension(:)      , pointer, private :: csou
  real   , dimension(:,:,:)  , pointer, private :: vmod
  real   , dimension(:,:)    , pointer, private :: vref
  real   , dimension(:)      , pointer, private :: kx2, ky2
  complex, dimension(:,:)    , pointer, private :: wfld_slice_d
  real   , dimension(:,:)    , pointer, private :: vmod_slice
  real   , dimension(:)      , pointer, private :: vref_slice
  integer                             , private :: nws_bak
  integer                             , private :: nzs_bak

  type(pro_grid)                      , private :: pgrid
  type(vel_grid)                      , private :: vgrid
  type(gre_grid)                      , private :: ggrid
  type(block)                         , private :: wblock
  type(block)                         , private :: zblock

  logical                             , private :: lloyd
  logical                             , private :: decon
  logical                             , private :: radian
  logical                             , private :: upward
  logical                             , private :: spike
  logical                             , private :: taper
  real                                , private :: taper_pct
  integer                             , private :: taper_type
  real                                , private :: maxang
  integer                             , private :: nref
  integer                             , private :: opsign

  integer                       , private, save :: num_report
  integer                       , private, save :: istep
  integer                       , private, save :: rstep
  integer                       , private, save :: astep
  integer                       , private, save :: nstep


  real   , parameter                  , private :: PI = 3.14159265359
  logical, parameter                  , private :: T = .true., F = .false.
  integer, external                   , private :: sreed_window, srite_window, getch
  character(len=256)                  , private :: sou_="sou", vel_="vel", greens_="greens"
contains


  integer function wei_greens_compute_init() result(stat)
    call wei_greens_param_init()
    call wei_greens_input_history()
    call wei_greens_output_history()
    call wei_greens_radians_calc()
    !call wei_greens_io_init()
    call wei_greens_block_calc(wblock, ggrid%aw)
    call wei_greens_propagation_grid()
    call wei_greens_velocity_grid()
    call wei_greens_k2_map()
    call wei_greens_wfld_d_alloc(wblock%ns)
    call wei_greens_velmod_alloc(zblock%ns)

    stat = wei_vel_init(lloyd, nref)
    stat = wei_ssf_init(kx2, ky2, pgrid%az%d, maxang)

    if (taper)  stat = wei_taper_2d_init(taper_type, pgrid%nxpad_beg, pgrid%nxpad_end, pgrid%nypad_beg, pgrid%nypad_end)

    nws_bak = wblock%ns
    nzs_bak = zblock%ns

    stat = 0

  end function
!========================================================================
  integer function wei_greens_compute_oper() result(stat)
    integer       :: imodx, imody, imodz
    real          :: modx, mody, modz
    type(axis)    :: amodz

    do imodz=1, ggrid%az_modl%n
write(0,*)"imodz = ", imodz
       modz = ggrid%az_modl%o + (imodz-1)*ggrid%az_modl%d
       ggrid%index_modpro_modz = 1.5 + (modz - vgrid%az%o) / vgrid%az%d
       !!.. propagate wavefields upward starting from the selected model points
       if (upward) then
          amodz%o = vgrid%az%o; amodz%d = vgrid%az%d; amodz%n = ggrid%index_modpro_modz
       !!.. propagate wavefields downward starting from the selected model points
       else
          amodz%o = modz; amodz%d = vgrid%az%d; amodz%n = vgrid%az%n - ggrid%index_modpro_modz + 1
       end if
       call wei_greens_block_calc(zblock, amodz)
       do imody=1, ggrid%ay_modl%n
          mody = ggrid%ay_modl%o + (imody-1)*ggrid%ay_modl%d
          ggrid%index_modpro_mody = 1.5 + (mody - pgrid%ay%o) / pgrid%ay%d

          do imodx=1, ggrid%ax_modl%n
             modx = ggrid%ax_modl%o + (imodx-1)*ggrid%ax_modl%d
             ggrid%index_modpro_modx = 1.5 + (modx - pgrid%ax%o) / pgrid%ax%d

             stat = wei_greens_wfld(imodx, imody, imodz)     
          end do
       end do
    end do

  end function
!=======================================================================
  integer function wei_greens_compute_clos() result(stat)
    stat = wei_vel_clos()
    stat = wei_ssf_clos()
    call wei_greens_wfld_d_dealloc()
    call wei_greens_velmod_dealloc()
  end function
!-------------------------------------------------------------------------
! wavefield propagation for one point in the model space
!-------------------------------------------------------------------------
  integer function wei_greens_wfld(imodx, imody, imodz) result(stat)
    real       :: modx, mody, modz
    integer    :: imodx, imody, imodz
    integer    :: iwb, iws, izs, izb
    type(axis) :: aws, azs
    real       :: ws, zs
    integer    :: nzb_beg, nzb_end, nzb_step
    integer    :: nzs_beg, nzs_end, nzs_step
    !!.. upword propagation or downward propatation
    if (upward) then
       nzb_beg = zblock%nb; nzb_end = 1; nzb_step = -1
    else
       nzb_beg = 1; nzb_end = zblock%nb; nzb_step = +1
    end if  
    

    do iwb=1, wblock%nb
       aws%o = wblock%os(iwb); aws%d = wblock%ds; aws%n = wblock%blen(iwb)

       if (aws%n /= nws_bak) then
           call wei_greens_wfld_d_dealloc()
           call wei_greens_wfld_d_alloc(aws%n)
           nws_bak = aws%n
       end if

       if (spike) then 
          csou = 1.0
       else
          call wei_greens_io_csou(wblock%indx(iwb), aws%n)
       end if
       call wei_greens_wfld_d_init(aws%n)
       !if (taper) stat = wei_taper_2d_oper(wfld_d, taper_pct)

       do izb=nzb_beg, nzb_end, nzb_step
          azs%o = zblock%os(izb); azs%d = zblock%ds; azs%n = zblock%blen(izb)
          if (azs%n /= nzs_bak) then
             call wei_greens_velmod_dealloc()
             call wei_greens_velmod_alloc(azs%n)
             nzs_bak = azs%n
          end if

          call wei_greens_io_velo(zblock%indx(izb), azs%n)

          if (upward) then
             nzs_beg = azs%n; nzs_end = 1; nzs_step = -1
          else
             nzs_beg = 1; nzs_end = azs%n; nzs_step = +1
          end if

          do izs=nzs_beg, nzs_end, nzs_step
             zs = azs%o + (izs-1)*azs%d
             vmod_slice => vmod(:,:,izs)
             vref_slice => vref(:,izs)
             stat = wei_vel_oper(vmod_slice, vref_slice)
             if (taper) stat = wei_taper_2d_oper(wfld_d, taper_pct)

             do iws=1, aws%n
                ws = aws%o + (iws-1)*aws%d
                wfld_slice_d => wfld_d(:,:,iws)
                stat = wei_ssf_oper(T, +1, ws, vref_slice(1), vmod_slice, wfld_slice_d) 
             end do

          end do
       end do 
       call wei_greens_io_wfld("rite", wfld_d, aws%n, wblock%indx(iwb), imodx, imody, imodz)
    end do

  end function

!=========================================================================
  subroutine wei_greens_param_init()
    call sep_init(SOURCE)
    call getch_add_string("head=/dev/null")
    call getch_add_string("noheader=y")
    call set_no_putch()
    call sep_begin_prog()

    call from_param("nx_modl", ggrid%ax_modl%n,  1)
    call from_param("dx_modl", ggrid%ax_modl%d, 0.)
    call from_param("ox_modl", ggrid%ax_modl%o, 0.)

    call from_param("ny_modl", ggrid%ay_modl%n,  1)
    call from_param("dy_modl", ggrid%ay_modl%d, 0.)
    call from_param("oy_modl", ggrid%ay_modl%o, 0.)

    call from_param("nz_modl", ggrid%az_modl%n,  1)
    call from_param("dz_modl", ggrid%az_modl%d, 0.)
    call from_param("oz_modl", ggrid%az_modl%o, 500.)

    call from_param("nref"      ,       nref,       1)
    call from_param("lloyd"     ,      lloyd, .false.)
    call from_param("maxang"    ,     maxang,     90.)
    call from_param("radian"    ,     radian, .false.)
    call from_param("decon"     ,      decon, .false.)
    call from_param("upward"    ,     upward,  .true.)
    call from_param("taper"     ,      taper,  .true.)
    call from_param("taper_type", taper_type,       1)
    call from_param("taper_pct" ,  taper_pct,     0.1)
  end subroutine

!===============================================================
  subroutine wei_greens_input_history()
    real    :: surface_xmin_default
    real    :: surface_xmax_default
    real    :: surface_ymin_default
    real    :: surface_ymax_default
    integer :: num_match
    character(len=256) :: string

    spike = .true.

    num_match = getch(sou_, "s", string)

    if (num_match == 1) spike = .false.

    if (spike) then
       call from_param("nw", ggrid%aw%n,  1)
       call from_param("dw", ggrid%aw%d, 1.)
       call from_param("ow", ggrid%aw%o, 0.)
    else
       call auxin(sou_)
       call from_aux(sou_, "n3", ggrid%aw%n)
       call from_aux(sou_, "d3", ggrid%aw%d)
       call from_aux(sou_, "o3", ggrid%aw%o)
    end if

    call auxin(vel_)
    call from_aux(vel_, "n1", vgrid%ax%n)
    call from_aux(vel_, "d1", vgrid%ax%d)
    call from_aux(vel_, "o1", vgrid%ax%o)

    call from_aux(vel_, "n2", vgrid%ay%n)
    call from_aux(vel_, "d2", vgrid%ay%d)
    call from_aux(vel_, "o2", vgrid%ay%o)

    call from_aux(vel_, "n3", vgrid%az%n)
    call from_aux(vel_, "d3", vgrid%az%d)
    call from_aux(vel_, "o3", vgrid%az%o)

    surface_xmin_default = vgrid%ax%o
    surface_xmax_default = vgrid%ax%o + (vgrid%ax%n-1)*vgrid%ax%d

    surface_ymin_default = vgrid%ay%o
    surface_ymax_default = vgrid%ay%o + (vgrid%ay%n-1)*vgrid%ay%d

!! currently assumes the receiver grids to be the same as the velocity grid
    call from_param("surface_xmin", ggrid%surface_xmin, surface_xmin_default)
    call from_param("surface_xmax", ggrid%surface_xmax, surface_xmax_default)
    call from_param("surface_ymin", ggrid%surface_ymin, surface_ymin_default)
    call from_param("surface_ymax", ggrid%surface_ymax, surface_ymax_default)


    call from_param("nxpad_beg", pgrid%nxpad_beg, 0)
    call from_param("nxpad_end", pgrid%nxpad_end, 0)
    call from_param("nypad_beg", pgrid%nypad_beg, 0)
    call from_param("nypad_end", pgrid%nypad_end, 0)

    call from_param("nws", wblock%ns, ggrid%aw%n)
    call from_param("nzs", zblock%ns, vgrid%az%n)
 end subroutine

!===============================================================
 subroutine wei_greens_output_history()

    ggrid%ax_surf%o = ggrid%surface_xmin
    ggrid%ax_surf%d = vgrid%ax%d
    ggrid%ax_surf%n = 1.5 + (ggrid%surface_xmax - ggrid%surface_xmin) / ggrid%ax_surf%d
    
    ggrid%ay_surf%o = ggrid%surface_ymin
    ggrid%ay_surf%d = vgrid%ay%d
    ggrid%ay_surf%n = 1.5 + (ggrid%surface_ymax - ggrid%surface_ymin) / ggrid%ay_surf%d


    call auxinout(greens_)
    call auxputch("n1", "i", ggrid%ax_surf%n, greens_)
    call auxputch("d1", "r", ggrid%ax_surf%d, greens_)
    call auxputch("o1", "r", ggrid%ax_surf%o, greens_)

    call auxputch("n2", "i", ggrid%ay_surf%n, greens_)
    call auxputch("d2", "r", ggrid%ay_surf%d, greens_)
    call auxputch("o2", "r", ggrid%ay_surf%o, greens_)

    call auxputch("n3", "i",      ggrid%aw%n, greens_)
    call auxputch("d3", "r",      ggrid%aw%d, greens_)
    call auxputch("o3", "r",      ggrid%aw%o, greens_)

    call auxputch("n4", "i", ggrid%ax_modl%n, greens_)
    call auxputch("d4", "r", ggrid%ax_modl%d, greens_)
    call auxputch("o4", "r", ggrid%ax_modl%o, greens_)

    call auxputch("n5", "i", ggrid%ay_modl%n, greens_)
    call auxputch("d5", "r", ggrid%ay_modl%d, greens_)
    call auxputch("o5", "r", ggrid%ay_modl%o, greens_)

    call auxputch("n6", "i", ggrid%az_modl%n, greens_)
    call auxputch("d6", "r", ggrid%az_modl%d, greens_)
    call auxputch("o6", "r", ggrid%az_modl%o, greens_)

    call auxputch("esize" , "i",         8   , greens_)
    call auxputch("label1", "s", "surface_x" , greens_)
    call auxputch("label2", "s", "surface_y" , greens_)
    call auxputch("label3", "s", "frequency" , greens_)
    call auxputch("label4", "s", "model_x"   , greens_)
    call auxputch("label5", "s", "model_y"   , greens_)
    call auxputch("label6", "s", "model_z"   , greens_)

  end subroutine
!-----------------------------------------------------------------------
!------------------------------------------------------------------------
  subroutine wei_greens_propagation_grid()
    real    :: pgrid_xmin, pgrid_xmax
    real    :: pgrid_ymin, pgrid_ymax

    pgrid%ax%d  = vgrid%ax%d
    pgrid%ay%d  = vgrid%ay%d
    pgrid%az%d  = vgrid%az%d
 
    pgrid%gxmin = ggrid%surface_xmin
    pgrid%gxmax = ggrid%surface_xmax
    pgrid%gymin = ggrid%surface_ymin
    pgrid%gymax = ggrid%surface_ymax


    pgrid_xmin = pgrid%gxmin; pgrid_xmax = pgrid%gxmax
    pgrid_ymin = pgrid%gymin; pgrid_ymax = pgrid%gymax
    !!.. add the padding on both x and y direction
    pgrid_xmin = pgrid_xmin - pgrid%nxpad_beg*pgrid%ax%d
    pgrid_xmax = pgrid_xmax + pgrid%nxpad_end*pgrid%ax%d

    pgrid_ymin = pgrid_ymin - pgrid%nypad_beg*pgrid%ay%d
    pgrid_ymax = pgrid_ymax + pgrid%nypad_end*pgrid%ay%d

    !!.. make sure the propagation grid is within the velocity model 
    pgrid_xmin = max(pgrid_xmin, vgrid%ax%o)
    pgrid_xmax = min(pgrid_xmax, vgrid%ax%o + (vgrid%ax%n-1)*vgrid%ax%d)

    pgrid_ymin = max(pgrid_ymin, vgrid%ay%o)
    pgrid_ymax = min(pgrid_ymax, vgrid%ay%o + (vgrid%ay%n-1)*vgrid%ay%d)

    !!.. propagation grid 
    pgrid%ax%o = pgrid_xmin
    pgrid%ay%o = pgrid_ymin
    pgrid%ax%n = 1.5 + (pgrid_xmax - pgrid_xmin)/pgrid%ax%d
    pgrid%ay%n = 1.5 + (pgrid_ymax - pgrid_ymin)/pgrid%ay%d

    !!.. the first and last index of surface x and y on the propagation grid (for output use)
    ggrid%index_surpro_xmin = 1.5 + (ggrid%surface_xmin - pgrid%ax%o) / pgrid%ax%d
    ggrid%index_surpro_xmax = 1.5 + (ggrid%surface_xmax - pgrid%ax%o) / pgrid%ax%d

    ggrid%index_surpro_ymin = 1.5 + (ggrid%surface_ymin - pgrid%ay%o) / pgrid%ay%d
    ggrid%index_surpro_ymax = 1.5 + (ggrid%surface_ymax - pgrid%ay%o) / pgrid%ay%d

    return
  end subroutine

!=====================================================================================
  subroutine wei_greens_velocity_grid()
    real  :: vgrid_xmin, vgrid_xmax
    real  :: vgrid_ymin, vgrid_ymax
    real  :: pgrid_xmin, pgrid_xmax
    real  :: pgrid_ymin, pgrid_ymax

    vgrid_xmin = vgrid%ax%o
    vgrid_ymin = vgrid%ay%o
    vgrid_xmax = vgrid%ax%o + (vgrid%ax%n-1)*vgrid%ax%d
    vgrid_ymax = vgrid%ay%o + (vgrid%ay%n-1)*vgrid%ay%d

    pgrid_xmin = pgrid%ax%o
    pgrid_ymin = pgrid%ay%o
    pgrid_xmax = pgrid%ax%o + (pgrid%ax%n-1)*pgrid%ax%d
    pgrid_ymax = pgrid%ay%o + (pgrid%ay%n-1)*pgrid%ay%d

    if(pgrid_xmin < vgrid_xmin) then
       write(0,*)"pgrid_xmin, vgrid_xmin", pgrid_xmin, vgrid_xmin
       call seperr("propagation grid out of velocity grid: pgrid_xmin < vgrid_xmin")
    end if
    if(pgrid_xmax > vgrid_xmax) then
       write(0,*)"pgrid_xmax, vgrid_xmax", pgrid_xmax, vgrid_xmax
       call seperr("propagation grid out of velocity grid: pgrid_xmax > vgrid_xmax")
    end if
    if(pgrid_ymin < vgrid_ymin) then
       write(0,*)"pgrid_ymin, vgrid_ymin", pgrid_ymin, vgrid_ymin
       call seperr("propagation grid out of velocity grid: pgrid_ymin < vgrid_ymin")
    end if
    if(pgrid_ymax > vgrid_ymax) then
       write(0,*)"pgrid_ymax, vgrid_ymax", pgrid_ymax, vgrid_ymax
       call seperr("propagation grid out of velocity grid: pgrid_ymax > vgrid_ymax")
    end if

    !!.. computing the first and the last x,y index of the input velocity on the whole velocity grid
    vgrid%index_provel_xmin = 1.5 + (pgrid%ax%o - vgrid%ax%o)/vgrid%ax%d
    vgrid%index_provel_ymin = 1.5 + (pgrid%ay%o - vgrid%ay%o)/vgrid%ay%d
    vgrid%index_provel_xmax = vgrid%index_provel_xmin + vgrid%ax%n - 1
    vgrid%index_provel_ymax = vgrid%index_provel_ymin + vgrid%ay%n - 1
  end subroutine
!============================================================================================
  subroutine wei_greens_velmod_alloc(nzs_in)
    integer :: nzs_in
    allocate(vmod(pgrid%ax%n, pgrid%ay%n, nzs_in))
    allocate(vref(nref, nzs_in))
  end subroutine
!============================================================================================
  subroutine wei_greens_wfld_d_alloc(nws_in)
    integer  :: nws_in
    allocate(csou(nws_in))
    allocate(wfld_d(pgrid%ax%n, pgrid%ay%n, nws_in))
  end subroutine
!============================================================================================
  subroutine wei_greens_wfld_d_dealloc()
    deallocate(csou)
    deallocate(wfld_d)
  end subroutine
!============================================================================================
  subroutine wei_greens_velmod_dealloc()
    deallocate(vmod)
    deallocate(vref)
  end subroutine
!============================================================================================
  subroutine wei_greens_wfld_d_init(nws_in)
    integer  :: nws_in
    integer  :: iws, itr, i1, i2, ind1, ind2
    wfld_d = 0.
    !!.. locate the source wavefield onto the propagation grid 
    do iws=1, nws_in
       wfld_d(ggrid%index_modpro_modx, ggrid%index_modpro_mody, iws) = csou(iws)
    end do
  end subroutine
!============================================================================================
  subroutine wei_greens_io_init()
    complex, allocatable, dimension(:,:,:) :: ctmp
    integer                                :: i, err
    allocate(ctmp(ggrid%ax_surf%n, ggrid%ay_surf%n, ggrid%aw%n), stat=err)
    if (err /= 0) call seperr("trouble allocating ctmp in wei_greens_io_init()")
    do i=1, ggrid%ax_modl%n*ggrid%ay_modl%n*ggrid%az_modl%n
       call sep_write(ctmp,greens_, esize=8)
    end do
    deallocate(ctmp)
  end subroutine
!=============================================================================================
  subroutine wei_greens_io_csou(indw_in, nws_in)
    integer               :: indw_in, nws_in
    integer, dimension(1) :: n_cube, j_wind, n_wind, f_wind
    integer               :: n_dim_cube, err
    n_dim_cube = 1
    n_cube = ggrid%aw%n
    j_wind = 1
    n_wind = nws_in
    f_wind = indw_in - 1
    csou = 0.
    err = sreed_window(sou_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,csou)
    if (err.lt.0) call seperr("wei_greens_io_csou: trouble reading source wavefield")
  end subroutine
!==============================================================================================
  subroutine wei_greens_io_velo(indz_in, nzs_in)
    integer               :: indz_in
    integer               :: nzs_in
    integer, dimension(3) :: n_cube, j_wind, n_wind, f_wind
    integer               :: n_dim_cube, err, indz
    real, allocatable, dimension(:,:,:) :: tmp_vmod

    n_dim_cube = 3
    n_cube = (/vgrid%ax%n, vgrid%ay%n, vgrid%az%n /)
    j_wind = (/1, 1, 1/)
    n_wind = (/pgrid%ax%n, pgrid%ay%n, nzs_in/)
    f_wind = (/vgrid%index_provel_xmin-1, vgrid%index_provel_ymin-1, indz_in-1/)

    vmod = 0.
    err = sreed_window(vel_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,vmod)
    if (err .lt. 0) call seperr("sreed_window: trouble reading velocity file")
  end subroutine

!=============================================================================================
  subroutine wei_greens_io_wfld(flag, wfld, nws_in, ind_ws, ind_modx, ind_mody, ind_modz)
    character(len=4)          :: flag
    integer                   :: nws_in
    integer                   :: ind_ws, ind_modx, ind_mody, ind_modz
    integer, dimension(6)     :: n_cube, j_wind, n_wind, f_wind
    integer                   :: n_dim_cube, err
    integer                   :: indx_xmin, indx_xmax
    integer                   :: indx_ymin, indx_ymax
    complex, dimension(:,:,:), pointer :: wfld
    complex, dimension(:,:,:), pointer :: pwfld
    n_dim_cube = 6
    n_cube = (/ggrid%ax_surf%n, ggrid%ay_surf%n, ggrid%aw%n, ggrid%ax_modl%n, ggrid%ay_modl%n, ggrid%az_modl%n/)
    j_wind = (/1, 1, 1, 1, 1, 1/)
    n_wind = (/ggrid%ax_surf%n, ggrid%ay_surf%n, nws_in, 1, 1, 1/)
    f_wind = (/0, 0, ind_ws-1, ind_modx-1, ind_mody-1, ind_modz-1/)


    indx_xmin = ggrid%index_surpro_xmin; indx_xmax = ggrid%index_surpro_xmax
    indx_ymin = ggrid%index_surpro_ymin; indx_ymax = ggrid%index_surpro_ymax

    pwfld => wfld(indx_xmin:indx_xmax, indx_ymin:indx_ymax, :)
    
    if (flag == "rite") then
       err = srite_window(greens_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,pwfld)
    else if(flag == "reed") then
       pwfld = 0.
       err = sreed_window(greens_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,pwfld)
    else
       call seperr("you have to specify the flag for wei_io_greens()")
    end if
    nullify(pwfld)
  end subroutine
!================================================================================================================
  subroutine wei_greens_block_calc(xblock, ax)
    type(block)  :: xblock
    type(axis)   :: ax
    integer      :: ixb

    xblock%bres = .false.
    xblock%ds = ax%d
    xblock%nb = ax%n / xblock%ns
    xblock%nb_res = mod(ax%n, xblock%ns)

    if(xblock%nb_res /= 0) then
       xblock%nb = xblock%nb + 1
       xblock%bres = .true.
    end if

    if(allocated(xblock%blen)) deallocate(xblock%blen)
    if(allocated(xblock%os))   deallocate(xblock%os)
    if(allocated(xblock%indx)) deallocate(xblock%indx)

    allocate(xblock%blen(xblock%nb), xblock%os(xblock%nb), xblock%indx(xblock%nb))

    xblock%blen = xblock%ns

    if(xblock%bres) xblock%blen(xblock%nb) = xblock%nb_res

    do ixb=1, xblock%nb
       xblock%os(ixb)   = ax%o + (ixb-1)*xblock%ns*ax%d
       xblock%indx(ixb) = 1 + (ixb-1)*xblock%ns
    end do
  end subroutine
!==============================================================================================
  subroutine wei_greens_radians_calc()
  !!.. convert frequencies to radians
  if (.not.radian) then
     ggrid%aw%d  = ggrid%aw%d * 2.0 * PI
     ggrid%aw%o  = ggrid%aw%o * 2.0 * PI
  end if
  end subroutine
!===========================================================================
  subroutine wei_greens_k2_map()

    if(associated(kx2)) deallocate(kx2)
    if(associated(ky2)) deallocate(ky2)

    allocate(kx2(pgrid%ax%n))
    allocate(ky2(pgrid%ay%n))

    call k2_map(kx2, pgrid%ax%n, pgrid%ax%d)
    call k2_map(ky2, pgrid%ay%n, pgrid%ay%d)

  end subroutine
!============================================================================
!
  subroutine k2_map(k,n,d)
    real,  dimension(:) :: k
    integer, intent(in) :: n
    real,    intent(in) :: d
    integer             :: i,m
    k=0.
    if(n/=1) then
       m = n/2
       do i=1,m+1
          k(i)=i-1
       end do
       do i=m+2,n
          k(i)=i-n-1
       end do
    else
       k=0. !! make sure k=0 if n=1
    end if

    k = k / (n*d)      !! 1/length
    k = 2 * PI * k     !! radians
    k = k * k          !! k squared
  end subroutine

end module
