module wei_hessian_off2ang2d_mod2
  use wei_type_mod
  use wei_off2ang2d_sd_mod
  use sep
  implicit none
  real, dimension(:,:,:), pointer, private :: hessian_off
  real, dimension(:,:,:), pointer, private :: hessian_ang

  real, dimension(:,:)  , pointer, private :: phessian_off
  real, dimension(:,:)  , pointer, private :: phessian_ang

  type(hes_grid)                 , private :: hgrid
  integer                        , private :: style
  integer                        , private :: step
  
  character(len=256)             , private :: hessian_off_="hessian_off"
  character(len=256)             , private :: hessian_ang_="hessian_ang"
  logical, parameter             , private :: T=.true., F=.false.
  integer, external              , private :: sreed_window, srite_window
contains

  integer function wei_hessian_off2ang2d_init() result(stat)

    call wei_hessian_off2ang2d_input_history()
    call wei_hessian_off2ang2d_param_init()
    call wei_hessian_off2ang2d_output_history()
    call wei_hessian_off2ang2d_array_alloc()
    if (step == 1) then
       write(0,*)"==============================================="
       write(0,*)"         2-D HESSIAN OFF2ANG TRANSFORM         "
       write(0,*)"         STEP 1: off-diagonal dimension        "
       write(0,*)"-----------------------------------------------"
    else
       write(0,*)"==============================================="
       write(0,*)"         2-D HESSIAN OFF2ANG TRANSFORM         " 
       write(0,*)"         STEP 2: main-diagonal dimension       "
       write(0,*)"-----------------------------------------------"
    end if
    stat = 0
  end function
!=========================================================================================================
  integer function wei_hessian_off2ang2d_clos() result(stat)
    call wei_hessian_off2ang2d_array_dealloc()
    stat = 0
  end function
!=========================================================================================================
  integer function wei_hessian_off2ang2d_oper() result(stat)
    integer    :: ix_modl, ihx_modl, iz_modl, ix_offd
    real       :: x_modl, hx_modl, z_modl
    real       :: oz_offd, ohx_offd
    type(axis) :: az_offd, ahx_offd, aax_offd 

    aax_offd = hgrid%aax_offd

    do iz_modl=1, hgrid%az_modl%n
       oz_offd = hgrid%az_offd%o
       if (step == 1) then
          z_modl  = hgrid%az_modl%o + (iz_modl-1)*hgrid%az_modl%d
          oz_offd = z_modl + hgrid%az_offd%o
       end if

       az_offd%n = hgrid%az_offd%n; az_offd%o = oz_offd; az_offd%d = hgrid%az_offd%d
write(0,*)"iz_modl", iz_modl, oz_offd

       do ihx_modl=1, hgrid%ahx_modl%n
          ohx_offd = hgrid%ahx_offd%o
          if (step == 1) then
             hx_modl  = hgrid%ahx_modl%o + (ihx_modl-1)*hgrid%ahx_modl%d
             ohx_offd = hx_modl + hgrid%ahx_offd%o
          end if

          ahx_offd%n = hgrid%ahx_offd%n; ahx_offd%o = ohx_offd; ahx_offd%d = hgrid%ahx_offd%d

          stat = wei_off2ang2d_sd_init(az_offd , ahx_offd, aax_offd, style)

          do ix_modl=1, hgrid%ax_modl%n
             call wei_hessian_off2ang2d_io_angle("reed", hessian_off_, hessian_off, ix_modl, iz_modl, ihx_modl)
             do ix_offd=1, hgrid%ax_offd%n
                phessian_off => hessian_off(:,:,ix_offd)
                phessian_ang => hessian_ang(:,:,ix_offd)
                stat = wei_off2ang2d_sd_oper(T, F, phessian_ang, phessian_off)
             end do
             call wei_hessian_off2ang2d_io_angle("rite", hessian_ang_, hessian_ang, ix_modl, iz_modl, ihx_modl)
          end do
          stat = wei_off2ang2d_sd_clos()
       end do  
    end do
    stat = 0
  end function 
!============================================================================================================
  subroutine wei_hessian_off2ang2d_io_angle(flag, filetag, hessian, index_x_modl, index_z_modl, index_hx_modl)
    character(len=4)                    :: flag
    character(len=256)                  :: filetag
    integer                             :: index_x_modl, index_z_modl, index_hx_modl
    real, dimension(:,:,:), pointer     :: hessian
    integer, dimension(6)               :: n_cube, j_wind, n_wind, f_wind
    integer                             :: n_dim_cube, err

    n_dim_cube = 6
    !n_cube = (/hgrid%az_offd%n, hgrid%ahx_offd%n, hgrid%ax_offd%n, &
    !           hgrid%az_modl%n, hgrid%ahx_modl%n, hgrid%ax_modl%n/)
    j_wind = (/1, 1, 1, 1, 1, 1/)
    !n_wind = (/hgrid%az_offd%n, hgrid%ahx_offd%n, hgrid%ax_offd%n, 1, 1, 1/)
    f_wind = (/0, 0, 0, index_z_modl-1, index_hx_modl-1, index_x_modl-1/)
    if (flag == "rite") then
       n_cube = (/hgrid%az_offd%n, hgrid%aax_offd%n, hgrid%ax_offd%n, &
                  hgrid%az_modl%n, hgrid%ahx_modl%n, hgrid%ax_modl%n/)
       n_wind = (/hgrid%az_offd%n, hgrid%aax_offd%n, hgrid%ax_offd%n, 1, 1, 1/)
       err = srite_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,hessian)

    else if (flag == "reed") then

       n_cube = (/hgrid%az_offd%n, hgrid%ahx_offd%n, hgrid%ax_offd%n, &
                  hgrid%az_modl%n, hgrid%ahx_modl%n, hgrid%ax_modl%n/)
       n_wind = (/hgrid%az_offd%n, hgrid%ahx_offd%n, hgrid%ax_offd%n, 1, 1, 1/)
       !f_wind = (/0, 0, 0, index_z_modl-1, index_hx_modl-1, index_x_modl-1/)
       hessian = 0.
       err = sreed_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,hessian)

    else
       call seperr("you have to specify the flag for wei_hessian_io_hessian()")

    end if
  end subroutine
!======================================================================
  subroutine wei_hessian_off2ang2d_param_init()
    call from_param("nang", hgrid%aax_offd%n, 31)
    call from_param("oang", hgrid%aax_offd%o, -30.)
    call from_param("dang", hgrid%aax_offd%d, 2.)
    call from_param("style", style,  1 )
    call from_param("step", step,    1 )

  end subroutine
!===========================================================================
  subroutine wei_hessian_off2ang2d_array_alloc()
    allocate(hessian_off(hgrid%az_offd%n, hgrid%ahx_offd%n, hgrid%ax_offd%n))
    allocate(hessian_ang(hgrid%az_offd%n, hgrid%aax_offd%n, hgrid%ax_offd%n))
  end subroutine
!===========================================================================
  subroutine wei_hessian_off2ang2d_array_dealloc()
    deallocate(hessian_off)
    deallocate(hessian_ang)
  end subroutine
!===========================================================================
  subroutine wei_hessian_off2ang2d_input_history()

    call auxin(hessian_off_)

    call from_aux(hessian_off_, "n1", hgrid%az_offd%n)
    call from_aux(hessian_off_, "d1", hgrid%az_offd%d)
    call from_aux(hessian_off_, "o1", hgrid%az_offd%o)
   
    call from_aux(hessian_off_, "n2", hgrid%ahx_offd%n)
    call from_aux(hessian_off_, "d2", hgrid%ahx_offd%d)
    call from_aux(hessian_off_, "o2", hgrid%ahx_offd%o)

    call from_aux(hessian_off_, "n3", hgrid%ax_offd%n)
    call from_aux(hessian_off_, "d3", hgrid%ax_offd%d)
    call from_aux(hessian_off_, "o3", hgrid%ax_offd%o)

    call from_aux(hessian_off_, "n4", hgrid%az_modl%n)
    call from_aux(hessian_off_, "d4", hgrid%az_modl%d)
    call from_aux(hessian_off_, "o4", hgrid%az_modl%o)

    call from_aux(hessian_off_, "n5", hgrid%ahx_modl%n)
    call from_aux(hessian_off_, "d5", hgrid%ahx_modl%d)
    call from_aux(hessian_off_, "o5", hgrid%ahx_modl%o)

    call from_aux(hessian_off_, "n6", hgrid%ax_modl%n)
    call from_aux(hessian_off_, "d6", hgrid%ax_modl%d)
    call from_aux(hessian_off_, "o6", hgrid%ax_modl%o)

    
  end subroutine
!===================================================================
  subroutine wei_hessian_off2ang2d_output_history()

    call auxout(hessian_ang_)

    call auxputch("n1", "i", hgrid%az_offd%n, hessian_ang_)
    call auxputch("d1", "r", hgrid%az_offd%d, hessian_ang_)
    call auxputch("o1", "r", hgrid%az_offd%o, hessian_ang_)

    call auxputch("n2", "i", hgrid%aax_offd%n, hessian_ang_)
    call auxputch("d2", "r", hgrid%aax_offd%d, hessian_ang_)
    call auxputch("o2", "r", hgrid%aax_offd%o, hessian_ang_)

    call auxputch("n3", "i", hgrid%ax_offd%n, hessian_ang_)
    call auxputch("d3", "r", hgrid%ax_offd%d, hessian_ang_)
    call auxputch("o3", "r", hgrid%ax_offd%o, hessian_ang_)

    call auxputch("n4", "i", hgrid%az_modl%n, hessian_ang_)
    call auxputch("d4", "r", hgrid%az_modl%d, hessian_ang_)
    call auxputch("o4", "r", hgrid%az_modl%o, hessian_ang_)

    call auxputch("n5", "i", hgrid%ahx_modl%n, hessian_ang_)
    call auxputch("d5", "r", hgrid%ahx_modl%d, hessian_ang_)
    call auxputch("o5", "r", hgrid%ahx_modl%o, hessian_ang_)

    call auxputch("n6", "i", hgrid%ax_modl%n, hessian_ang_)
    call auxputch("d6", "r", hgrid%ax_modl%d, hessian_ang_)
    call auxputch("o6", "r", hgrid%ax_modl%o, hessian_ang_)
    if (step == 1) then
       call auxputch("label1", "s", "z_offd", hessian_ang_)
       call auxputch("label2", "s", "ax_offd", hessian_ang_)
       call auxputch("label3", "s", "x_offd", hessian_ang_)
       call auxputch("label4", "s", "z_modl", hessian_ang_)
       call auxputch("label5", "s", "hx_modl", hessian_ang_)
       call auxputch("label6", "s", "x_modl", hessian_ang_) 
    else
       call auxputch("label1", "s", "z_modl", hessian_ang_)
       call auxputch("label2", "s", "ax_modl", hessian_ang_)
       call auxputch("label3", "s", "x_modl", hessian_ang_)
       call auxputch("label4", "s", "z_offd", hessian_ang_)
       call auxputch("label5", "s", "ax_offd", hessian_ang_)
       call auxputch("label6", "s", "x_offd", hessian_ang_)
    end if 
  end subroutine  


end module








