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

  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, 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:: 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 of constituents (1)
  real(DP), allocatable:: xyzf_QMixS(:,:,:,:)
                              ! $ q (t-\Delta t) $ .   混合比. Mass mixing ratio of constituents (1)

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

  real(DP) :: U0               ! 剛体回転流の速度振幅
  real(DP) :: alpha            ! 剛体回転流の角度パラメター
  real(DP) :: rotrad = 0.0D0   ! 回転角度
  
  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_horadv_test_nml/  lam0, phi0, delta

  ! 判定誤差設定
  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_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

  ! 1. 回転パラメータの設定 (数値の後ろに _dp を付与)
  !
  ! 極軸回転
!!$  lam0 = 180.0_dp * pi / 180.0_dp
!!$  phi0 = 90.0_dp * pi / 180.0_dp
!!$  delta = 10.0_dp * pi / 180.0_dp
  !
  ! 赤道軸回転
!!$  lam0 = 180.0_dp * pi / 180.0_dp
!!$  phi0 = 0.0_dp * pi / 180.0_dp
!!$  delta = 5.0_dp * pi / 180.0_dp
  !
  ! 一般軸
  !
  lam0  = lam0 * pi / 180.0_dp
  phi0  = phi0 * pi / 180.0_dp
  delta = delta * pi / 180.0_dp

  U0    = 4*RPlanet * delta/(EndTime)
  alpha = pi/2 - phi0

  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( 'QMixS' ,        &  ! (in)
        & (/ 'lon ', 'lat ', 'sig ', 'time'/),  &  ! (in)
        & 'QMixS', 'kg kg-1' )                     ! (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( 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) )
  allocate( xyzf_QMixS(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )

  rotrad = - U0 * DelTime/Rplanet
  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)
        !
        do n=1,ncmax
           do k=1,kmax
              xyzf_QMixB(:,:,k,n) = calc_analyticsolution( lam0, phi0, rotrad )
           end do
        end do
     end do
  end do

  rotrad = 0.0D0
  do j = 1, jmax
     do i = 0, imax-1
        xyzf_QMixN(i,j,:,:) = f_distribution(x_Lon(i),y_Lat(j))
     end do
  end do

  call HistoryAutoPut( TimeN, 'U', xyz_U )
  call HistoryAutoPut( TimeN, 'V', xyz_V )
  call HistoryAutoPut( TimeN, 'QMixS', xyzf_QMixN(:,:,:,1) )
  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_QMixA = SLTTHorAdv( xyzf_QMixB, xyz_U, xyz_V,    &
       &                      xyzf_QMixLinA = xyzf_QMixLinA )

     rotrad = rotrad + U0 * DelTime/Rplanet
     
     do n=1,ncmax
        do k=1,kmax
           xyzf_QMixS(:,:,k,n) = calc_analyticsolution( lam0, phi0, rotrad )
        end do
     end do

     call HistoryAutoPut( TimeA, 'U', xyz_U )
     call HistoryAutoPut( TimeA, 'V', xyz_V )
     call HistoryAutoPut( TimeA, 'QMixS', xyzf_QMixS(:,:,:,1) )
     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
  
  !----------------------------------------------------------------------
  ! 結果チェック
  call AssertEqual(&
       message='QMix',                                             &
       answer = xyzf_QMixS(:,:,:,1),                               &
       check = xyzf_QMixA(:,:,:,1),                                &
       significant_digits = check_digits, ignore_digits = ignore   &
       )

  ! Finalize MPI
  !
  call MPIWrapperFinalize

contains

  real(dp) function f_distribution(lam, phi)
    real(dp), intent(in) :: lam, phi
    character(2), save :: disttype='CB'  ! 分布タイプ
                                         ! CB : cosine Bell
                                         ! 10,11,20,21,22 : 球面調和函数
    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, radius)
    else
       f_distribution = Ylm(lam, phi, disttype)
    end if

  end function f_distribution

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

    center_lam = PI/2
    center_phi = 0.0_dp

    ! 角度距離の計算 (単純化のためユークリッド距離、厳密には大円距離も可)
    dist = sqrt((lam - center_lam)**2 + (phi - center_phi)**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 Ylm(lam, phi, lm)
    real(dp), intent(in)     :: lam, phi
    character(2), intent(in) :: lm

    if ( lm == '10') then
       Ylm = sin(phi)
    else if ( lm == '11') then
       Ylm = cos(phi)*sin(lam)
    else if ( lm == '20') then
       Ylm = 0.5D0*(3.0D0*sin(phi)**2-1.0D0)
    else if ( lm == '21') then
       Ylm = -3.0D0*sin(phi)*cos(phi)*sin(lam)
    else if ( lm == '22') then
       Ylm = 3.0D0*cos(phi)**2*sin(2.0D0*lam)
    else
       call MessageNotify('E','Ylm','l,m values not supported')
    end if

  end function Ylm

  function calc_analyticsolution( lam0, phi0, delta ) result(xy_fsol)
    
    real(dp), intent(IN) :: lam0, phi0, delta
    real(dp)             :: xy_fsol( 0:imax-1, 1:jmax )

    real(dp) :: axis(3), rot_mat(3,3)
    real(dp) :: xp, yp, zp, x_old, y_old, z_old
    real(dp) :: old_lam, old_phi
    integer :: i, j

   
    ! 移流計算 (解析解)
    !
    ! 2. 回転軸ベクトルと逆回転行列の生成
    axis(1) = cos(phi0) * cos(lam0)
    axis(2) = cos(phi0) * sin(lam0)
    axis(3) = sin(phi0)
    call get_rotation_matrix(axis, -delta, rot_mat)

    ! 3. 各グリッド点での計算
    do j = 1, jmax
       do i = 0, imax-1
          xp = cos(y_Lat(j)) * cos(x_Lon(i))
          yp = cos(y_Lat(j)) * sin(x_Lon(i))
          zp = sin(y_Lat(j))

          ! 逆回転 R^T * v' (回転行列の逆行列は転置)
          x_old = rot_mat(1,1)*xp + rot_mat(1,2)*yp + rot_mat(1,3)*zp
          y_old = rot_mat(2,1)*xp + rot_mat(2,2)*yp + rot_mat(2,3)*zp
          z_old = rot_mat(3,1)*xp + rot_mat(3,2)*yp + rot_mat(3,3)*zp

          old_phi = asin(max(-1.0_dp, min(1.0_dp, z_old)))
          old_lam = atan2(y_old, x_old)

          xy_fsol(i,j) = f_distribution(old_lam, old_phi)
       end do
    end do

  end function calc_analyticsolution
  
  ! --- 回転行列生成ルーチン (倍精度) ---
  subroutine get_rotation_matrix(n, d, R)
    real(dp), intent(in) :: n(3), d
    real(dp), intent(out) :: R(3,3)
    real(dp) :: c, s, t
    c = cos(d); s = sin(d); t = 1.0_dp - c
    R(1,1) = t*n(1)*n(1) + c
    R(1,2) = t*n(1)*n(2) - n(3)*s
    R(1,3) = t*n(1)*n(3) + n(2)*s
    R(2,1) = t*n(2)*n(1) + n(3)*s
    R(2,2) = t*n(2)*n(2) + c
    R(2,3) = t*n(2)*n(3) - n(1)*s
    R(3,1) = t*n(3)*n(1) - n(2)*s
    R(3,2) = t*n(3)*n(2) + n(1)*s
    R(3,3) = t*n(3)*n(3) + c
  end subroutine get_rotation_matrix
  
end program dcpam_main_sltt_test