!= thermal_forcing_mod
!
! This file provides thermal_forcing module.
!
module thermal_forcing_mod

!== Dependency
use type_mod,      only : STRING, DBKIND, INTKIND
use gt4_history, only:  HistoryAddVariable, GT_HISTORY_AXIS
implicit none

!== Public Interface
private
public :: thermal_forcing_sub     !subroutine

contains

subroutine thermal_forcing_sub
  !== Dependency
  use type_mod,      only : STRING, DBKIND, INTKIND
  use nmlfile_mod,   only: nmlfile_init, nmlfile_open, nmlfile_close
  use gt4_history, only:  HistoryAddVariable, GT_HISTORY_AXIS
  use constants_mod  , only: constants_init, pi
  use gt4_history,     only: HistoryAxisInquire, HistoryAxisCreate
  use gt4_history, only: HistoryCreate, HistoryPut, HistoryAddVariable, &
      & GT_HISTORY, GT_HISTORY_AXIS
  use axis_type_mod, only : AXISINFO
  use time_mod,    only : InitTime, DelTime, StepInterval, &
      & tvar, ttype, tname, tunit
  use dycore_type_mod, only: DYCORE_DIMS
  use dycore_grid_mod, only: dycore_grid_init, im, jm, km, nm
  implicit none

  !!==  変数を定義する ==============================================
    real(DBKIND)      :: LonLat_Radius_Rad = 0.174  ! 半径 (ラジアン)
    real(DBKIND)      :: Lat_Center_Rad = 0.0 ! 緯度の中心位置 (ラジアン)
    real(DBKIND)      :: Lon_Center_Rad = 3.14 ! 経度の中心位置 (ラジアン)
    real(DBKIND)      :: Thermal_forcing = 10.0d0 / 86400.0d0 ! dT/dt (Max value)
    real(DBKIND)      :: RadDegFact     ! ラジアンを度数に変換する係数
    real(DBKIND), allocatable :: xyz_Lon(:,:,:) ! 経度座標 (ラジアン表記)
    real(DBKIND), allocatable :: xyz_Lat(:,:,:) ! 緯度座標 (ラジアン表記)
    real(DBKIND), save, allocatable :: xyz_Temp_phy(:,:,:)       ! 温度の加熱散逸効果
    character(STRING) :: axis_units

  !!#== gt4io =============================================================
    type(GT_HISTORY_AXIS) :: axes(1:4)
    axes(1) = Dims % x_Lon % axisinfo
    axes(2) = Dims % y_Lat % axisinfo

    !-----------------------------------------------------------------
    !  Initialize Dependent Modules
    !-----------------------------------------------------------------
    call constants_init
    call dycore_grid_init
    call hs94forcing_init(   &
         & Dims%x_Lon      , & ! intent(in): 経度座標
         & Dims%y_Lat      , & ! intent(in): 緯度座標
         & Dims%z_Sigma    )   ! intent(in): σレベル(整数)座標

  namelist /init_nml/ &
       & LonLat_Radius_Deg   , & ! 半径 (度数)
       & LonLat_Radius_Rad   , & ! 半径 (ラジアン)
       & Lat_Center_Deg      , & ! 緯度の中心位置 (度数)
                                                                 !=begin
  !==== Input
  !
    type(DYCORE_DIMS), intent(in)   :: Dims   ! 次元データ全種
                                                                 !=end
    !----- 作業用内部変数 -----
    integer :: i, j
    type(GT_HISTORY_AXIS) :: axes(1:4)


    !-----------------------------------------------------------------
    !   Allocate variables
    !-----------------------------------------------------------------
    allocate( xyz_Temp_phy(im,jm,km) )    ! 温度の加熱散逸効果
    allocate( xyz_Lon(im,jm,km) )
    allocate( xyz_Lat(im,jm,km) )


    !---- 度 or ラジアン ---------------------------------------------
    call HistoryAxisInquire(Dims%x_Lon % axisinfo, units=axis_units)
    if (  StrHead( 'radians', trim(LChar(axis_units)) ) .or.&
      &   StrHead( 'rad.', trim(LChar(axis_units)) ) ) then
      RadDegFact = 1.
    else
      RadDegFact = 180. / pi
    end if

    do i = 1, im
      xyz_Lon(i,:,:) = Dims%x_Lon%a_Dim(i) / RadDegFact
    end do

    call HistoryAxisInquire(Dims%y_Lat % axisinfo, units=axis_units)
    if (  StrHead( 'radians', trim(LChar(axis_units)) ) .or.&
      &   StrHead( 'rad.', trim(LChar(axis_units)) ) ) then
      RadDegFact = 1.
    else
      RadDegFact = 180. / pi
    end if

    do j = 1, jm
      xyz_Lat(:,j,:) = Dims%y_Lat%a_Dim(j) / RadDegFact
    end do


    !---- Thermal_forcing ---------------------------------------------
    Thermal_forcing = 10.0d0 / 86400.0d0

    xyz_Temp_phy = Thermal_forcing &
      &            *  ( 1.0D0 &
      &                 - sqrt(  &
      &                        min(  LonLat_Radius_Rad**2 ,              &
      &                              ( xyz_Lon - Lon_Center_Rad )**2     &
      &                                + ( xyz_Lat - Lat_Center_Rad )**2 &
      &                           )     &
      &                       ) / LonLat_Radius_Rad  &
      &               )

end subroutine thermal_forcing_sub
end module thermal_forcing_mod
