!= dcpam 主プログラム
!
!= dcpam main program
!
! Authors::   Yasuhiro Morikawa, Satoshi Noda, Yoshiyuki O. Takahashi, Shin-ichi Takehiro
! Version::   $Id: dcpam_main.f90,v 1.67 2026/05/21 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
  !
  ! <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, r_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, SLTTHorAdv, 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:: 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:: xyr_SigDot (:,:,:)
                              ! $ u (t-\Delta t) $ .   鉛直風速. vertial wind (m s-1)
  real(DP), allocatable:: xyzf_QMixA(:,:,:,:)
                              ! $ q (t+\Delta t) $ .   混合比. Mass mixing ratio of constituents (1)
  real(DP), allocatable:: xyzf_QMixLinA(:,:,:,:)
                              ! $ q (t+\Delta t) $ .   混合比. Mass mixing ratio
  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

  real(DP) :: U0 = 1.0D0      ! 東西流の速度振幅
  character(2) :: vtype='11'  ! 速度分布タイプ (m,n) : (11,21)
  
  integer :: i, j, k, n
  
  integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                            ! Unit number for NAMELIST file open
  integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                            ! IOSTAT of NAMELIST read

  ! NAMELIST 変数群
  ! NAMELIST group name
  !
  namelist /sltt_zonaladv_test_nml/  U0, Vtype

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

  ! 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_zonaladv_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_zonaladv_test_nml )
  end if

  !----------------------------------------------------------------------
  ! 変数割付
  !
  call HistoryAutoAddVariable( 'U' ,            &  ! (in)
        & (/ 'lon ', 'lat ', 'sig ','time'/),   &  ! (in)
        & 'U', 'm/s' )                             ! (in)
  call HistoryAutoAddVariable( 'V' ,            &  ! (in)
        & (/ 'lon ', 'lat ', 'sig ','time'/),   &  ! (in)
        & 'V', 'm/s' )                             ! (in)
  call HistoryAutoAddVariable( 'SigDot' ,       &  ! (in)
        & (/ 'lon ', 'lat ', 'sigm','time'/),   &  ! (in)
        & 'SigDot', '1/s' )                        ! (in)
  call HistoryAutoAddVariable( 'QVap',          &  ! (in)
        & (/ 'lon ', 'lat ', 'sig ','time'/),   &  ! (in)
        & 'QVap', 'kg kg-1' )                      ! (in)
  call HistoryAutoAddVariable( 'QLiq',          &  ! (in)
        & (/ 'lon ', 'lat ', 'sig ','time'/),   &  ! (in)
        & 'QLiq', 'kg kg-1' )                      ! (in)
  call HistoryAutoAddVariable( 'QSol',          &  ! (in)
        & (/ 'lon ', 'lat ', 'sig ','time'/),   &  ! (in)
        & 'QSol', 'kg kg-1' )                      ! (in)

  allocate( xyz_U (0:imax-1, 1:jmax, 1:kmax) )
  allocate( xyz_V (0:imax-1, 1:jmax, 1:kmax) )
  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_QMixLinA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
  allocate( xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )

  !----------------------------------------------------------------------
  ! 速度設定
  !
  if ( Vtype == '11') then
     do k = 1, kmax
        do j = 1, jmax
           do i = 0, imax-1
              xyz_U(i,j,k) = -U0*PI*cos(PI*z_Sigma(k))*cos(x_Lon(i))*cos(y_Lat(j))
              xyz_V(i,j,k) = 0.0D0
           end do
        end do
     end do
     do k = 0, kmax
        do i = 0, imax-1
           xyr_SigDot(i,:,k) = -U0/RPlanet * sin(PI*r_Sigma(k))*sin(x_Lon(i))
        end do
     end do
  else if ( Vtype == '21') then
     do k = 1, kmax
        do j = 1, jmax
           do i = 0, imax-1
              xyz_U(:,j,k) = -U0*PI*cos(PI*z_Sigma(k))*cos(2*x_Lon(i))*cos(y_Lat(j))
              xyz_V(:,j,k) = 0.0D0
           enddo
        enddo
     enddo
     do k = 0, kmax
        do j = 0, imax-1
           xyr_SigDot(i,:,k) = -2*U0/RPlanet * sin(PI*r_Sigma(k))*sin(x_Lon(i))
        enddo
     enddo
  else
     call MessageNotify('E','main','n,m values not supported')
  end if
  
  !----------------------------------------------------------------------
  ! トレーサー分布設定
  !
  do k = 1, kmax
     do j = 1, jmax
        do i = 0, imax-1
           xyzf_QMixB(i,j,k,:) = f_distribution(x_Lon(i), y_Lat(j),z_Sigma(k))
           xyzf_QMixN(i,j,k,:) = f_distribution(x_Lon(i), y_Lat(j),z_Sigma(k))
        end do
     end do
  end do

  call HistoryAutoPut( TimeN, 'U', xyz_U )
  call HistoryAutoPut( TimeN, 'V', xyz_V )
  call HistoryAutoPut( TimeN, 'SigDot', xyr_SigDot )
  call HistoryAutoPut( TimeN, 'QVap',  xyzf_QMixN(:,:,:,1) )
  call HistoryAutoPut( TimeN, 'QLiq', xyzf_QMixN(:,:,:,2) )
  call HistoryAutoPut( TimeN, 'QSol', xyzf_QMixN(:,:,:,3) )
  
  !----------------------------------------------------------------------
  ! 時間積分
  ! Time integration
  !
  loop_time : do while ( TimeB < EndTime )

     xyzf_QMixLinA = xyzf_QMixB
     
     xyzf_QMixB = SLTTHorAdv( xyzf_QMixB, xyz_U, xyz_V,    &
       &                      xyzf_QMixLinA = xyzf_QMixLinA )

     xyzf_QMixA = SLTTVerAdv( xyr_SigDot, xyzf_QMixB )
     
     call HistoryAutoPut( TimeA, 'U', xyz_U )
     call HistoryAutoPut( TimeA, 'V', xyz_V )
     call HistoryAutoPut( TimeN, 'QVap',  xyzf_QMixN(:,:,:,1) )
     call HistoryAutoPut( TimeN, 'QLiq', xyzf_QMixN(:,:,:,2) )
     call HistoryAutoPut( TimeN, 'QSol', xyzf_QMixN(:,:,:,3) )

     ! 次の時間ステップに向けて予報変数の入れ替え
     ! 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
  
  ! Finalize MPI
  !
  call MPIWrapperFinalize

contains
  
  !----------------------------------------------------------------------
  ! トレーサー分布
  ! 
  real(dp) function f_distribution(lam, phi, sig)
    real(dp), intent(in) :: lam, phi, sig 
    character(2), save :: disttype='CB'  ! 分布タイプ
                                         ! CB : cosine Bell
                                         ! 11,21 : (l,n)
    real(dp), save  :: radius            ! cosine Bell 型の半径

    namelist /fdist_nml/ disttype, radius

    logical, save :: first = .true.
    
    
    if ( first ) then
       first = .false.
       radius = 1.0D0/3.0D0         ! デフォルト値
       
       ! 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 = fdist_nml, &  ! (out)
               & iostat = iostat_nml )        ! (out)
          close( unit_nml )

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

    if ( disttype == 'CB' ) then
       f_distribution = CosBell(lam, phi, sig, radius)
    else
       f_distribution = Distln(lam, phi, sig, disttype)
    end if

  end function f_distribution

  ! --- 元の分布関数 (cos Bell) ---
  real(dp) function CosBell(lam, phi, sig, radius)
    real(dp), intent(in) :: lam, phi, sig, radius
    real(dp) :: dist, center_lam, center_phi, center_sig

    center_lam = PI
    center_sig = 0.5_dp
    center_phi = 0.0_dp

    ! 角度距離の計算 (単純化のためユークリッド距離、厳密には大円距離も可)
    dist = sqrt((sig - center_sig)**2 + (phi - center_phi)**2 + (lam - center_lam)**2)

    if (dist < radius) then
       CosBell = cos(0.5_dp * pi * dist / radius)**2
    else
       CosBell = 0.0_dp
    end if
    
  end function CosBell

  ! --- 元の分布関数 (Y_l^m) ---
  real(dp) function Distln(lam, phi, sig, ln)
    real(dp), intent(in)     :: lam, phi, sig
    character(2), intent(in) :: ln

    if ( ln == '11') then
       distln = sin(PI*sig) * cos(lam) * cos(phi)
    else if ( ln == '21') then
       distln = sin(PI*sig) * cos(2*lam) * cos(phi)
    else
       call MessageNotify('E','Distln','n,m values not supported')
    end if
    
  end function Distln

end program dcpam_main_sltt_test