module wei_hessian_off2ang2d_mod
  use wei_type_mod
  use 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
  real                           , private :: axmin_modl, axmax_modl
  real                           , private :: axmin_offd, axmax_offd
  integer                        , private :: style
  logical                        , private :: offd
  
  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()
    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)




  subroutine wei_hessian_off2ang2d_step_one()
    integer  :: ix_modl, ihx_modl, iz_modl
    integer  :: ix_offd, ihx_offd, iz_offd
    integer  :: nx_modl, nhx_modl, nz_modl
    integer  :: nx_offd, nhx_offd, nz_offd
    real     :: x_modl, hx_modl, z_modl, oz_offd, ohx_offd

!======== step 1 =======================================
! off2ang transform in the off-diagonal dimensions
!-------------------------------------------------------
write(0,*)"hessian_off2ang: STEP 1"
    nx_modl  = 1; nz_modl = 1; nhx_modl = 1
    nx_offd  = hgrid%ax_offd%n
    nz_offd  = hgrid%az_offd%n
    nhx_offd = hgrid%ahx_offd%n

    call wei_hessian_off2ang2d_array_alloc(nz_offd, nhx_offd, nx_offd)

    do iz_modl=1, hgrid%az_modl%n
       z_modl = hgrid%az_modl%o + (iz_modl-1)*hgrid%az_modl%d
       oz_offd = z_modl + hgrid%az_offd%o
write(0,*)"iz_modl=", iz_modl

       do ihx_modl=1, hgrid%ahx_modl%n
          hx_modl = hgrid%ahx_modl%o + (ihx_modl-1)*hgrid%ahx_modl%d
          ohx_offd = hx_modl + hgrid%ahx_offd%o
          stat = off2ang2d_sd_init(oz_offd , ohx_offd, hgrid%az_offd%d, hgrid%ahx_offd%d, &
                                   hgrid%az_offd%n, hgrid%ahx_offd%n, axmin_offd, axmax_offd, style)

          do ix_modl=1, hgrid%ax_modl%n
             x_modl = hgrid%ax_modl%o + (ix_modl-1)*hgrid%ax_modl%d
!write(0,*)"iz_modl, ihx_modl, ix_modl", iz_modl, ihx_modl, ix_modl
             call wei_hessian_off2ang2d_io_angle("reed", hessian_off_, hessian_off, &
                                                  nx_modl, nz_modl, nhx_modl,       &
                                                  nx_offd, nz_offd, nhx_offd,       &
                                                  ix_modl, iz_modl, ihx_modl,       &
                                                        1,       1,        1)

             do ix_offd=1, hgrid%ax_offd%n
                phessian_off => hessian_off(:,:,ix_offd)
                phessian_ang => hessian_ang(:,:,ix_offd)
                stat = off2ang2d_sd_oper(T, F, phessian_ang, phessian_off)
             end do
             call wei_hessian_off2ang2d_io_angle("rite", hessian_ang_, hessian_ang, &
                                                  nx_modl, nz_modl, nhx_modl,       &
                                                  nx_offd, nz_offd, nhx_offd,       & 
                                                  ix_modl, iz_modl, ihx_modl,       &
                                                        1,       1,        1)
          end do
       end do  
    end do
    call wei_hessian_off2ang2d_array_dealloc()

!============ step 2 ===================================
! off2ang transform in the model dimension
!-------------------------------------------------------
write(0,*)"hessian_off2ang: STEP 2"
    nx_modl  = hgrid%ax_modl%n
    nz_modl  = hgrid%az_modl%n
    nhx_modl = hgrid%ahx_modl%n
    nx_offd  = 1; nz_offd = 1; nhx_offd = 1

    call wei_hessian_off2ang2d_array_alloc(nz_modl, nhx_modl, nx_modl)

    stat = off2ang2d_sd_init(hgrid%az_modl%o, hgrid%ahx_modl%o, hgrid%az_modl%d, hgrid%ahx_modl%d, &
                             hgrid%az_modl%n, hgrid%ahx_modl%n, axmin_modl, axmax_modl, style)

    do iz_offd =1, hgrid%az_offd%n
    do ihx_offd=1, hgrid%ahx_offd%n
    do ix_offd =1, hgrid%ax_offd%n
write(0,*)"iz_offd,ihx_offd,ix_offd", iz_offd,ihx_offd,ix_offd
!       call wei_hessian_off2ang2d_io_angle("reed", hessian_ang_, hessian_off, &
!                                            nx_modl, nz_modl, nhx_modl,       &
!                                            nx_offd, nz_offd, nhx_offd,       & 
!                                                  1,       1,        1,       &
!                                            ix_offd, iz_offd, ihx_offd)
!write(0,*)"111111111111"
        do ix_modl=1, hgrid%ax_modl%n
           phessian_off => hessian_off(:,:,ix_modl)
           phessian_ang => hessian_ang(:,:,ix_modl)
           !stat = off2ang2d_sd_oper(T, F, phessian_ang, phessian_off)
        end do
!write(0,*)"222222222222222"
!        call wei_hessian_off2ang2d_io_angle("rite", hessian_ang_, hessian_ang, &
!                                             nx_modl, nz_modl, nhx_modl,       &
!                                             nx_offd, nz_offd, nhx_offd,       &
!                                                   1,       1,        1,       &
!                                             ix_offd, iz_offd, ihx_offd)
!write(0,*)"3333333333333"
    end do
    end do
    end do
    stat = 0
  end function 
!============================================================================================================
  subroutine wei_hessian_off2ang2d_io_step_one(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
       err = srite_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,hessian)

    else if (flag == "reed") then
       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("amin_modl", axmin_modl, -60.)
    call from_param("amax_modl", axmax_modl,  60.)
    call from_param("amin_offd", axmin_offd, -60.)
    call from_param("amax_offd", axmax_offd,  60.)
    call from_param("style", style,  1 )
    call from_param("offd", offd, .true.)
    hgrid%aax_offd%n = hgrid%ahx_offd%n
    hgrid%aax_offd%o = axmin_offd
    if (hgrid%aax_offd%n == 1) then
      hgrid%aax_offd%d = 0
    else
      hgrid%aax_offd%d = (axmax_offd - axmin_offd) / (hgrid%aax_offd%n - 1)
    end if
    
    hgrid%aax_modl%n = hgrid%ahx_modl%n
    hgrid%aax_modl%o = axmin_modl
    if (hgrid%aax_modl%n == 1) then
       hgrid%aax_modl%d = 0
    else
       hgrid%aax_modl%d = (axmax_modl - axmin_modl) / (hgrid%aax_modl%n - 1)
    end if 
  end subroutine
!===========================================================================
  subroutine wei_hessian_off2ang2d_array_alloc(nz, nhx, nx)
    integer  :: nz, nhx, nx
    allocate(hessian_off(nz, nhx, nx))
    allocate(hessian_ang(nz, nhx, nx))
  end subroutine
!===========================================================================
  subroutine wei_hessian_off2ang2d_array_dealloc()
    deallocate(hessian_off)
    deallocate(hessian_ang)
  end subroutine
!===========================================================================
  subroutine wei_hessian_off2ang2d_input_history()

    call auxinout(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 auxinout(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%aax_modl%n, hessian_ang_)
    call auxputch("d5", "r", hgrid%aax_modl%d, hessian_ang_)
    call auxputch("o5", "r", hgrid%aax_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_)

    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", "ax_modl", hessian_ang_)
    call auxputch("label6", "s", "x_modl", hessian_ang_) 

  end subroutine  


end module








