!---------------------------------------------------------------------
!   Copyright (C) GFD Dennou Club, 2005.  All rights reserved
!---------------------------------------------------------------------
!=begin
!= Module ((* physics_hs94forcing_mod *))
!
!== Overview
!
!"Held-Suarez" の乾燥大気 GCM ベンチマーク用の加熱と散逸を計算する。
!
!=end

module physics_hs94forcing_mod

  !=begin
  !== Dependency
  use type_mod,    only : INTKIND, STRING, TOKEN, REKIND, DBKIND
  !=end

  private
  !=begin
  != Public Interface
  public:: physics_hs94forcing_init, physics_hs94forcing
  !=end

!    use physics_hs94forcing_mod [, only: physics_hs94forcing_init, &
!         &                               physics_hs94forcing,      ]
!
!: physics_hs94forcing_init
!    モジュールを初期化します。
!
!: physics_hs94forcing
!    Tendency を返します。
!
!
!= PUBLIC ROUTINE
!
!    call physics_hs94forcing_init
!
!    call physics_hs94forcing(day, kappa, im, jm, km, sinlat, coslat, &
!         gup, gvp, gtp, gphp, gpp, gduphy, gdvphy, gdtphy 
!
!: level
!  * 文字型変数 or 整数型変数。
!    * 'W' : 警告を出力。
!    * 'M' : メッセージを出力。
!    * 'E' : StoreError(USR_ECHAR, ...) を呼ぶ。
!    * その他の文字 : 'M' と同じ動作。
!    * 数値: StoreError にそのまま渡す。
!
!: where
!  * 文字型変数
!    * サブルーチンの名前
!
!: cause_c
!  * 文字型変数
!    * 出力するメッセージ
!
!= NAMELIST
!
!
!= ERROR
!
!= NOTE
!
!= FUTURE PLANS


  subroutine physics_hs94forcing_init
    use grid_3d_mod, only : grid_3d_init

  end subroutine physics_hs94forcing_init

  subroutine physics_hs94forcing( day, kappa, sinlat, coslat, gup&
       &, gvp, gtp, gphp, gpp, gduphy, gdvphy, gdtphy )

    use type_mod,    only : INTKIND, STRING, TOKEN, REKIND, DBKIND
    use grid_3d_mod, only : im, jm, km

    implicit none

    real(DBKIND), intent(in) :: day       ! １日間の秒数
    real(DBKIND), intent(in) :: kappa     ! κ = Ｒ/Ｃp
    real(DBKIND), intent(in) :: sinlat( jm ), coslat( jm )
    real(DBKIND), intent(in) :: gup ( im, jm, km   ), gvp( im, jm, km ), gtp( im, jm, km )
    real(DBKIND), intent(in) :: gphp( im, jm, km+1 ), gpp( im, jm, km )
    real(DBKIND), intent(out) :: gduphy( im, jm, km ), gdvphy( im, jm, km ), gdtphy( im, jm, km )


    ! Local variables
    !
    real(DBKIND) :: gps( im, jm )

    integer(INTKIND) :: i, j, k


    gps( :, : ) = gphp( :, :, km+1 )

    call hs94tforcing( day, kappa, im, jm, km, sinlat, coslat, gpp, gps, &
         gtp, gdtphy )

    call hs94uvforcing( day, im, jm, km, gpp, gps, gup, gvp, gduphy, gdvphy )

  end subroutine hs94forcing

    !**************************************************************************

  subroutine hs94tforcing( day, kappa, im, jm, km, sinlat, coslat, &
       gpp, gps, gt, gdtphy )

    use vtype

    implicit none

    real(DBKIND), intent(in) :: day       ! １日間の秒数
    real(DBKIND), intent(in) :: kappa     ! κ = Ｒ/Ｃp
    integer(INTKIND), intent(in) :: im, jm, km
    real(DBKIND), intent(in) :: sinlat( jm ), coslat( jm )
    real(DBKIND), intent(in) :: gpp( im, jm, km ), gps( im, jm )
    real(DBKIND), intent(in) :: gt( im, jm, km )
    real(DBKIND), intent(out) :: gdtphy( im, jm, km )


    ! Local variables
    !
    real(DBKIND) :: gteq( im, jm, km )
    real(DBKIND) :: p0, kt, ka, ks, sig, sigb, dty, dthetaz

    integer(INTKIND) :: i, j, k


    p0 = 1000.0d0 * 1.0d2

    sigb    = 0.7d0
    ka      = 1.0d0 / ( 40.0d0 * day )
    ks      = 1.0d0 / ( 4.0d0 * day )
    dty     = 60.0d0
    dthetaz = 10.0d0

    do k = 1, km
       do j = 1, jm
          do i = 1, im
             gteq( i, j, k ) = max( 200.0d0, &
                  ( 315.0d0 - dty * sinlat( j ) * sinlat( j ) &
                  - dthetaz * log( gpp( i, j, k ) / p0 ) &
                  * coslat( j ) * coslat( j ) ) &
                  *( gpp( i, j, k ) / p0 )**kappa )
          enddo
       enddo
    enddo

    do k = 1, km
       do j = 1, jm
          do i = 1, im
             sig = gpp( i, j, k ) / gps( i, j )
             kt = ka &
                  + ( ks - ka ) &
                  * max( 0.0d0, &
                  (sig-sigb)/(1.0d0-sigb) ) &
                  *coslat(j)*coslat(j)*coslat(j)*coslat(j)

             gdtphy( i, j, k ) = -kt * ( gt( i, j, k ) - gteq( i, j, k ) )
          enddo
       enddo
    enddo


  end subroutine hs94tforcing

    !**************************************************************************

  subroutine hs94uvforcing( day, im, jm, km, gpp, gps, gu, gv, &
       gduphy, gdvphy )

    use vtype

    implicit none

    real(DBKIND), intent(in) :: day
    integer(INTKIND), intent(in) :: im, jm, km
    real(DBKIND), intent(in) :: gpp( im, jm, km ), gps( im, jm )
    real(DBKIND), intent(in) :: gu( im, jm, km ), gv( im, jm, km )
    real(DBKIND), intent(out) :: gduphy( im, jm, km ), gdvphy( im, jm, km )


    ! Local variables
    !
    real(DBKIND) :: p0, kv, kf, sig, sigb

    integer(INTKIND) :: i, j, k


    p0 = 1000.0d0 * 1.0d2

    sigb = 0.7d0
    kf   = 1.0d0 / day

    do k = 1, km
       do j = 1, jm
          do i = 1, im
             sig = gpp( i, j, k ) / gps( i, j )
             kv = kf * max( 0.0d0, (sig-sigb)/(1.0d0-sigb) )

             gduphy( i, j, k ) = - kv * gu( i, j, k )
             gdvphy( i, j, k ) = - kv * gv( i, j, k )
          enddo
       enddo
    enddo


  end subroutine hs94uvforcing

end module physics_hs94forcing_mod
