!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Module hs94forcing_mod
!
!   * Developers: Morikawa Yasuhiro
!   * Version: $Id: hs94forcing.f90,v 1.9 2005/01/22 09:23:41 morikawa Exp $
!   * Tag Name: $Name: dcpam2-20050522-1 $
!   * Change History: 
!
!== Overview
!
!This module compute heating and dissipation 
!for Held and Suarez(1994) benchmark integration of a dry GCM.
!
!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 hs94forcing_mod
                                                                 !=begin
  !== Dependency
  use type_mod,    only : STRING, DBKIND
                                                                 !=end
  implicit none
                                                                 !=begin
  !== Public Interface
  private
  public:: hs94forcing_init, hs94forcing, hs94forcing_end  ! subroutines
                                                                 !=end

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


  logical, save :: hs94forcing_initialized = .false.
  character(STRING),parameter:: version = &
       & '$Id: hs94forcing.f90,v 1.9 2005/01/22 09:23:41 morikawa Exp $'
  character(STRING),parameter:: tagname = '$Name: dcpam2-20050522-1 $'

contains

  !== Procedure Interface
  !
  !=== Initialize module
  !
  !以降のサブルーチンで用いる変数の allocate を行なう。
  !
  !本来はここで、実験設定と ((< constants_mod >)) の値との比較を
  !行ないたいが、まだ実装されていない。
  !(実数同士のそのままの比較は難しいので、有効数字数桁と指数部分を
  !比較するようにすべきだろう)。
  !
  subroutine hs94forcing_init( &
       & x_Lon         , y_Lat         , z_Sigma )

  !==== Dependency
    use type_mod,      only: INTKIND, STRING, TOKEN, REKIND, DBKIND
    use constants_mod, only: constants_init, pi
    use grid_3d_mod  , only: grid_3d_init, im, jm, km
    use axis_type_mod, only: AXISINFO
    use io_gt4_out_mod,only: io_gt4_out_init, io_gt4_out_SetVars
    use dc_string    , only: LChar, StrHead
    use dc_trace     , only: DbgMessage, BeginSub, EndSub, DataDump
                                                                 !=end
    implicit none
                                                                 !=begin
    !==== Input
    !
    type(AXISINFO), intent(in) :: &
         & x_Lon                , & ! 経度座標
         & y_Lat                , & ! 緯度座標
         & z_Sigma                  ! σレベル(整数)座標
                                                                 !=end
    integer(INTKIND) :: i, j, k
    real(DBKIND)     :: RadDegFact     ! ラジアンと度数の変換係数

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

    !----------------------------------------------------------------
    !   Check Initialization
    !----------------------------------------------------------------
    call BeginSub(subname)
    if (hs94forcing_initialized) then
       call EndSub( subname, '%c is already called.', c1=trim(subname) )
       return
    else
       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
    call io_gt4_out_init

    !-------------------------------------------------------------------
    !   Setting Ouput Data by io_gt4_out_set_Vars
    !-------------------------------------------------------------------
    call io_gt4_out_SetVars('xyz_VelLon_phy')
    call io_gt4_out_SetVars('xyz_VelLat_phy')
    call io_gt4_out_SetVars('xyz_Temp_phy')
    call io_gt4_out_SetVars('xyz_TempEQ')
    call io_gt4_out_SetVars('xyz_PotTempEQ')
    call io_gt4_out_SetVars('1.0/xyz_kt')
    call io_gt4_out_SetVars('1.0/max(xyz_kv, 1.0/86400d2)/86400.0')
    call io_gt4_out_SetVArs('xyz_kv')


    !-----------------------------------------------------------------
    !   Allocate variables
    !-----------------------------------------------------------------
    allocate( &
         & xyz_Sigma(im, jm, km)    , & ! σ
         & xyz_kv(im, jm, km)       , & ! k_v
         & xyz_Press_b(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

    !----------------------------------------------------------------
    !   全体で用いる変数
    !----------------------------------------------------------------
    do k = 1, km
       xyz_Sigma(:,:,k) = z_Sigma%a_Dim(k)
    enddo

    !----------------------------------------------------------------
    !   xyz_Temp_phy を求めるための変数
    !----------------------------------------------------------------
    if (  StrHead( 'radians', trim(LChar(y_Lat%axisinfo%units)) ) .or.&
         & StrHead( 'rad.', trim(LChar(y_Lat%axisinfo%units)) ) ) then
       RadDegFact = 1.
    else
       RadDegFact = pi / 180.
    end if

    do j = 1, jm
       xyz_SinLat(:,j,:) = sin( y_Lat%a_Dim(j) * RadDegFact )
       xyz_CosLat(:,j,:) = cos( y_Lat%a_Dim(j) * RadDegFact )
    enddo

    call DataDump('xyz_SinLat', xyz_SinLat, strlen=70)
    call DataDump('xyz_CosLat', xyz_CosLat, strlen=70)

    call EndSub(subname)
  end subroutine hs94forcing_init

                                                                 !=begin
  !=== 加熱と散逸の計算
  !
  !加熱と散逸の結果は xyz_VelLon_phy, xyz_VelLat_phy, xyz_Temp_phy 
  !として返る。
  !
  subroutine hs94forcing( &
       & xyz_VelLon_b  , xyz_VelLat_b  , xyz_Temp_b  , xy_Ps_b    , &
       & 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 io_gt4_out_mod,only: io_gt4_out_Put
    use dc_trace     , only: DbgMessage, BeginSub, EndSub, DataDump
                                                                 !=end
    implicit none
                                                                 !=begin
    !==== Input
    !
    real(DBKIND), intent(in) :: &
         & xyz_VelLon_b(:,:,:) , & ! 速度経度成分 (t-Δt)
         & xyz_VelLat_b(:,:,:) , & ! 速度緯度成分 (t-Δt)
         & xyz_Temp_b(:,:,:)   , & ! 温度         (t-Δt)
         & xy_Ps_b(:,:)            ! 地表面気圧   (t-Δ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 = "hs94forcing"
  continue

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

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

    !----------------------------------------------------------------
    !   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_b
    xyz_VelLat_phy = - xyz_kv * xyz_VelLat_b

    call DbgMessage('kf=<%f>', d=(/kf/))
    call io_gt4_out_Put('xyz_VelLon_phy', xyz_VelLon_phy)
    call io_gt4_out_Put('xyz_VelLat_phy', xyz_VelLat_phy)
    call io_gt4_out_Put('xyz_kv', xyz_kv)

    !----------------------------------------------------------------
    !   xyz_Temp_phy を求める
    !----------------------------------------------------------------
    do k = 1, km
       xyz_Press_b(:,:,k) = xyz_Sigma(:,:,k) * xy_Ps_b(:,:)
    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_b / Press0 ) &
         &                          * xyz_CosLat**2 &
         &      ) &
         &      * ( xyz_Press_b / 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_b - xyz_TempEQ )

    call io_gt4_out_Put('xyz_Temp_phy', xyz_Temp_phy)

    call io_gt4_out_Put('xyz_TempEQ', xyz_TempEQ)
    call io_gt4_out_Put('xyz_PotTempEQ', &
         & xyz_TempEQ * ( xyz_Sigma )**( -Kappa )  )

!!$    call io_gt4_out_Put('1.0/max(xyz_kv, 1.0/86400d2)/86400.0', &
!!$         & 1.0/max(xyz_kv, 1.0/86400d2)/86400.0 )
!!$    call io_gt4_out_Put('1.0/xyz_kt', 1.0/xyz_kt)

!!$    call DataDump('1.0/xyz_kt', 1.0/xyz_kt, strlen=60)
!!$
!!$    call DataDump('xyz_Press_b', xyz_Press_b, strlen=60)
!!$
!!$    call DataDump('xyz_TempEQ1', &
!!$         &      ( 315.0d0 - DelTempY * xyz_SinLat**2 &
!!$         &        - DelPotTempZ * log( xyz_Press_b / Press0 ) &
!!$         &                          * xyz_CosLat**2 &
!!$         &      ) &
!!$         &      * ( xyz_Press_b / Press0 )**Kappa  &
!!$         & , strlen=60 )
!!$
!!$    call DataDump('xyz_TempEQ2', &
!!$         &       - DelTempY * xyz_SinLat**2 &
!!$         & , strlen=60 )
!!$
!!$    call DataDump('xyz_TempEQ3', &
!!$         &      ( - DelPotTempZ * log( xyz_Press_b / Press0 ) &
!!$         &                          * xyz_CosLat**2 &
!!$         &      ) &
!!$         & , strlen=60 )

    call EndSub(subname)
  end subroutine hs94forcing



                                                                 !=begin
  !=== Terminate module
  !
  !((< hs94forcing_init>)) で allocate した変数を deallocate する。
  !
  subroutine hs94forcing_end
  !==== Dependency
    use type_mod, only: STRING, DBKIND, INTKIND
    use dc_trace, only: BeginSub, EndSub, DbgMessage
                                                                 !=end
    implicit none

    !-----------------------------------------------------------------
    !   変数定義
    !-----------------------------------------------------------------
    !----- 作業用内部変数 -----
    character(STRING),  parameter:: subname = "hs94forcing_end"

  continue

    !-----------------------------------------------------------------
    !   Check Initialization
    !-----------------------------------------------------------------
    call BeginSub(subname)
    if ( .not. hs94forcing_initialized) then
       call EndSub( subname, 'hs94forcing_init was not called', &
            &       c1=trim(subname) )
       return
    else
       hs94forcing_initialized = .false.
    endif

    !-----------------------------------------------------------------
    !   Deallocate Variables
    !-----------------------------------------------------------------
    deallocate( &
         & xyz_Sigma    , & ! σ
         & xyz_kv       , & ! k_v
         & xyz_Press_b  , & ! 気圧
         & xyz_SinLat   , & ! sinφ
         & xyz_CosLat   , & ! cosφ
         & xyz_TempEQ   , & ! T_eq
         & xyz_kt      )    ! k_t


    call EndSub(subname)
  end subroutine hs94forcing_end


end module hs94forcing_mod
