!= dcpam 主プログラム
!
!= dcpam main program
!
! Authors::   Yasuhiro Morikawa, Satoshi Noda, Yoshiyuki O. Takahashi, Shin-ichi Takehiro
! Version::   $Id: dcpam_main.f90,v 1.67 2025/09/17 22:00:00 takepiro Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2025. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

program dcpam_main_sltt_test_dp
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! モデルの使い方については {チュートリアル}[link:../../../doc/tutorial/rakuraku/] を
  ! 参照してください.
  !
  ! See {Tutorial}[link:../../../doc/tutorial/rakuraku/index.htm.en] for usage of the 
  ! model. 
  !

  use dc_message, only: MessageNotify

  ! メッセージ制御
  ! Message control
  !
  use mpi_messagecntl, only : DoesOutputMPIMessage, MPIMessageCntlInit
  
  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision. 
    &                 STDOUT, &  ! 標準出力の装置番号. Unit number of standard output
    &                 STRING, &  ! 文字列.       Strings. 
    &                 TOKEN      ! キーワード.   Keywords. 

  ! ファイル入出力補助
  ! File I/O support
  !
  use dc_iounit, only: FileOpen
  
  use dc_test, only : AssertEqual
  
  ! モジュール引用 ; USE statements
  !
  ! MPI
  !
  use mpi_wrapper, only : MPIWrapperInit, MPIWrapperFinalize, myrank

  ! コマンドライン引数処理
  ! Command line option parser
  !
  use option_parser, only: OptParseInit

  ! NAMELIST ファイル入力に関するユーティリティ
  ! Utilities for NAMELIST file input
  !
  use namelist_util, only: NmlutilInit, NmlutilMsg

  ! 物理定数設定
  ! Physical constants settings
  !
  use constants0, only: PI
  use constants, only: ConstantsInit, Rplanet

  ! 格子点設定
  ! Grid points settings
  !
#ifdef LIB_MPI  
  use gridset, only: GridSetInit, &
                     imax, &  ! 経度格子点数. 
                              ! Number of grid points in longitude
    &                jmax, &  ! 緯度格子点数. 
                              ! Number of grid points in latitude
    &                kmax, &  ! 鉛直層数. 
                              ! Number of vertical level
    &                imax_global ! 経度格子点数 (全球).
                                 ! Number of grid points in longitude on whole globe
  use ua_mpi_module_base, only : kc
#else
  use gridset, only: GridSetInit, &
                     imax, &  ! 経度格子点数. 
                              ! Number of grid points in longitude
    &                jmax, &  ! 緯度格子点数. 
                              ! Number of grid points in latitude
    &                kmax, &  ! 鉛直層数. 
                              ! Number of vertical level
    &           kc => kmax,& ! 並列鉛直層数.
                             ! Number of vertical level for parallel compuation
    &                imax_global ! 経度格子点数 (全球).
                                 ! Number of grid points in longitude on whole globe
#endif
  
  use axesset    , only : AxesSetInit, x_Lon, y_Lat

  use composition, only:  CompositionInit, ncmax

  use sltt_const , only : SLTTConstInit, &
       & jexmin, jexmax, jmaxh, jexglobalmin, jexglobalmax
  use sltt_extarr, only : SLTTExtArrInit, SLTTExtArrf, SlttExtArr
  use sltt_dp    , only : SLTTDPHor


  ! 宣言文 ; Declaration statements
  !
  implicit none

  character(*), parameter:: prog_name = 'dcpam_main'
                            ! 主プログラム名. 
                            ! Main program name

  character(STRING)      :: namelist_filename
                            ! NAMELIST ファイルの名称. 
                            ! NAMELIST file name

  ! 予報変数 (ステップ $ t-\Delta t $ , $ t $ , $ t+\Delta t $ )
  ! Prediction variables  (Step $ t-\Delta t $ , $ t $ , $ t+\Delta t $ )
  !
  real(DP), allocatable:: xyz_U (:,:,:)
                              ! $ u (t-\Delta t) $ .   東西風速. Eastward wind (m s-1)
  real(DP), allocatable:: xyz_V (:,:,:)
                              ! $ v (t-\Delta t) $ .   南北風速. Northward wind (m s-1)
  real(DP), allocatable:: xyzf_QMix(:,:,:,:)
                              ! $ q (t-\Delta t) $ .   混合比. Mass mixing ratio of constituents (1)
  real(DP)    , save, allocatable :: xyzf_QMix_dlon(:,:,:,:)
  real(DP)    , save, allocatable :: xyzf_QMix_dlat(:,:,:,:)
  real(DP)    , save, allocatable :: xyzf_QMix_dlonlat(:,:,:,:)


  real(DP)    , save, allocatable :: x_SinLon(:)
                              ! $\sin\lambda_S$
  real(DP)    , save, allocatable :: x_CosLon(:)
                              ! $\cos\lambda_S$
  real(DP)    , save, allocatable :: y_SinLat(:)
                              ! $\sin\varphai_S$
  real(DP)    , save, allocatable :: y_CosLat(:)
                              ! $\cos\varphai_S$
  real(DP)    , save, allocatable :: x_ExtLon(:)
                              ! $ x_LonSの拡張配列。
                              !Extended array of x_LonS.
  real(DP)    , save, allocatable :: y_ExtLat(:)
                              ! $ x_LatSの拡張配列。
                              !Extended array of x_LatS.

  real(DP)    , save, allocatable :: xyz_ExtU(:,:,:)
                              ! 東西風速の拡張配列（南半球）
                              ! Extended array (SH) of Zonal Wind        
  real(DP)    , save, allocatable :: xyz_ExtV(:,:,:)
                              ! 南北風速の拡張配列（南半球）
                              ! Extended array (SH) of Meridional Wind

  real(DP)    , save, allocatable :: xyzf_ExtQMix(:,:,:,:)
  real(DP)    , save, allocatable :: xyzf_ExtQMix_dlon(:,:,:,:)
  real(DP)    , save, allocatable :: xyzf_ExtQMix_dlat(:,:,:,:)
  real(DP)    , save, allocatable :: xyzf_ExtQMix_dlonlat(:,:,:,:)

  real(DP)    , save, allocatable :: xyz_DPLon(:,:,:)
  real(DP)    , save, allocatable :: xyz_DPLat(:,:,:)

  real(DP)    , save, allocatable :: xyz_DPLonAns(:,:,:)
  real(DP)    , save, allocatable :: xyz_DPLatAns(:,:,:)


  real(dp) :: lam0  = 180.0D0  ! 剛体回転流の回転ベクトル (経度, 度)
  real(dp) :: phi0  =  90.0D0  ! 剛体回転流の回転ベクトル (緯度, 度)
  real(dp) :: delta =   5.0D0  ! 剛体回転流の回転角

  ! 速度場パラメター
  real(DP) :: U0 = 1.0D1
!  real(DP) :: alpha = 0.0D0
  real(DP) :: alpha = 90.0D0

  real(DP) :: Deltime = 1.0D1 * 60.0D0

  ! NAMELIST 変数群
  ! NAMELIST group name
  !
  namelist /sltt_horadv_test_nml/  U0, lam0, phi0, delta, DelTime

  ! 判定誤差設定
!!$  integer, parameter :: check_digits = 3
!!$  integer, parameter :: ignore = -4
  integer, parameter :: check_digits = 5
  integer, parameter :: ignore = -6
!!$  integer, parameter :: check_digits = 10
!!$  integer, parameter :: ignore = -11
  
  real(DP) :: PM
  integer  :: i, j, k, n

  integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                            ! Unit number for NAMELIST file open
  integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                            ! IOSTAT of NAMELIST read
  
  ! 実行文 ; Executable statement
  !
  ! Initialize MPI
  !
  call MPIWrapperInit

  ! コマンドライン引数処理
  ! Command line option parser
  !
  call OptParseInit(       &
       & namelist_filename,    & ! (out)
       & prog_name            & ! (in )
       & )

  ! NAMELIST ファイル名入力
  ! Input NAMELIST file name
  !
  call NmlutilInit( &
       & namelist_filename  & ! (in)
       & )


  ! 格子点設定
  ! Grid points settings
  !
  call GridsetInit
  
  ! 組成に関わる配列の設定
  ! Settings of array for atmospheric composition
  !
  call CompositionInit

  ! 物理定数設定
  ! Physical constants settings
  !
  call ConstantsInit

  ! 座標データ設定
  ! Axes data settings
  !
  call AxessetInit
  
  call SLTTConstInit
  
  ! NAMELIST の読み込み
  ! NAMELIST is input
  !
  if ( trim(namelist_filename) /= '' ) then
     call FileOpen( unit_nml, &          ! (out)
          & namelist_filename, mode = 'r' ) ! (in)

     rewind( unit_nml )
     read( unit_nml, &                ! (in)
          & nml = sltt_horadv_test_nml, &  ! (out)
          & iostat = iostat_nml )        ! (out)
     close( unit_nml )

     call NmlutilMsg( iostat_nml, 'main' ) ! (in)
     if ( iostat_nml == 0 .AND. DoesOutputMPIMessage() ) write( STDOUT, nml = sltt_horadv_test_nml )
  end if

  !
  ! 一般軸
  !
  lam0  = lam0 * pi / 180.0_dp
  phi0  = phi0 * pi / 180.0_dp
  delta = delta * pi / 180.0_dp

  alpha = pi/2 - phi0

  allocate( xyz_U (0:imax-1, 1:jmax, 1:kmax) )
  allocate( xyz_V (0:imax-1, 1:jmax, 1:kmax) )
  allocate( xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
  allocate( xyzf_QMix_dlon(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
  allocate( xyzf_QMix_dlat(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
  allocate( xyzf_QMix_dlonlat(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )

  do j = 1, jmax
     do i = 0, imax-1
        xyz_U(i,j,:) &
             = U0 * (   cos(y_Lat(j))*cos(alpha) &
                      + sin(y_Lat(j))*cos(x_Lon(i)-lam0)*sin(alpha) )
        xyz_V(i,j,:) = U0 * sin(x_Lon(i)-lam0)*sin(alpha)
     end do
  end do
  xyzf_QMix = 0.0D0
  xyzf_QMix_dlon = 0.0D0
  xyzf_QMix_dlat = 0.0D0
  xyzf_QMix_dlonlat = 0.0D0


  allocate( x_SinLon(0:imax-1) )
  allocate( x_CosLon(0:imax-1) )
  allocate( y_SinLat(1:jmax) )
  allocate( y_CosLat(1:jmax) )
  do i = 0, imax-1
     x_SinLon(i) = sin( x_Lon(i) )
     x_CosLon(i) = cos( x_Lon(i) )
  end do
  do j = 1, jmax
     y_SinLat(j) = sin( y_Lat(j) )
     y_CosLat(j) = cos( y_Lat(j) )
  end do

  allocate( x_ExtLon( 0:imax_global-1 ) )
  allocate( y_ExtLat( jexmin:jexmax ) )

  allocate( xyz_ExtU(0:imax_global-1, jexmin:jexmax, 1:kc) )
  allocate( xyz_ExtV(0:imax_global-1, jexmin:jexmax, 1:kc) )

  allocate(xyzf_ExtQMix(0:imax_global-1, jexmin:jexmax, 1:kc, 1:ncmax))

  allocate(xyzf_ExtQMix_dlon(0:imax_global-1, jexmin:jexmax, 1:kc, 1:ncmax))
  allocate(xyzf_ExtQMix_dlat(0:imax_global-1, jexmin:jexmax, 1:kc, 1:ncmax))
  allocate(xyzf_ExtQMix_dlonlat(0:imax_global-1, jexmin:jexmax, 1:kc, 1:ncmax))

  allocate(xyz_DPLon(0:imax_global-1, 1:jmax, 1:kc))
                              ! 上流点経度（南半球）
                              ! Lon of the departure point (SH)
  allocate(xyz_DPLat(0:imax_global-1, 1:jmax, 1:kc))
                              ! 上流点緯度（南半球）
                              ! Lat of the departure point (SH)    

  allocate(xyz_DPLonAns(0:imax_global-1, 1:jmax, 1:kc))
                              ! 上流点経度（南半球）
                              ! Lon of the departure point (SH)

  allocate(xyz_DPLatAns(0:imax_global-1, 1:jmax, 1:kc))
                              ! 上流点緯度（南半球）
                              ! Lat of the departure point (SH)    


  ! 配列の分割と拡張
  ! Division and extension of arrays
  !
  call SLTTExtArrInit(            &
       & x_Lon, y_Lat,            & ! (in )
       & x_ExtLon, y_ExtLat       & ! (out)
       & )
  
  pm = -1.0_DP ! 配列拡張する際、極ごえ後に符号が変わる場合は -1.0を与える。そうでない場合は1.0を与える。
  ! -1.0 if the sign of value changes over the poles; if not 1.0. 
  !
  call SLTTExtArrf(                   &
       & xyzf_QMix_dlon, PM,          & ! (in)
       & xyzf_ExtQMix_dlon            & ! (out)
       & )


  pm = -1.0_DP ! 配列拡張する際、極ごえ後に符号が変わる場合は -1.0を与える。そうでない場合は1.0を与える。
  ! -1.0 if the sign of value changes over the poles; if not 1.0.
  !
  call SLTTExtArrf(                     &
       & xyzf_QMix_dlat, PM,            & ! (in)
       & xyzf_ExtQMix_dlat              & ! (out)
       & )

  pm = +1.0_DP ! 配列拡張する際、極ごえ後に符号が変わる場合は -1.0を与える。そうでない場合は1.0を与える。
  ! -1.0 if the sign of value changes over the poles; if not 1.0.
  !
  call SLTTExtArrf(                     &
       & xyzf_QMix_dlonlat, PM,         & ! (in)
       & xyzf_ExtQMix_dlonlat           & ! (out)
       & )

  PM = 1.0_DP
  do n = 1, ncmax
     call SLTTExtArrf(                         &
          & xyzf_QMix(:,:,:,n), 1.0D0,         & ! (in)
          & xyzf_ExtQMix(:,:,:,n)              & ! (out)
          & )
  end do

#ifdef AXISYMMETRY
  xyz_ExtU = 0.0D0
#else
  PM = -1.0_DP
  call SLTTExtArr(                 &
       & xyz_U, PM,                & ! (in)
       & xyz_ExtU                  & ! (out)
       & )
#endif
  PM = -1.0_DP
  call SLTTExtArr(                 &
       & xyz_V, PM,                & ! (in)
       & xyz_ExtV                  & ! (out)
       & )

  ! 上流点の計算
  ! estimation of departure point
  call SLTTDPHor(                                    &
       & DelTime, x_Lon, y_Lat, y_SinLat, y_CosLat,  & ! (in)
       & jexmin, jexmax,                             & ! (in)
       & x_ExtLon, y_ExtLat, xyz_ExtU, xyz_ExtV,     & ! (in)
       & xyz_DPLon, xyz_DPLat                        & ! (out)
       & )
  
  ! 上流点の計算 (解析解)
  ! estimation of departure point
  do k=1,kc
     do j=1,jmax
        do i=0,imax_global-1
           xyz_DPLonAns(i,j,k) &
                = x_ExtLon(i) - xyz_ExtU(i,j,k)/(RPlanet*cos(y_Lat(j))) * 2.0D0*DelTime
           xyz_DPLatAns(i,j,k) = y_Lat(j) - xyz_ExtV(i,j,k)/RPlanet * 2.0D0*DelTime
        end do
     end do
  end do
  xyz_DPLonAns = modulo(xyz_DPLonAns,2*PI)

  !----------------------------------------------------------------------
  ! 結果チェック
  call AssertEqual(&
       message='SLTTDPHor (DPLon)',                                &
       answer = xyz_DPLonAns,                                      &
       check = xyz_DPLon,                                          &
       significant_digits = check_digits, ignore_digits = ignore   &
       )
  call AssertEqual(&
       message='SLTTDPHor (DPLat)',                                &
       answer = xyz_DPLatAns,                                      &
       check = xyz_DPLat,                                          &
       significant_digits = check_digits, ignore_digits = ignore   &
       )

!!$  write(6,*) xyz_DPLatAns
!!$  write(6,*)
!!$  write(6,*) xyz_DPLat
  
  ! Finalize MPI
  !
  call MPIWrapperFinalize

  
end program dcpam_main_sltt_test_dp