!= 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_veradv
  !
  ! <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

  ! 時刻管理
  ! Time control
  !
  use timeset, only: TimesetInit, TimesetProgress, &
    & TimeB, &                ! ステップ $ t - \Delta t $ の時刻.
                              ! Time of step $ t - \Delta t $.
    & TimeN, &                ! ステップ $ t $ の時刻.
                              ! Time of step $ t $.
    & TimeA, &                ! ステップ $ t + \Delta t $ の時刻.
                              ! Time of step $ t + \Delta t $.
    & EndTime, &              ! 計算終了時刻.
                              ! End time of calculation
    & DelTime                 ! $ \Delta t $ [s]

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

  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: GridSetInit, &
                     imax, &  ! 経度格子点数. 
                              ! Number of grid points in longitude
    &                jmax, &  ! 緯度格子点数. 
                              ! Number of grid points in latitude
    &                kmax     ! 鉛直層数. 
                              ! Number of vertical level

  use axesset    , only : AxesSetInit, x_Lon, y_Lat, z_Sigma

  use composition, only:  CompositionInit, ncmax
  !
  !

  ! 出力ファイルの基本情報管理
  ! Management basic information for output files
  !
  use fileset, only: FilesetInit
  
  ! ヒストリデータファイルの初期化
  ! Initialization of history data files
  !
  use history_file_io, only: HistoryFileOpen, HistoryFileClose

  
  use sltt       , only: SLTTinit, SLTTVerAdv, SLTTCheckInit

  use gtool_historyauto, only: &
       & HistoryAutoAddVariable, HistoryAutoPut

  use gtool_history, only: HistoryClose

  ! 宣言文 ; 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:: xyr_SigDot (:,:,:)
                              ! $ u (t-\Delta t) $ .   東西風速. Eastward wind (m s-1)
  real(DP), allocatable:: xyzf_QMixA(:,:,:,:)
                              ! $ q (t-\Delta t) $ .   混合比. Mass mixing ratio of constituents (1)
   real(DP), allocatable:: xyzf_QMixN(:,:,:,:)
                              ! $ q (t) $ .            混合比. Mass mixing ratio of constituents (1)
   real(DP), allocatable:: xyzf_QMixB(:,:,:,:)
                              ! $ q (t-\Delta t) $ .   混合比. Mass mixing ratio of constituents (1)
  real(DP), allocatable:: xyzf_QMixS(:,:,:,:)
                              ! $ q (t-\Delta t) $ .   混合比. Mass mixing ratio of constituents (1)

  real(DP) :: sig               ! t=0 での Sigma 座標
  real(DP) :: Time0             ! 初期時刻
  
  integer :: i, j, k, n
  
  integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                            ! Unit number for NAMELIST file open
  integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                            ! IOSTAT of NAMELIST read

  !
  ! 物質分布パラメター (Gauss 分布)
  !
  real(dp) :: lon0 = 180.0D0   ! 中心経度
  real(dp) :: dlon = 20.0D0    ! 経度幅

  real(dp) :: lat0 = 0.0D0     ! 中心緯度
  real(dp) :: dlat = 10.0D0    ! 緯度幅

  real(dp) :: sig0 = 0.5D0     ! 中心 Sigma
  real(dp) :: dsig = 0.1D0     ! Sigma 幅

  real(DP) :: SigDot0  = 1.0D0 ! 平行流の振幅

  real(dp) :: lon0r, lat0r
  real(dp) :: dlonr, dlatr

  ! NAMELIST 変数群
  ! NAMELIST group name
  !
  namelist /sltt_veradv_test_nml/  &
       SigDot0, lon0, dlon, lat0, dlat, sig0, dsig

  ! 判定誤差設定
  integer, parameter :: check_digits = 3
  integer, parameter :: ignore = -3

  ! SLTT ルーチンチェック出力スイッチ
  logical :: sltt_check = .false.
  ! logical :: sltt_check = .true.
  
  !======================================================================
  ! 実行文 ; 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)
       & )

  ! 時刻管理
  ! Time control
  !
  call TimesetInit

  ! 出力ファイルの基本情報管理
  ! Management basic information for output files
  !
  call FilesetInit

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

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

  ! 座標データ設定
  ! Axes data settings
  !
  call AxessetInit

  ! ヒストリデータファイルの初期化
  ! Initialization of history data files
  !
  call HistoryFileOpen
  
  call SLTTInit
  if ( sltt_check ) call SLTTCheckInit
 
  ! 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_veradv_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_veradv_test_nml )
  end if

  lon0r = lon0 * PI/180.0D0
  dlonr = dlon * PI/180.0D0

  lat0r = lat0 * PI/180.0D0
  dlatr = dlat * PI/180.0D0

  call HistoryAutoAddVariable( 'SigDot' ,       &  ! (in)
        & (/ 'lon ', 'lat ', 'sigm','time'/),   &  ! (in)
        & 'SigDot', '1/s' )                        ! (in)
  call HistoryAutoAddVariable( 'QMix',          &  ! (in)
        & (/ 'lon ', 'lat ', 'sig ','time'/),   &  ! (in)
        & 'QMix', 'kg kg-1' )                      ! (in)
  call HistoryAutoAddVariable( 'QMixS' ,        &  ! (in)
        & (/ 'lon ', 'lat ', 'sig ', 'time'/),  &  ! (in)
        & 'QMixS', 'kg kg-1' )                     ! (in)
  
  allocate( xyr_SigDot (0:imax-1, 1:jmax, 0:kmax) )
  allocate( xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
  allocate( xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
  allocate( xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
  allocate( xyzf_QMixS(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )

  xyr_SigDot = SigDot0
  
  ! 初期値
  !
  do k = 1, kmax
     do j = 1, jmax
        do i = 0, imax-1
           xyzf_QMixN(i,j,k,:) = fdist( x_Lon(i), y_Lat(j), z_Sigma(k) )
        end do
     end do
  end do

  call HistoryAutoPut( TimeN, 'SigDot', xyr_SigDot )
  call HistoryAutoPut( TimeN, 'QMix',  xyzf_QMixN(:,:,:,1) )
  call HistoryAutoPut( TimeN, 'QMixS', xyzf_QMixN(:,:,:,1) )

  !
  ! DetTime 後
  !
  Time0 = TimeN
  do k = 1, kmax
     do j = 1, jmax
        do i = 0, imax-1
           sig = z_Sigma(k) - DelTime * (xyr_SigDot(i,j,k)+xyr_SigDot(i,j,k-1))/2.0D0
           xyzf_QMixB(i,j,k,:) = fdist( x_Lon(i), y_Lat(j), sig )
        end do
     end do
  end do
  
  ! 時間積分
  ! Time integration
  !
  loop_time : do while ( TimeB < EndTime )

     xyzf_QMixA = SLTTVerAdv( xyr_SigDot, xyzf_QMixB )

     do k = 1, kmax
        do j = 1, jmax
           do i = 0, imax-1
              sig = z_Sigma(k) - (TimeA -Time0)* (xyr_SigDot(i,j,k)+xyr_SigDot(i,j,k-1))/2.0D0
              xyzf_QMixS(i,j,k,:) = fdist( x_Lon(i), y_Lat(j), sig )
           end do
        end do
     end do

     call HistoryAutoPut( TimeA, 'SigDot', xyr_SigDot )
     call HistoryAutoPut( TimeA, 'QMix',  xyzf_QMixB(:,:,:,1) )
     call HistoryAutoPut( TimeA, 'QMixS', xyzf_QMixB(:,:,:,1) )

     ! 次の時間ステップに向けて予報変数の入れ替え
     ! Exchange prediction variables for the next time step
     !
     xyzf_QMixB = xyzf_QMixN ; xyzf_QMixN = xyzf_QMixA

     ! 時刻の進行
     ! Progress time
     !
     call TimesetProgress

  end do loop_time

  call HistoryFileClose
  if ( sltt_check ) call HistoryClose
  
  !----------------------------------------------------------------------
  ! 結果チェック
  call AssertEqual(&
       message='QMix',                                             &
       answer = xyzf_QMixS(:,:,:,1),                               &
       check = xyzf_QMixA(:,:,:,1),                                &
       significant_digits = check_digits, ignore_digits = ignore   &
       )

  ! Finalize MPI
  !
  call MPIWrapperFinalize

contains
  !
  ! 初期分布
  !
  function fdist( lon, lat, sig )

    use constants0, only : PI
    
    real(dp), intent(IN) :: lon, lat, sig
    real(dp)             :: fdist

    fdist = exp( -(lon-lon0r)**2/dlonr**2 &
                 -(lat-lat0r)**2/dlatr**2 &
                 -(sig-sig0)**2/dsig**2 )

  end function fdist
  
end program dcpam_main_sltt_test_veradv