!= dcpam 主プログラム
!
!= dcpam main program
!
! Authors::   Yasuhiro Morikawa, Satoshi Noda, Yoshiyuki O. Takahashi, Shin-ichi Takehiro
! Version::   $Id: dcpam_main_hs94.f90,v 1.67 2026/08/15 00: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_calc_tem
  !
  ! <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 statements
  !
  ! 時刻管理
  ! Time control
  !
  use timeset, only: 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]

!!$  ! リスタートデータ入出力
!!$  ! Restart data input/output
!!$  !
!!$  use restart_file_io, only: RestartFileOutPut

  ! ヒストリデータ出力
  ! History data output
  !
  use gtool_historyauto, only: HistoryAutoPut, HistoryAutoAllVarFix

  ! 組成に関わる配列の設定
  ! Settings of array for atmospheric composition
  !

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

  ! 座標データ設定
  ! Axes data settings
  !
  use axesset, only: &
       x_Lon, &
                            ! $ \lambda $ [rad.] . 経度. Longitude
       y_Lat, &
                            ! $ \varphi $ [rad.] . 緯度. Latitude
       z_Sigma, &
                            ! $ \sigma $ レベル (整数). 
                            ! Full $ \sigma $ level
       r_Sigma
                            ! $ \sigma $ レベル (半整数). 
                            ! Half $ \sigma $ level
  !
  !
  use dc_message, only: MessageNotify

  ! 物理・数学定数
  use constants0, only : PI
                            ! $ \pi $ .
                            ! 円周率.  Circular constant


  ! 物理定数設定
  ! Physical constants settings
  !
  use constants, only: ConstantsInit, &
    & RPlanet,&
                            ! $ a $ [m]. 
                            ! 惑星半径. 
                            ! Radius of planet
    & GasRDry, &
                            ! $ R $ [J kg-1 K-1]. 
                            ! 乾燥大気の気体定数. 
                            ! Gas constant of air
    & CpDry,&
                            ! $ C_p $ [J kg-1 K-1]. 
                            ! 乾燥大気の定圧比熱. 
                            ! Specific heat of air at constant pressure
    & Omega,&
                            ! $ \Omega $ [s-1]. 
                            ! 回転角速度. 
                            ! Angular velocity
    & Grav
                            ! $ g $ [m s-2]. 
                            ! 重力加速度. 
                            ! Gravitational acceleration

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision. 
    &                 STRING, &  ! 文字列.       Strings. 
    &                 TOKEN      ! キーワード.   Keywords. 

  !
  ! Calc Mass stream function
  !
  use mass_streamfunc, only :  yr_MassStreamFunc

  !
  ! TEM analysis
  !
  use tem, only : TEMOutput, TEMOutputInit

  ! 宣言文 ; Declaration statements
  !
  implicit none

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

  ! 予報変数
  !
  real(DP), allocatable :: xyz_U (:,:,:)
                              ! $ u $ .   東西風速. Eastward wind (m s-1)
  real(DP), allocatable :: xyz_V (:,:,:)
                              ! $ v $ .   南北風速. Northward wind (m s-1)
  real(DP), allocatable :: xyr_SigDot (:,:,:)
                              ! $ \dot\sigma $ . 鉛直流
                              ! vertical flow
  real(DP), allocatable :: xyz_Temp (:,:,:)
                              ! $ T $ .   温度. Temperature (K)
  real(DP), allocatable :: xyz_QVap(:,:,:)
                              ! $ q $ .   混合比. Mass mixing ratio of constituents (1)
  real(DP), allocatable :: xy_Ps (:,:)
                              ! $ p_s $ . 地表面気圧. Surface pressure (Pa)

  real(DP), allocatable :: xyz_EpLat (:,:,:)
               ! $ Ep_lat $ . EP フラックス緯度成分. 
               ! latitudinal EP flux

  real(DP), allocatable :: xyr_EpSig (:,:,:)
               ! $ Ep_lat $ . EP フラックス鉛直成分. 
               ! vertical EP flux

  real(DP), allocatable :: xyz_PsVRes (:,:,:)
               ! $ vres $ . 残差南北風速. 
               ! Residual northward wind

  real(DP), allocatable :: xyr_PsSigDotRes (:,:,:)
               ! $ \dot\sigma_res $ . 残差鉛直風速. 
               ! Residual vertical wind

  real(DP), allocatable :: xyr_StrRes (:,:,:)
               ! $ \psi_res $ . 残差質量流線関数. 
               ! Residual mass stream fkuction

  real(DP), allocatable :: yr_Msf    (:,:)
               ! $ \psi $ . 質量流線関数. 
               ! Mass stream fkuction


  integer  :: ntime    ! 時間ステップ数
  real(DP) :: ctime    ! 時刻
    
  ! 実行文 ; Executable statement
  !

  ! 主プログラムの初期化 (内部サブルーチン)
  ! Initialization for the main program (Internal subroutine)
  !
  call MainInit

  ! 時間積分
  ! Time integration
  !
  loop_time : do while ( TimeB < EndTime )
     
     ! データ読み込み
     ! Reading data
     !
     call read_history( TimeN )
  
     ! 診断変数計算
     ! Calculation of diagnostic variables
     !
     yr_MSF = yr_MassStreamFunc( xy_Ps, xyz_V )
     call HistoryAutoPut( TimeN, 'Msf', spread(yr_MSF,1,imax) )

     !
     ! TEM 計算
     !
     call TEMOutput(                                               &
          & xyz_U,  xyz_V, xyr_SigDot, xyz_QVap, xyz_Temp, xy_Ps   & ! (in)
          & )

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

  ! 時間積分終了
  ! Time integration is finished
  !
  end do loop_time
 
  ! 主プログラムの終了処理 (内部サブルーチン)
  ! Termination for the main program (Internal subroutine)
  !
  call MainTerminate

contains

  !-------------------------------------------------------------------

  subroutine MainInit
    !
    ! 主プログラムの初期化手続き. 
    !
    ! Initialization procedure for the main program. 
    !

    ! MPI
    !
    use mpi_wrapper, only : MPIWrapperInit

    !
    ! MPI 計算での Message 制御
    ! Message control for MPI calculation
    !
    use mpi_messagecntl, only : MPIMessageCntlInit

    use dc_message, only: MessageNotify

    ! コマンドライン引数処理
    ! 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, TimesetDelTimeHalf, &
      & TimeN                 ! ステップ $ t $ の時刻. Time of step $ t $. 

    ! 出力ファイルの基本情報管理
    ! Management basic information for output files
    ! 
    use fileset, only: FilesetInit

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: GridsetInit, &
      &                  nmax, & ! 最大全波数. 
                                 ! Maximum truncated wavenumber
      &                  imax, & ! 経度格子点数. 
                                 ! Number of grid points in longitude
      &                  jmax, & ! 緯度格子点数. 
                                 ! Number of grid points in latitude
      &                  kmax, & ! 鉛直層数. 
                                 ! Number of vertical level
      &                kslmax, & ! 地下の鉛直層数. 
                                 ! Number of subsurface vertical level
      &                ksimax    ! 海氷の鉛直層数. 
                                 ! Number of sea ice vertical level

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only: CompositionInit


    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: AxessetInit

    ! ヒストリデータ出力
    ! History data output
    !
    use history_file_io, only: HistoryFileOpen
    use gtool_history, only: HistoryGetPointer
    use gtool_historyauto, only: &
      & HistoryAutoAddVariable, HistoryAutoAddAttr, HistoryAutoPut

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: &
      & STDOUT               ! 標準出力の装置番号. Unit number of standard output

    ! 時系列データの読み込み
    ! Reading time series
    !
    use read_time_series, only : ReadTimeSeriesInit

    !
    ! Calc Mass stream function
    !
    use mass_streamfunc, only :  yr_MassStreamFunc

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level,
    ! and calculate pressure and height
    !
    use auxiliary, only: AuxVarsInit

    ! メッセージ制御
    ! Message control
    !
    use mpi_messagecntl, only : DoesOutputMPIMessage

    ! 宣言文 ; Declaration statements
    !
    implicit none

    character(*), parameter:: version = &
      & '$Name:  $' // &
      & '$Id: dcpam_main_hs94.f90,v 1.67 2026/08/15 00:00:00 takepiro Exp $'
                              ! 主プログラムのバージョン
                              ! Main program version

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

    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                              ! IOSTAT of NAMELIST read

    integer:: j,k             ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! 実行文 ; 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)
      & )
    !
    ! MPI 計算での Message 制御
    ! Message control for MPI calculation
    !
    call MPIMEssageCntlInit

    ! Initialization of modules used in this module
    !

    ! 時刻管理
    ! Time control
    !
    call TimesetInit

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

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

    ! 格子点設定
    ! Grid points settings
    !
    call GridsetInit

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

    ! 時系列データの読み込み
    ! Reading time series
    !
    call ReadTimeSeriesInit

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level,
    ! and calculate pressure and height
    !
    call AuxVarsInit

    ! 予報変数の割付
    ! Allocation of prediction variables
    !
    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( xyz_Temp   (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_QVap   (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xy_Ps      (0:imax-1, 1:jmax) )

    allocate( yr_Msf    (1:jmax,0:kmax) )
    !
    ! ヒストリデータファイルの初期化
    ! Initialization of history data files
    !
    call HistoryFileOpen
    !
    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'Msf' , &       ! (in)
      & (/ 'lon ', 'lat ', 'sigm', 'time' /), &  ! (in)
      & 'mass stream function', 'kg s-1' )       ! (in)

    !
    ! TEM 計算
    ! TEM calculations
    !
    call TEMOutputInit
    
  end subroutine MainInit

  !-------------------------------------------------------------------

  subroutine MainTerminate
    !
    ! 主プログラムの終了処理手続き. 
    !
    ! Termination procedure for the main program. 
    !

    ! モジュール引用 ; USE statements
    !

    ! MPI
    !
    use mpi_wrapper, only : MPIWrapperFinalize

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: AxessetFinalize

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimesetClose

    ! ヒストリデータ出力
    ! History data output
    !
    use history_file_io, only: HistoryFileClose

    ! 宣言文 ; Declaration statements
    !
    implicit none

    ! 実行文 ; Executable statement
    !

    ! ヒストリデータファイルクローズ
    ! Close history data files
    !
    call HistoryFileClose

    ! 予報変数の割付解除
    ! Deallocation of prediction variables
    !
    deallocate( xyz_U     )
    deallocate( xyz_V     )
    deallocate( xyz_Temp  )
    deallocate( xyz_QVap )
    deallocate( xy_Ps     )

    ! 診断変数の割付解除
    ! Dellocation of diagnostic variables
    !
    deallocate( yr_MSF )

    ! 各モジュール内の変数の割付解除
    ! Dellocation of variables in modules
    !
    call AxessetFinalize

    ! 時刻管理終了処理
    ! Termination of time control
    !
    call TimesetClose

    ! Finalize MPI
    !
    call MPIWrapperFinalize

  end subroutine MainTerminate
  !
  !----------------------------------------------------------------------
  !
  subroutine read_history( time_sec )
    !
    ! ヒストリデータ入力
    ! History data output
    !
    use gtool_history, only : HistoryGet, HistoryGetAttr
    use dc_string,     only : cprintf
    use dc_calendar,   only : DCCalConvertByUnit
    use mpi_wrapper,   only : nprocs, myrank, flag_mpi_init
    
    real(DP), intent(IN) :: time_sec    ! time (sec)

    real(DP)          :: time
    character(STRING) :: Unit
    character(STRING) :: infile_U
    character(STRING) :: infile_V
    character(STRING) :: infile_Temp
    character(STRING) :: infile_SigDot
    character(STRING) :: infile_Ps
    
    if ( flag_mpi_init ) then
       infile_U      = cprintf("U_rank%06d.nc",      i=(/myrank/) )
       infile_V      = cprintf("V_rank%06d.nc",      i=(/myrank/) )
       infile_Temp   = cprintf("Temp_rank%06d.nc",   i=(/myrank/) )
       infile_SigDot = cprintf("SigDot_rank%06d.nc", i=(/myrank/) )
       infile_Ps     = cprintf("Ps_rank%06d.nc",     i=(/myrank/) )
    else
       infile_U      = 'U.nc'
       infile_V      = 'V.nc'
       infile_Temp   = 'Temp.nc'
       infile_SigDot = 'SigDot.nc'
       infile_Ps     = 'Ps.nc'
    endif
    
    call HistoryGetAttr( infile_U, 'time', 'units', Unit )
    Time = DCCalConvertByUnit( time_sec, 'sec', Unit) 

    call HistoryGet( infile_U,      'U',      xyz_U,      time )
    call HistoryGet( infile_V,      'V',      xyz_V,      time )
    call HistoryGet( infile_Temp,   'Temp',   xyz_Temp,   time )
    call HistoryGet( infile_SigDot, 'SigDot', xyr_SigDot, time )
    call HistoryGet( infile_Ps,     'Ps',    xy_Ps,      time )

  end subroutine read_history
    
end program dcpam_calc_tem