!= Held and Suarez (1994) による动扩と欢帮
!
!= Forcing and dissipation suggested by Held and Suarez (1994)
!
! Authors::   Yasuhiro MORIKAWA
! Version::   $Id: held_suarez_1994.f90,v 1.1.1.1 2010-08-17 05:24:51 takepiro Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module held_suarez_1994
  !
  != Held and Suarez (1994) による动扩と欢帮
  !
  != Forcing and dissipation suggested by Held and Suarez (1994)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! Held and Suarez (1994) で捏捌される触羚络丹 GCM ベンチマ〖ク脱の
  ! 动扩と欢帮を纷换します.
  ! 涂える动扩と欢帮として, 补刨眷の掠觉滦疚眷への词帽なニュ〖トン武笛と,
  ! 董肠霖肃护を山附する布霖慎のレイリ〖肃护を脱います.
  ! 拒嘿を笆布に淡します.
  !
  ! Forcing and dissipation for dry air GCM benchmark
  ! suggested by Held and Suarez (1994) are caluclate.
  ! We use simple Newtonian relaxation of the temperature field to a
  ! zonally symmetric state and Rayleigh damping of low-level winds to
  ! represent boundary-layer friction.
  ! Their specifications are detailed as follows.
  !
  ! \[
  !    \left( \DP{\Dvect{v}}{t} \right)_{\mathrm{HS94}} =
  !        - k_v (\sigma) \Dvect{v}, \] \[
  !    \left( \DP{T}{t} \right)_{\mathrm{HS94}} =
  !        - k_T (\phi, \sigma) [T - T_{eq} (\phi,p)], \] \[
  !    T_{eq} = \mathrm{max}
  !     \left\{
  !        200 \mathrm{K},
  !        \left[
  !          315 \mathrm{K} - (\Delta T)_y \sin^2\phi
  !                         - (\Delta \theta)_z
  !                           \log \left(\frac{p}{p_0}\right) \cos^2\phi
  !        \right] \left(\frac{p}{p_0}\right)^\kappa
  !     \right\}, \] \[
  !    k_T = k_a + (k_s - k_a)
  !          \mathrm{max}
  !          \left(0, \frac{\sigma - \sigma_b}{1 - \sigma_b}\right) \cos^4\phi,
  !     \] \[
  !    k_v = k_f
  !          \mathrm{max}
  !          \left(0, \frac{\sigma - \sigma_b}{1 - \sigma_b}\right),
  !     \] \[
  !    \sigma_b = 0.7, \qquad
  !    k_f = 1 \mathrm{day}^{-1}, \qquad
  !    k_a = \Dinv{40} \mathrm{day}^{-1}, \qquad
  !    k_s = \Dinv{4} \mathrm{day}^{-1}, \] \[
  !    (\Delta T)_y = 60 \mathrm{K}, \qquad
  !    (\Delta \theta)_z = 10 \mathrm{K}, \qquad
  !    p_0 = 1000 \mathrm{hPa}, \qquad
  !    \kappa = \frac{R}{c_p}.
  ! \]
  !
  ! Forcing では, 涂えられた庐刨や补刨 ( $ t+\Delta t$ を鳞年)
  ! に滦して笆布のように动扩と欢帮を努脱します.
  !
  ! By Forcing, forcing and dissipation are applied to
  ! given wind and temperature ($ t+\Delta t$ is expected) as follows.
  !
  ! \[
  !    \hat{\Dvect{v}}^{t+\Delta t} =
  !      \Dvect{v}^{t+\Delta t}
  !      + 2 \Delta t \left( \DP{\Dvect{v}}{t} \right)_{\mathrm{HS94}} \] \[
  !    \hat{T}^{t+\Delta t} =
  !      T^{t+\Delta t}
  !      + 2 \Delta t \left( \DP{T}{t} \right)_{\mathrm{HS94}}
  ! \]
  !
  !== Procedures List
  !
  ! Hs94Forcing   :: 动扩と欢帮の纷换
  ! Hs94Finalize  :: 姜位借妄 (モジュ〖ル柒婶の恃眶の充り烧け豺近)
  ! ------------  :: ------------
  ! Hs94Forcing   :: Calculate forcing and dissipation
  ! Hs94Finalize  :: Termination (deallocate variables in this module)
  !
  !--
  !== NAMELIST
  !
  ! NAMELIST#held_suarez_1994_nml
  !++
  !== References
  !
  ! * Held, I. M., Suarez, M. J., 1994:
  !   A proposal for the intercomparison of the dynamical cores of
  !   atmospheric general circuation models.
  !   <i>Bull. Am. Meteor. Soc.</i>, <b>75</b>, 1825--1830.
  !

  ! モジュ〖ル苞脱 ; USE statements
  !

  ! 呈灰爬肋年
  ! 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

  ! 硷侍房パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 擒篮刨悸眶房. Double precision.
    &                 STRING     ! 矢机误.       Strings.

  ! メッセ〖ジ叫蜗
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 离咐矢 ; Declaration statements
  !
  implicit none
  private

  ! 给倡缄鲁き
  ! Public procedure
  !
  public:: Hs94Forcing, Hs94Finalize

  ! 给倡恃眶
  ! Public variables
  !
  logical, save, public:: held_suarez_1994_inited = .false.
                              ! 介袋肋年フラグ.
                              ! Initialization flag

  ! 润给倡恃眶
  ! Private variables
  !

  real(DP):: Kappa            ! $ \kappa = R/C_p $ .
                              ! 丹挛年眶の年暗孺钱に滦する孺. Ratio of gas constant to specific heat
  real(DP):: P0               ! $ p_0 $ .
  real(DP):: DelTempY         ! $ (\Delta T)_y $ .
  real(DP):: DelPotTempZ      ! $ (\Delta \theta)_z $ .
  real(DP), allocatable:: z_kv (:)
                              ! $ k_v $ .
  real(DP), allocatable:: yz_kt (:,:)
                              ! $ k_T $ .

  character(*), parameter:: module_name = 'held_suarez_1994'
                              ! モジュ〖ルの叹疚.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: held_suarez_1994.f90,v 1.1.1.1 2010-08-17 05:24:51 takepiro Exp $'
                              ! モジュ〖ルのバ〖ジョン
                              ! Module version

  ! INTERFACE 矢 ; INTERFACE statements
  !
  interface Hs94Forcing
    module procedure Hs94Forcing
  end interface

  interface Hs94Finalize
    module procedure Hs94Finalize
  end interface

contains

  subroutine Hs94Forcing( &
    & xyz_U,    xyz_V,    xyz_Temp, xy_Ps, & ! (in)
    & xyz_DUDt, xyz_DVDt, xyz_DTempDt &      ! (out)
    & )
    !
    ! 苞眶として涂えられた澎谰慎庐 xyz_U, 祁颂慎庐 xyz_V,
    ! 补刨 xyz_Temp から,
    ! 补刨眷の掠觉滦疚眷への词帽なニュ〖トン武笛と
    ! 董肠霖肃护を山附する布霖慎のレイリ〖肃护による
    ! 慎庐と补刨の恃步唯を滇め,
    ! xyz_DUDt, xyz_DVDt, xyz_DTempDt に手します.
    !
    ! Tendencies by simple Newtonian relaxation of the temperature field to a
    ! zonally symmetric state and Rayleigh damping of low-level winds to
    ! represent boundary-layer friction are calculated
    ! from eastward wind "xyz_U", northward wind "xyz_V",
    ! temperature "xyz_Temp".
    ! And the tencencies are returned as
    ! "xyz_DUDt", "xyz_DVDt", "xyz_DTempDt".
    !
    !

    ! モジュ〖ル苞脱 ; USE statements
    !

    ! 郝筛デ〖タ肋年
    ! Axes data settings
    !
    use axesset, only: &
      & y_Lat, &              ! $ \varphi $ [rad.] . 稗刨. Latitude
      & z_Sigma               ! $ \sigma $ レベル (腊眶).
                              ! Full $ \sigma $ level

    ! 箕癸瓷妄
    ! Time control
    !
    use timeset, only: &
      & DelTime, &            ! $ \Delta t $
      & TimeN, &              ! ステップ $ t $ の箕癸. Time of step $ t $.
      & TimesetClockStart, TimesetClockStop

    ! ヒストリデ〖タ叫蜗
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 离咐矢 ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u $ . 澎谰慎庐.
                              ! Eastward wind
    real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v $ . 祁颂慎庐.
                              ! Northward wind
    real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ . 补刨.
                              ! Temperature
    real(DP), intent(in):: xy_Ps (0:imax-1, 1:jmax)
                              ! $ p_s $ . 孟山烫丹暗.
                              ! Surface pressure
    real(DP), intent(out):: xyz_DUDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{u}{t} $ . 澎谰慎庐恃步.
                              ! Eastward wind tendency
    real(DP), intent(out):: xyz_DVDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{v}{t} $ . 祁颂慎庐恃步.
                              ! Northward wind tendency
    real(DP), intent(out):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{T}{t} $ . 补刨恃步.
                              ! Temperature tendency

    ! 侯度恃眶
    ! Work variables
    !
    real(DP):: xyz_TempEQ (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T_{eq} $ . 士拐补刨.
                              ! Equilibrium temperature
    real(DP):: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ . 暗蜗.
                              ! Pressure

    integer:: j               ! 稗刨数羹に搀る DO ル〖プ脱侯度恃眶
                              ! Work variables for DO loop in latitude
    integer:: k               ! 濒木数羹に搀る DO ル〖プ脱侯度恃眶
                              ! Work variables for DO loop in vertical direction

    ! 悸乖矢 ; Executable statement
    !

    ! 纷换箕粗纷卢倡幌
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )

    ! 介袋步
    ! Initialization
    !
    if ( .not. held_suarez_1994_inited ) call Hs94Init

    ! 澎谰慎庐 $ u $ と祁颂慎庐 $ v $ へレイリ〖肃护を努脱
    ! Apply Rayleigh damping to eastward wind $ u $ and northward wind $ v $
    !
    do k = 1, kmax
      xyz_DUDt (:,:,k) = - z_kv (k) * xyz_U (:,:,k)
      xyz_DVDt (:,:,k) = - z_kv (k) * xyz_V (:,:,k)
    end do

    ! 补刨 $ T $ へニュ〖トン武笛を努脱
    ! Apply Newtonian relaxation to temperature $ T $
    !
    do k = 1, kmax
       xyz_Press(:,:,k) = z_Sigma(k) * xy_Ps
    enddo

    do j = 1, jmax
      xyz_TempEQ(:,j,:) = &
        & max( 200.0_DP, &
        &      (   315.0_DP &
        &        - DelTempY * sin( y_Lat(j) ) ** 2 &
        &        - DelPotTempZ * log( xyz_Press(:,j,:) / P0 ) &
        &                                 * cos( y_Lat(j) ) ** 2 &
        &      ) &
        &      * ( xyz_Press(:,j,:) / P0 ) ** Kappa &
        &     )
    end do

    do k = 1, kmax
      do j = 1, jmax
        xyz_DTempDt (:,j,k) = &
          & - yz_kt (j,k) * ( xyz_Temp (:,j,k) - xyz_TempEQ (:,j,k) )
      end do
    end do

    ! ヒストリデ〖タ叫蜗
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'DUDtHs94',    xyz_DUDt )
    call HistoryAutoPut( TimeN, 'DVDtHs94',    xyz_DVDt )
    call HistoryAutoPut( TimeN, 'DTempDtHs94', xyz_DTempDt )
    call HistoryAutoPut( TimeN, 'TempEQHs94',  xyz_TempEQ )

    ! 纷换箕粗纷卢办箕匿贿
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine Hs94Forcing

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

  subroutine Hs94Init
    !
    ! held_suarez_1994 モジュ〖ルの介袋步を乖います.
    ! NAMELIST#held_suarez_1994_nml の粕み哈みはこの缄鲁きで乖われます.
    !
    ! "held_suarez_1994" module is initialized.
    ! "NAMELIST#held_suarez_1994_nml" is loaded in this procedure.
    !

    ! モジュ〖ル苞脱 ; USE statements
    !

    ! 湿妄年眶肋年
    ! Physical constants settings
    !
    use constants, only: &
      & 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

    ! 郝筛デ〖タ肋年
    ! Axes data settings
    !
    use axesset, only: &
      & y_Lat, &              ! $ \varphi $ [rad.] . 稗刨. Latitude
      & z_Sigma               ! $ \sigma $ レベル (腊眶).
                              ! Full $ \sigma $ level

    ! NAMELIST ファイル掐蜗に簇するユ〖ティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! ファイル掐叫蜗输锦
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! 硷侍房パラメタ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT ! 筛洁叫蜗の刘弥戎规. Unit number of standard output

    ! 矢机误拎侯
    ! Character handling
    !
    use dc_string, only: StoA

    ! ヒストリデ〖タ叫蜗
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoAddVariable

    ! 离咐矢 ; Declaration statements
    !
    implicit none

    real(DP):: day_seconds
                              ! 1 泣の擅眶.
                              ! Seconds in day.
    real(DP):: SigmaB         ! $ \sigma_b $ .
    real(DP):: kf             ! $ k_f $ .
    real(DP):: ka             ! $ k_a $ .
    real(DP):: ks             ! $ k_s $ .

    integer:: j               ! 稗刨数羹に搀る DO ル〖プ脱侯度恃眶
                              ! Work variables for DO loop in latitude
    integer:: k               ! 濒木数羹に搀る DO ル〖プ脱侯度恃眶
                              ! Work variables for DO loop in vertical direction

!!$    integer:: unit_nml        ! NAMELIST ファイルオ〖プン脱刘弥戎规.
!!$                              ! Unit number for NAMELIST file open
!!$    integer:: iostat_nml      ! NAMELIST 粕み哈み箕の IOSTAT.
!!$                              ! IOSTAT of NAMELIST read

    ! NAMELIST 恃眶凡
    ! NAMELIST group name
    !
!!$    namelist /held_suarez_1994_nml/
!!$          !
!!$          ! デフォルト猛については介袋步缄鲁 "held_suarez_1994#Hs94Init"
!!$          ! のソ〖スコ〖ドを徊救のこと.
!!$          !
!!$          ! Refer to source codes in the initialization procedure
!!$          ! "held_suarez_1994#Hs94Init" for the default values.
!!$          !

    ! 悸乖矢 ; Executable statement
    !

    if ( held_suarez_1994_inited ) return
    call InitCheck

    ! デフォルト猛の肋年
    ! Default values settings
    !

    ! NAMELIST の粕み哈み
    ! NAMELIST is input
    !
!!$    if ( trim(namelist_filename) /= '' ) then
!!$      call FileOpen( unit_nml, &          ! (out)
!!$        & namelist_filename, mode = 'r' ) ! (in)
!!$
!!$      rewind( unit_nml )
!!$      read( unit_nml, &                ! (in)
!!$        & nml = held_suarez_1994_nml, &  ! (out)
!!$        & iostat = iostat_nml )        ! (out)
!!$      close( unit_nml )
!!$
!!$      call NmlutilMsg( iostat_nml, module_name ) ! (in)
!!$!$      if ( iostat_nml == 0 ) write( STDOUT, nml = held_suarez_1994_nml )
!!$    end if

    ! 犯眶の肋年
    ! Configure coefficients
    !
    Kappa = GasRDry / CpDry

    P0          = 1000.0e2_DP
    DelTempY    = 60.0_DP
    DelPotTempZ = 10.0_DP

    day_seconds = 86400.0_DP

    SigmaB = 0.7_DP
    kf     = 1.0_DP / day_seconds
    ka     = 1.0_DP / ( 40.0_DP * day_seconds )
    ks     = 1.0_DP / (  4.0_DP * day_seconds )

    allocate( z_kv (1:kmax) )
    z_kv = kf * max( 0.0_DP, &
      &              ( z_Sigma - SigmaB ) / ( 1.0_DP - SigmaB ) )

    allocate( yz_kt (1:jmax, 1:kmax) )

    do k = 1, kmax
      do j = 1, jmax
        yz_kt(j,k) = &
          & ka + ( ks - ka ) &
          &  * max( 0.0_DP, &
          &         ( z_Sigma(k) - SigmaB ) / ( 1.0_DP - SigmaB ) &
          &       ) * cos( y_Lat(j) ) ** 4
      end do
    end do

    ! ヒストリデ〖タ叫蜗のためのへの恃眶判峡
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'DUDtHs94', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'eastward wind tendency', 'm s-2' )
    call HistoryAutoAddVariable( 'DVDtHs94', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'northward wind tendency', 'm s-2' )
    call HistoryAutoAddVariable( 'DTempDtHs94', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'temperature tendency', 'K s-1' )
    call HistoryAutoAddVariable( 'TempEQHs94', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'equilibrium temperature', 'K' )

    ! 磅机 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    held_suarez_1994_inited = .true.
  end subroutine Hs94Init

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

  subroutine Hs94Finalize
    !
    ! モジュ〖ル柒婶の恃眶の充り烧け豺近を乖います.
    !
    ! Deallocate variables in this module.
    !

    ! 离咐矢 ; Declaration statements
    !
    implicit none

    ! 悸乖矢 ; Executable statement
    !

    if ( .not. held_suarez_1994_inited ) return

    ! 充り烧け豺近
    ! Deallocation
    !
    if ( allocated( z_kv  ) ) deallocate( z_kv  )
    if ( allocated( yz_kt ) ) deallocate( yz_kt )

    held_suarez_1994_inited = .false.
  end subroutine Hs94Finalize

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

  subroutine InitCheck
    !
    ! 巴赂モジュ〖ルの介袋步チェック
    !
    ! Check initialization of dependency modules

    ! モジュ〖ル苞脱 ; USE statements
    !

    ! NAMELIST ファイル掐蜗に簇するユ〖ティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_util_inited

    ! 呈灰爬肋年
    ! Grid points settings
    !
    use gridset, only: gridset_inited

    ! 湿妄年眶肋年
    ! Physical constants settings
    !
    use constants, only: constants_inited

    ! 郝筛デ〖タ肋年
    ! Axes data settings
    !
    use axesset, only: axesset_inited

    ! 箕癸瓷妄
    ! Time control
    !
    use timeset, only: timeset_inited

    ! 悸乖矢 ; Executable statement
    !

    if ( .not. namelist_util_inited ) &
      & call MessageNotify( 'E', module_name, '"namelist_util" module is not initialized.' )

    if ( .not. gridset_inited ) &
      & call MessageNotify( 'E', module_name, '"gridset" module is not initialized.' )

    if ( .not. constants_inited ) &
      & call MessageNotify( 'E', module_name, '"constants" module is not initialized.' )

    if ( .not. axesset_inited ) &
      & call MessageNotify( 'E', module_name, '"axesset" module is not initialized.' )

    if ( .not. timeset_inited ) &
      & call MessageNotify( 'E', module_name, '"timeset" module is not initialized.' )

  end subroutine InitCheck

end module held_suarez_1994