program dyn_matrix_test
  use dyn_matrix, only: DYNMTX, Create, Close, PutLine, Solve, initialized
  use constants, only: CONST, Create, Get
  use dc_test, only: AssertEqual, AssertGreaterThan, AssertLessThan
  use dc_args, only: ARGS, Open, HelpMsg, Option, Debug, Help, Strict, Close
  use dc_types, only: DP, STRING
  use dc_string, only: toChar
  use dc_message, only: MessageNotify
  use gt4_history, only: GT_HISTORY, HistoryGet
  implicit none

  !---------------------------------------------------------
  !  実験の表題, モデルの名称, 所属機関名
  !  Title of a experiment, name of model, sub-organ
  !---------------------------------------------------------
  character(*), parameter:: title = 'dyn_matrix_test $Name: dcpam4-20070731 $ :: ' // 'Test program of "dyn_matrix" module'
  character(*), parameter:: source = 'dcpam4 ' // '(See http://www.gfd-dennou.org/library/dcpam)'
  character(*), parameter:: institution = 'GFD Dennou Club (See http://www.gfd-dennou.org)'

  !---------------------------------------------------------
  !  作業変数
  !  Work variables
  !---------------------------------------------------------
  type(ARGS):: arg            ! コマンドライン引数. 
                              ! Command line arguments
!!$  logical:: OPT_namelist      ! -N, --namelist オプションの有無. 
!!$                              ! Existence of '-N', '--namelist' option
!!$  character(STRING):: VAL_namelist
!!$                              ! -N, --namelist オプションの値. 
!!$                              ! Value of '-N', '--namelist' option
  type(CONST):: const_earth
  real(DP):: RPlanet    ! $ a $ .      惑星半径.       Radius of planet
  type(DYNMTX):: dyn_mtx00, dyn_mtx01
  real(DP):: DelTime
  integer, allocatable:: nmo(:,:,:)
  real(DP), allocatable:: wz_DiffVorDiv (:,:)
  real(DP), allocatable:: wz_DiffTherm (:,:)
  real(DP), allocatable:: zz_siMtxWH(:,:)
  real(DP), allocatable:: zz_siMtxGCt(:,:)
  real(DP), allocatable:: wzz_siMtxLU(:,:,:)
  integer, allocatable:: wz_siMtxPiv(:,:)
  real(DP), allocatable:: wz_DivAvrTime(:,:)
  real(DP), allocatable:: wz_DivAvrTimeAns(:,:)
  real(DP), allocatable:: wz_DivAvrTimeAnsL(:,:)
  real(DP), allocatable:: wz_DivAvrTimeAnsG(:,:)
  integer:: m, l, k, kk
  real:: r00
  real(DP):: d00
  real(DP), allocatable:: d20(:,:)
  real(DP), allocatable:: d10(:)
  real, allocatable:: r10(:)
  logical:: err
  character(*), parameter:: subname = 'dyn_matrix_test'
  character(*), parameter:: spdata_file = 'dyn_matrix_test00.nc'
  character(*), parameter:: create_answer_file = 'dyn_matrix_test01.nc'
  character(*), parameter:: solve_src_file = 'dyn_matrix_test02.nc'
  character(*), parameter:: solve_ans_file = 'dyn_matrix_test03.nc'
continue

  !---------------------------------------------------------
  !  コマンドライン引数の処理
  !  Command line arguments handling
  !---------------------------------------------------------
  call Open( arg )
  call HelpMsg( arg, 'Title', title )
  call HelpMsg( arg, 'Usage', './dyn_matrix_test [Options]' )
  call HelpMsg( arg, 'Source', source )
  call HelpMsg( arg, 'Institution', institution )
!!$  call Option( arg, StoA('-N', '--namelist'), &
!!$    & OPT_namelist, VAL_namelist, help = "NAMELIST filename" )
  call Debug( arg ) 
   call Help( arg ) 
   call Strict( arg, severe = .true. )
  call Close( arg )

  !---------------------------------------------------------
  !  物理定数の準備
  !  Prepare physical constants
  !---------------------------------------------------------
  call Create( const_earth ) ! (inout)

  DelTime = 1200.0_DP

  call Get( constant = const_earth, RPlanet = RPlanet )             ! (out)


  !---------------------------------------------------------
  !  スペクトルの添字順番の準備
  !  Prepare spectral subscript expression
  !---------------------------------------------------------
  call MessageNotify('M', subname, 'Input ' // spdata_file // '  ...' )

  if ( allocated(nmo) ) deallocate(nmo)
  allocate( nmo(1:2, 0:10, 0:10) )
  if ( allocated(wz_DiffVorDiv) ) deallocate(wz_DiffVorDiv)
  allocate( wz_DiffVorDiv(11**2, 0:1) )
  if ( allocated(wz_DiffTherm) ) deallocate(wz_DiffTherm)
  allocate( wz_DiffTherm(11**2, 0:1) )

  do l = 0, 10
    do m = 0, 10
      call HistoryGet(spdata_file, 'nmo1', r00, range='m='//trim(toChar(m))//',l='//trim(toChar(l)), quiet=.true.)
      nmo(1,m,l) = int(r00)
    end do
  end do

  do l = 0, 10
    do m = 0, 10
      call HistoryGet(spdata_file, 'nmo2', r00, range='m='//trim(toChar(m))//',l='//trim(toChar(l)), quiet=.true.)
      nmo(2,m,l) = int(r00)
    end do
  end do

!!$  do l = 0, 10
!!$    do m = 0, 10
!!$      write(*,*) 'nmo(1,', m, ',', l, ')= ', nmo(1,m,l)
!!$    end do
!!$  end do
!!$
!!$  do l = 0, 10
!!$    do m = 0, 10
!!$      write(*,*) 'nmo(2,', m, ',', l, ')= ', nmo(2,m,l)
!!$    end do
!!$  end do

  !---------------------------------------------------------
  !  初期設定に必要なデータの準備
  !  Prepare data for initialization
  !---------------------------------------------------------
  wz_DiffVorDiv = 0.0_DP
  do l = 0, 10
    do m = 0, 10
      if (nmo(1,m,l) < 1) cycle
      call HistoryGet(spdata_file, 'DiffVorDiv1', d00, range='m='//trim(toChar(m))//',l='//trim(toChar(l)), quiet=.true.)
      wz_DiffVorDiv(nmo(1,m,l),:) = d00
    end do
  end do

  do l = 0, 10
    do m = 0, 10
      if (nmo(2,m,l) < 1) cycle
      call HistoryGet(spdata_file, 'DiffVorDiv2', d00, range='m='//trim(toChar(m))//',l='//trim(toChar(l)), quiet=.true.)
      wz_DiffVorDiv(nmo(2,m,l),:) = d00
    enddo
  enddo

  wz_DiffTherm = 0.0_DP
  do l = 0, 10
    do m = 0, 10
      if (nmo(1,m,l) < 1) cycle
      call HistoryGet(spdata_file, 'DiffTherm1', d00, range='m='//trim(toChar(m))//',l='//trim(toChar(l)), quiet=.true.)
      wz_DiffTherm(nmo(1,m,l),:) = d00
    enddo
  enddo

  do l = 0, 10
    do m = 0, 10
      if (nmo(2,m,l) < 1) cycle
      call HistoryGet(spdata_file, 'DiffTherm2', d00, range='m='//trim(toChar(m))//',l='//trim(toChar(l)), quiet=.true.)
      wz_DiffTherm(nmo(2,m,l),:) = d00
    enddo
  enddo

!!$  do l = 0, 10
!!$    do m = 0, 10
!!$      if (m+l > 10) cycle
!!$      write(*,*) 'wz_DiffVorDiv(', 'nmo(1,', m, ',', l, '),0)= ', &
!!$        & wz_DiffVorDiv(nmo(1,m,l),0)
!!$    end do
!!$  end do
!!$
!!$  do l = 0, 10
!!$    do m = 0, 10
!!$      if (m+l > 10) cycle
!!$      write(*,*) 'wz_DiffVorDiv(', 'nmo(2,', m, ',', l, '),0)= ', &
!!$        & wz_DiffVorDiv(nmo(2,m,l),0)
!!$    end do
!!$  end do
!!$
!!$  do l = 0, 10
!!$    do m = 0, 10
!!$      if (m+l > 10) cycle
!!$      write(*,*) 'wz_DiffTherm(', 'nmo(1,', m, ',', l, '),0)= ', &
!!$        & wz_DiffTherm(nmo(1,m,l),0)
!!$    end do
!!$  end do
!!$
!!$  do l = 0, 10
!!$    do m = 0, 10
!!$      if (m+l > 10) cycle
!!$      write(*,*) 'wz_DiffTherm(', 'nmo(2,', m, ',', l, '),0)= ', &
!!$        & wz_DiffTherm(nmo(2,m,l),0)
!!$    end do
!!$  end do

  allocate( zz_siMtxWH(0:1, 0:1) )
  allocate( zz_siMtxGCt(0:1, 0:1) )
  zz_siMtxWH = reshape((/   0.000000000000_DP, 7741.623281371381_DP, 7741.623281371375_DP, 33725.56984179635_DP/), (/2,2/))
  zz_siMtxGCt = reshape((/43056.00000000003_DP, 43056.00000000001_DP, 43056.00000000003_DP, 43056.00000000001_DP/), (/2,2/))


  !---------------------------------------------------------
  !  初期設定テスト
  !  Initializationテスト
  !---------------------------------------------------------
  call Create( dyn_mtx01, nmax = -10, kmax = 2, nmo = nmo, RPlanet = RPlanet, DelTime = DelTime, zz_siMtxWH = zz_siMtxWH, zz_siMtxGCt = zz_siMtxGCt, w_DiffVorDiv = wz_DiffVorDiv(:,0), w_DiffTherm = wz_DiffTherm(:,0), err = err )                             ! (out)
  call AssertEqual('initialization test 1', answer = .true., check = err )
  call PutLine( dyn_mtx01 )

  call Create( dyn_mtx01, nmax = 10, kmax = 2, nmo = nmo, RPlanet = RPlanet, DelTime = DelTime, zz_siMtxWH = zz_siMtxWH, zz_siMtxGCt = zz_siMtxGCt, w_DiffVorDiv = wz_DiffVorDiv(:,0), w_DiffTherm = wz_DiffTherm(:,0) )       ! (in)
  call AssertEqual('initialization test 2', answer = .true., check = initialized(dyn_mtx01) )
  call PutLine( dyn_mtx01 )

  call Close( dyn_mtx01 )
  call AssertEqual('termination test 1', answer = .false., check = initialized(dyn_mtx01) )
  call PutLine( dyn_mtx01 )

  call Close( dyn_mtx01, err = err )
  call AssertEqual('termination test 2', answer = .true., check = err )
  call PutLine( dyn_mtx01 )

  call Create( dyn_mtx00, nmax = 10, kmax = 2, nmo = nmo, RPlanet = RPlanet, DelTime = DelTime, zz_siMtxWH = zz_siMtxWH, zz_siMtxGCt = zz_siMtxGCt, w_DiffVorDiv = wz_DiffVorDiv(:,0), w_DiffTherm = wz_DiffTherm(:,0) )       ! (in)
  call AssertEqual('initialization test 3', answer = .true., check = initialized(dyn_mtx00) )
  call PutLine(dyn_mtx00)


  !---------------------------------------------------------
  !  初期設定値のテスト
  !  Test initial setting values
  !---------------------------------------------------------
  call MessageNotify('M', subname, 'Input ' // create_answer_file // '  ...' )

  if ( allocated(wzz_siMtxLU) ) deallocate(wzz_siMtxLU)
  allocate( wzz_siMtxLU(11**2, 0:1, 0:1) )
  allocate( d20(2,2) )

  wzz_siMtxLU = 0.0_DP
  do l = 0, 10
    do m = 0, 10
      if (nmo(1,m,l) < 1) cycle
      call HistoryGet(create_answer_file, 'siMtxLU1', d20, range='m='//trim(toChar(m))//',l='//trim(toChar(l)), quiet=.true.)
      wzz_siMtxLU(nmo(1,m,l),:,:) = d20
    end do
  end do

  do l = 0, 10
    do m = 0, 10
      if (nmo(2,m,l) < 1) cycle
      call HistoryGet(create_answer_file, 'siMtxLU2', d20, range='m='//trim(toChar(m))//',l='//trim(toChar(l)), quiet=.true.)
      wzz_siMtxLU(nmo(2,m,l),:,:) = d20
    enddo
  end do

  call AssertLessThan( 'LU decomposition for semi-implicit matrix test 1', answer = (wzz_siMtxLU + 1.0e-10_DP) * 1.001_DP, check = dyn_mtx00 % wzz_siMtxLU)

  call AssertGreaterThan( 'LU decomposition for semi-implicit matrix test 1', answer = (wzz_siMtxLU - 1.0e-10_DP) * 0.999_DP, check = dyn_mtx00 % wzz_siMtxLU)

!!$  do l = 0, 10
!!$    do m = 0, 10
!!$      if (m+l > 10) cycle
!!$      do k = 0, 1
!!$        do kk = 0, 1
!!$          write(*,*) 'dyn_mtx00 % wzz_siMtxLU(', 'nmo(1,', m, ',', l, '),', k+1, ',', kk+1, ')= ', &
!!$            & dyn_mtx00 % wzz_siMtxLU(nmo(1,m,l),k,kk)
!!$          write(*,*) 'wzz_siMtxLU(', 'nmo(1,', m, ',', l, '),', k+1, ',', kk+1, ')= ', &
!!$            & wzz_siMtxLU(nmo(1,m,l),k,kk)
!!$        end do
!!$      end do
!!$    end do
!!$  end do
!!$
!!$  do l = 0, 10
!!$    do m = 0, 10
!!$      if (m+l > 10) cycle
!!$      do k = 0, 1
!!$        do kk = 0, 1
!!$          write(*,*) 'wzz_siMtxLU(', 'nmo(2,', m, ',', l, '),', k+1, ',', kk+1, ')= ', &
!!$            & dyn_mtx00 % wzz_siMtxLU(nmo(2,m,l),k,kk)
!!$        end do
!!$      end do
!!$    end do
!!$  end do

  if ( allocated(wz_siMtxPiv) ) deallocate(wz_siMtxPiv)
  allocate( wz_siMtxPiv(11**2, 0:1) )
  allocate( r10(2) )

  wz_siMtxPiv = 0.0_DP
  do l = 0, 10
    do m = 0, 10
      if (nmo(1,m,l) < 1) cycle
      call HistoryGet(create_answer_file, 'siMtxPiv1', r10, range='m='//trim(toChar(m))//',l='//trim(toChar(l)), quiet=.true.)
      wz_siMtxPiv(nmo(1,m,l),:) = int(r10)
    end do
  end do

  do l = 0, 10
    do m = 0, 10
      if (nmo(2,m,l) < 1) cycle
      call HistoryGet(create_answer_file, 'siMtxPiv2', r10, range='m='//trim(toChar(m))//',l='//trim(toChar(l)), quiet=.true.)
      wz_siMtxPiv(nmo(2,m,l),:) = int(r10)
    enddo
  end do

  call AssertEqual( 'Pivot for semi-implicit matrix test 1', answer = wz_siMtxPiv, check = dyn_mtx00 % wz_siMtxPiv)

!!$  do l = 0, 10
!!$    do m = 0, 10
!!$!      if (m+l > 10) cycle
!!$      do k = 0, 1
!!$        write(*,*) 'wz_siMtxPiv(', 'nmo(1,', m, ',', l, '),', k+1, ')= ', &
!!$          & dyn_mtx00 % wz_siMtxPiv(nmo(1,m,l),k)
!!$      end do
!!$    end do
!!$  end do
!!$
!!$  do l = 0, 10
!!$    do m = 0, 10
!!$!      if (m+l > 10) cycle
!!$      do k = 0, 1
!!$        write(*,*) 'wz_siMtxPiv(', 'nmo(2,', m, ',', l, '),', k+1, ')= ', &
!!$          & dyn_mtx00 % wz_siMtxPiv(nmo(2,m,l),k)
!!$      end do
!!$    end do
!!$  end do

  !---------------------------------------------------------
  !  Solve サブルーチンのテスト
  !  Test 'Solve' subroutine
  !---------------------------------------------------------
  call MessageNotify('M', subname, 'Input ' // solve_src_file // '  ...' )

  allocate( wz_DivAvrTime(11**2, 0:1) )
  allocate( d10(2) )

  wz_DivAvrTime = 0.0_DP
  do l = 0, 10
    do m = 0, 10
      if (nmo(1,m,l) < 1) cycle
      call HistoryGet(solve_src_file, 'DivAvrTime1', d10, range='m='//trim(toChar(m))//',l='//trim(toChar(l)), quiet=.true.)
      wz_DivAvrTime(nmo(1,m,l),:) = d10
    end do
  end do

  do l = 0, 10
    do m = 0, 10
      if (nmo(2,m,l) < 1) cycle
      call HistoryGet(solve_src_file, 'DivAvrTime2', d10, range='m='//trim(toChar(m))//',l='//trim(toChar(l)), quiet=.true.)
      wz_DivAvrTime(nmo(2,m,l),:) = d10
    enddo
  end do

  call Solve( dyn_mtx00, wz_DivAvrTime ) ! (inout)

  call MessageNotify('M', subname, 'Input ' // solve_ans_file // '  ...' )

  allocate( wz_DivAvrTimeAns(11**2, 0:1) )
  allocate( wz_DivAvrTimeAnsL(11**2, 0:1) )
  allocate( wz_DivAvrTimeAnsG(11**2, 0:1) )

  wz_DivAvrTimeAns = 0.0_DP
  do l = 0, 10
    do m = 0, 10
      if (nmo(1,m,l) < 1) cycle
      call HistoryGet(solve_ans_file, 'DivAvrTimeAns1', d10, range='m='//trim(toChar(m))//',l='//trim(toChar(l)), quiet=.true.)
      wz_DivAvrTimeAns(nmo(1,m,l),:) = d10
    end do
  end do

  do l = 0, 10
    do m = 0, 10
      if (nmo(2,m,l) < 1) cycle
      call HistoryGet(solve_ans_file, 'DivAvrTimeAns2', d10, range='m='//trim(toChar(m))//',l='//trim(toChar(l)), quiet=.true.)
      wz_DivAvrTimeAns(nmo(2,m,l),:) = d10
    enddo
  end do

  where ( wz_DivAvrTimeAns < 0 )
    wz_DivAvrTimeAnsL = (wz_DivAvrTimeAns - 1.0e-13_DP) * 1.0001_DP
  elsewhere
    wz_DivAvrTimeAnsL = (wz_DivAvrTimeAns + 1.0e-13_DP) * 1.0001_DP
  end where

  call AssertLessThan( 'Solve test 1', answer = wz_DivAvrTimeAnsL, check = wz_DivAvrTime)

  where ( wz_DivAvrTimeAns < 0 )
    wz_DivAvrTimeAnsG = (wz_DivAvrTimeAns + 1.0e-13_DP) * 0.9999_DP
  elsewhere
    wz_DivAvrTimeAnsG = (wz_DivAvrTimeAns - 1.0e-13_DP) * 0.9999_DP
  end where

  call AssertGreaterThan( 'Solve test 2', answer = wz_DivAvrTimeAnsG, check = wz_DivAvrTime)

end program dyn_matrix_test