module grn_glocal_mod
  use  base_struct_mod
  implicit none

contains

  subroutine grn_wfs_init(wfds, csou, ggrid, pgrid)
    complex, dimension(:,:,:,:) :: wfds
    complex, dimension(:)       :: csou
    type(gre_grid)              :: ggrid
    type(pro_grid)              :: pgrid

    integer                     :: ix, iy
    real                        :: wtx, wty, PI
    PI = 3.14159265359

    wfds = 0.

    do iy=max(ggrid%index_surpro_sy-2,1), min(ggrid%index_surpro_sy+2,pgrid%ay%n); wty=0.5+0.5*cos(PI*real(iy-ggrid%index_surpro_sy)/3)
    do ix=max(ggrid%index_surpro_sx-2,1), min(ggrid%index_surpro_sx+2,pgrid%ax%n); wtx=0.5+0.5*cos(PI*real(ix-ggrid%index_surpro_sx)/3)
       wfds(ix,iy,1,:) = wtx*wty*csou(:)
    end do
    end do
    return
  end subroutine

  subroutine grn_progrd_gen3d(pgrid, vgrid, ggrid)
    type(pro_grid) :: pgrid
    type(vel_grid) :: vgrid
    type(gre_grid) :: ggrid 
    real           :: pgrid_xmin, pgrid_xmax
    real           :: pgrid_ymin, pgrid_ymax
    real           :: pgrid_zmin, pgrid_zmax

    pgrid%ax  = vgrid%ax
    pgrid%ay  = vgrid%ay
    pgrid%az  = vgrid%az

    pgrid%gxmin = ggrid%surf_xmin
    pgrid%gxmax = ggrid%surf_xmax
    pgrid%gymin = ggrid%surf_ymin
    pgrid%gymax = ggrid%surf_ymax

    pgrid_xmin = min(ggrid%surf_xmin,ggrid%modl_xmin); pgrid_xmax = max(ggrid%surf_xmax,ggrid%modl_xmax)
    pgrid_ymin = min(ggrid%surf_ymin,ggrid%modl_ymax); pgrid_ymax = max(ggrid%surf_ymax,ggrid%modl_ymax)
    pgrid_zmin = vgrid%az%o                          ; pgrid_zmax = ggrid%modl_zmax


    write(0,*)"pgrid_x", pgrid_xmin, pgrid_xmax, ggrid%surf_xmin, ggrid%surf_xmax, ggrid%modl_xmin, ggrid%modl_xmax
    write(0,*)"pgrid_y", pgrid_ymin, pgrid_ymax
    write(0,*)"pgrid_z", pgrid_zmin, pgrid_zmax

    !!.. 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 
    if (pgrid_xmin<vgrid%ax%o .or. pgrid_xmax>vgrid%ax%o + (vgrid%ax%n-1)*vgrid%ax%d) then
       write(0,*)"ERROR: pgrid_xmin < vgrid_xmin .or. pgrid_max > vgrid_xmax"
       stop
    end if
    if (pgrid_ymin<vgrid%ay%o .or. pgrid_ymax>vgrid%ay%o + (vgrid%ay%n-1)*vgrid%ay%d) then
       write(0,*)"ERROR: pgrid_ymin < vgrid_ymin .or. pgrid_max > vgrid_ymax"
       stop
    end if

    !!.. propagation grid 
    pgrid%ax%o = pgrid_xmin; pgrid%ax%n = 1.5+(pgrid_xmax-pgrid_xmin)/pgrid%ax%d
    pgrid%ay%o = pgrid_ymin; pgrid%ay%n = 1.5+(pgrid_ymax-pgrid_ymin)/pgrid%ay%d
    pgrid%az%o = pgrid_zmin; pgrid%az%n = 1.5+(pgrid_zmax-pgrid_zmin)/pgrid%az%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%surf_xmin - pgrid%ax%o) / pgrid%ax%d
    ggrid%index_surpro_xmax = 1.5 + (ggrid%surf_xmax - pgrid%ax%o) / pgrid%ax%d

    ggrid%index_surpro_ymin = 1.5 + (ggrid%surf_ymin - pgrid%ay%o) / pgrid%ay%d
    ggrid%index_surpro_ymax = 1.5 + (ggrid%surf_ymax - pgrid%ay%o) / pgrid%ay%d

    !!.. the first and last index of model x, y, and z on the pro grid
    ggrid%index_modpro_xmin = 1.5 + (ggrid%modl_xmin - pgrid%ax%o) / pgrid%ax%d
    ggrid%index_modpro_xmax = 1.5 + (ggrid%modl_xmax - pgrid%ax%o) / pgrid%ax%d

    ggrid%index_modpro_ymin = 1.5 + (ggrid%modl_ymin - pgrid%ay%o) / pgrid%ay%d
    ggrid%index_modpro_ymax = 1.5 + (ggrid%modl_ymax - pgrid%ay%o) / pgrid%ay%d

    ggrid%index_modpro_zmin = 1.5 + (ggrid%modl_zmin - pgrid%az%o) / pgrid%az%d
    ggrid%index_modpro_zmax = 1.5 + (ggrid%modl_zmax - pgrid%az%o) / pgrid%az%d
    return
  end subroutine

end module
