!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Module physics_hs94forcing_mod
!
!   * Developers: Morikawa Yasuhiro
!   * Version: $Id$
!   * Tag Name: $Name$
!   * Change History: 
!
!== Overview
!
!Held and Suarez(1994) の乾燥大気 GCM ベンチマーク用の加熱と散逸を計算する。
!
!== Reference
!
!* Held, I. M., Suarez, M. J., 1994: 
!  A proposal for the intercomparison of the dynamical cores of
!  atmospheric general circuation models.
!  Bull. Am. Meteor. Soc., 75, 1825--1830.
!
!== Error Handling
!
!== Known Bugs
!
!* 地球半径や回転角速度、重力加速度などの物理定数が ((<constant_mod>))
!  に依存するようになっている。正しくテストしたいのならば、それらが
!  Held and Suarez(1994) の実験設定と整合的かどうかチェックするように
!  すべきである。
!
!== Note
!
!== Future Plans
!
                                                                 !=end
module physics_hs94forcing_mod
                                                                 !=begin
  !== Dependency
  use type_mod,    only : STRING, DBKIND
                                                                 !=end
  implicit none
                                                                 !=begin
  !== Public Interface
  private
  public:: physics_hs94forcing_init, physics_hs94forcing ! subroutines
                                                                 !=end

  !----------------------------------------------------------------
  !   Multi-Dimensional Data
  !----------------------------------------------------------------
  real(DBKIND), allocatable, save  :: &
       & xyz_Sigma(:,:,:)           , & ! σ
       & xyz_kv(:,:,:)              , & ! k_v
       & xyz_Press(:,:,:)           , & ! 気圧
       & xyz_SinLat(:,:,:)          , & ! sinφ
       & xyz_CosLat(:,:,:)          , & ! cosφ
       & xyz_TempEQ(:,:,:)          , & ! T_eq
       & xyz_kt(:,:,:)                  ! k_t


  logical, save :: physics_hs94forcing_initialized = .false.
  character(STRING),parameter:: version = &
       & '$Id$'
  character(STRING),parameter:: tagname = '$Name:  $'

contains

  !== Procedure Interface
  !
  !=== Initialize module
  !
  !以降のサブルーチンで用いる変数の allocate を行なう。
  !
  !本来はここで、実験設定と ((< constants_mod >)) の値との比較を
  !行ないたいが、まだ実装されていない。
  !(実数同士のそのままの比較は難しいので、有効数字数桁と指数部分を
  !比較するようにすべきだろう)。
  !
  subroutine physics_hs94forcing_init
  !==== Dependency
    use type_mod,      only: INTKIND, STRING, TOKEN, REKIND, DBKIND
    use constants_mod, only: constants_init
    use grid_3d_mod  , only: grid_3d_init, im, jm, km
    use dc_trace     , only: DbgMessage, BeginSub, EndSub
                                                                 !=end
    implicit none


    character(STRING),  parameter:: subname = "physics_hs94forcing_init"
  continue

    !----------------------------------------------------------------
    !   Check Initialization
    !----------------------------------------------------------------
    call BeginSub(subname)
    if (physics_hs94forcing_initialized) then
       call EndSub( subname, '%c is already called.', c1=trim(subname) )
       return
    else
       physics_hs94forcing_initialized = .true.
    endif

    !----------------------------------------------------------------
    !   Version identifier
    !----------------------------------------------------------------
    call DbgMessage('%c :: %c', c1=trim(version), c2=trim(tagname))

    !-------------------------------------------------------------------
    !   Initialize Dependent modules
    !-------------------------------------------------------------------
    call constants_init
    call grid_3d_init

    !-------------------------------------------------------------------
    !   Allocate variables
    !-------------------------------------------------------------------
    allocate( &
         & xyz_Sigma(im, jm, km)  , & ! σ
         & xyz_kv(im, jm, km)     , & ! k_v
         & xyz_Press(im, jm, km)  , & ! 気圧
         & xyz_SinLat(im, jm, km) , & ! sinφ
         & xyz_CosLat(im, jm, km) , & ! cosφ
         & xyz_TempEQ(im,jm,km)   , & ! T_eq
         & xyz_kt(im,jm,km)      )    ! k_t

    call EndSub(subname)
  end subroutine physics_hs94forcing_init

                                                                 !=begin
  !=== 加熱と散逸の計算
  !
  !加熱と散逸の結果は xyz_VelLon_phy, xyz_VelLat_phy, xyz_Temp_phy 
  !として返る。
  !
  subroutine physics_hs94forcing( &
       & x_Lon         , y_Lat         , z_Sigma     , z_SigmaHalf, &
       & xyz_VelLon    , xyz_VelLat    , xyz_Temp    , xy_Ps      , &
       & xyz_VelLon_phy, xyz_VelLat_phy, xyz_Temp_phy            )

  !==== Dependency
    use type_mod,      only: INTKIND, STRING, TOKEN, REKIND, DBKIND
    use constants_mod, only: RAir, Cp, SecPerDay
    use grid_3d_mod  , only: im, jm, km
    use dc_trace     , only: DbgMessage, BeginSub, EndSub
                                                                 !=end
    implicit none
                                                                 !=begin
    !==== Input
    !
    real(DBKIND), intent(in) :: &
         & x_Lon(:)          , & ! 経度座標
         & y_Lat(:)          , & ! 緯度座標
         & z_Sigma(:)        , & ! σレベル(整数)座標
         & z_SigmaHalf(:)    , & ! σレベル(半整数)座標
         &
         & xyz_VelLon(:,:,:) , & ! 速度経度成分 (t)
         & xyz_VelLat(:,:,:) , & ! 速度緯度成分 (t)
         & xyz_Temp(:,:,:)   , & ! 温度         (t)
         & xy_Ps(:,:)            ! 地表面気圧   (t)

    !==== Output
    !
    real(DBKIND), intent(out) :: &
         & xyz_VelLon_phy(:,:,:) , & ! 速度経度成分の加熱散逸効果
         & xyz_VelLat_phy(:,:,:) , & ! 速度緯度成分の加熱散逸効果
         & xyz_Temp_phy(:,:,:)       ! 温度の加熱散逸効果
                                                                 !=end

    !----------------------------------------------------------------
    !   全体で用いる変数
    !----------------------------------------------------------------
    real(DBKIND) :: SigmaB  ! σ_b

    !----------------------------------------------------------------
    !   xyz_VelLon_phy, xyz_VelLat_phy を求めるための変数
    !----------------------------------------------------------------
    real(DBKIND) :: kf     ! k_f

    !----------------------------------------------------------------
    !   xyz_Temp_phy を求めるための変数
    !----------------------------------------------------------------
    real(DBKIND) ::            &
         & Kappa             , & ! κ = Ｒ/Ｃp
         & Press0            , & ! Ｐ_0 (Pa)
         & ka                , & ! k_a
         & ks                , & ! k_s
         & DelTempY          , & ! (ΔT)_y
         & DelPotTempZ           ! (Δθ)_z

    !----------------------------------------------------------------
    !   汎用変数
    !----------------------------------------------------------------
    integer(INTKIND) :: i, j, k

    character(STRING),  parameter:: subname = "physics_hs94forcing"
  continue

    !----------------------------------------------------------------
    !   Check Initialization
    !----------------------------------------------------------------
    call BeginSub(subname)
    if (.not. physics_hs94forcing_initialized) then
       call EndSub( subname, 'Call physics_hs94forcing_init before call %c',  &
            &       c1=trim(subname) )
       return
    endif

    !----------------------------------------------------------------
    !   全体で用いる変数
    !----------------------------------------------------------------
    SigmaB    = 0.7d0

    do k = 1, km
       xyz_Sigma(:,:,k) = z_Sigma(k)
    enddo

    !----------------------------------------------------------------
    !   xyz_VelLon_phy, xyz_VelLat_phy を求める
    !----------------------------------------------------------------
    kf        = 1.0d0 / SecPerDay

    xyz_kv = kf * max(  0.0d0, ( xyz_Sigma - SigmaB ) /( 1.0d0 - SigmaB )  )

    xyz_VelLon_phy = - xyz_kv * xyz_VelLon
    xyz_VelLat_phy = - xyz_kv * xyz_VelLat


    !----------------------------------------------------------------
    !   xyz_Temp_phy を求める
    !----------------------------------------------------------------
    do j = 1, jm
       xyz_SinLat(:,j,:) = sin( y_Lat(j) )
       xyz_CosLat(:,j,:) = cos( y_Lat(j) )
    enddo

    do k = 1, km
       xyz_Press(:,:,k) = z_Sigma(k) / xy_Ps(:,:)
    enddo

    Kappa = RAir / Cp    ! κ=Ｒ／Ｃｐ (気体定数/定圧比熱)
    Press0    = 1000.0d0 * 1.0d2

    ka        = 1.0d0 / ( 40.0d0 * SecPerDay )
    ks        = 1.0d0 / ( 4.0d0 * SecPerDay )
    DelTempY    = 60.0d0
    DelPotTempZ = 10.0d0


    xyz_TempEQ = &
         & max( 200.0d0, &
         &      ( 315.0d0 - DelTempY * xyz_SinLat**2 &
         &        - DelPotTempZ * log( xyz_Press / Press0 ) &
         &                          * xyz_CosLat**2 &
         &      ) &
         &      * ( xyz_Press / Press0 )**Kappa &
         &     )


    xyz_kt = ka + ( ks - ka ) &
         &           * max(  0.0d0, &
         &                   ( xyz_Sigma - SigmaB ) / ( 1.0d0 - SigmaB ) &
         &                ) &
         &           * xyz_CosLat**4

    xyz_Temp_phy = - xyz_kt * ( xyz_Temp - xyz_TempEQ )


    call EndSub(subname)
  end subroutine physics_hs94forcing

end module physics_hs94forcing_mod
