!= Gravity wave drag by McFarlane (1987)
!
!= Gravity wave drag by McFarlane (1987)
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: gwd_m1987.f90,v 1.2 2015/03/11 04:54:29 yot Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module gwd_m1987
  !
  != Gravity wave drag by McFarlane (1987)
  !
  != Gravity wave drag by McFarlane (1987)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! Calculate tendency by gravity wave drag
  !
  !== References
  !
  !  McFarlane, N. A., 
  !    The effect of orographically excited gravity wave drag on the general
  !    circulation of the lower stratosphere and troposphsere,
  !    J. Atmos. Sci., 44, 1775-1800, 1987.
  !
  !== Procedures List
  !
  ! GWDM1987     :: Calculation of gravity wave drag tendency
  ! GWDM1987Init :: Initialization
  ! ------------ :: ------------
  ! GWDM1987     :: Calculation of gravity wave drag tendency
  ! GWDM1987Init :: Initialization
  !
  !== NAMELIST
  !
  ! NAMELIST#gwd_m1987_nml
  !

  ! モジュール引用 ; 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. 

  ! NAMELIST ファイル入力に関するユーティリティ
  ! Utilities for NAMELIST file input
  !
  use namelist_util, only: MaxNmlArySize
                              ! NAMELIST から読み込む配列の最大サイズ. 
                              ! Maximum size of arrays loaded from NAMELIST

  ! メッセージ出力
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 物理・数学定数設定
  ! Physical and mathematical constants settings
  !
  use constants0, only: &
    & PI                    ! $ \pi $.
                            ! 円周率. Circular constant


  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

!!$  logical , save :: FlagUse
!!$                              ! 使用フラグ
!!$                              ! flag for use of this scheme

  logical , save :: FlagDetermineRefLevByStd
  real(DP), save :: SigmaRef
                              ! Sigma at reference level
  real(DP), save :: Efficiency
                              ! "Efficiency"
  real(DP), save :: OrogEffWaveLength
                              ! Orography effective wave length
  real(DP), save :: CrtlFNumSq
                              ! Critical Froude number squared
  real(DP), save :: MomFluxFactor
                              ! Factor for momentum flux

!!$  logical, save  :: FlagGWDDamp
!!$  real(DP), save :: GWDDampPeriod

  ! 公開手続き
  ! Public procedure
  !
  public :: GWDM1987
  public :: GWDM1987Init

  ! 公開変数
  ! Public variables
  !

  ! 非公開変数
  ! Private variables
  !
  logical, save :: gwd_m1987_inited = .false.
                              ! 初期設定フラグ. 
                              ! Initialization flag

  character(*), parameter:: module_name = 'gwd_m1987'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: gwd_m1987.f90,v 1.2 2015/03/11 04:54:29 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

  subroutine GWDM1987(                                 &
    & xyz_U, xyz_V, xyz_Temp,                          & ! (in)
    & xyz_Press, xyr_Press, xyz_Exner, xyz_Height,     & ! (in)
    & xy_SurfHeight, xy_SurfHeightStd,                 & ! (in)
    & xyz_DUDt, xyz_DVDt                               & ! (out)
    & )
    !
    ! 
    !
    ! Tendency of gravity wave drag based on McFarlane (1987)
    !

    ! モジュール引用 ; USE statements
    !

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Grav, & 
                              ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration
      & GasRDry
                              ! $ R $ [J kg-1 K-1]. 
                              ! 乾燥大気の気体定数. 
                              ! Gas constant of air

    ! 時刻管理
    ! 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
    !

    real(DP), intent(in ) :: xyz_U       (0:imax-1, 1:jmax, 1:kmax)
                              ! Zonal wind
    real(DP), intent(in ) :: xyz_V       (0:imax-1, 1:jmax, 1:kmax)
                              ! Meridional wind
    real(DP), intent(in ) :: xyz_Temp    (0:imax-1, 1:jmax, 1:kmax)
                              ! Temperature
    real(DP), intent(in ) :: xyz_Press   (0:imax-1, 1:jmax, 1:kmax)
                              ! Pressure
    real(DP), intent(in ) :: xyr_Press   (0:imax-1, 1:jmax, 0:kmax)
                              ! Pressure
    real(DP), intent(in ) :: xyz_Exner   (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner function
    real(DP), intent(in ) :: xyz_Height   (0:imax-1, 1:jmax, 1:kmax)
                              ! Height
    real(DP), intent(in ) :: xy_SurfHeight   (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xy_SurfHeightStd(0:imax-1, 1:jmax)
    real(DP), intent(out) :: xyz_DUDt    (0:imax-1, 1:jmax, 1:kmax)
                              ! 東西風変化率. 
                              ! Zonal wind tendency
    real(DP), intent(out) :: xyz_DVDt    (0:imax-1, 1:jmax, 1:kmax)
                              ! 南北風変化率. 
                              ! Meridional wind tendency


    ! 作業変数
    ! Work variables
    !
    real(DP) :: xy_OrogEffHeight(0:imax-1, 1:jmax)

    real(DP) :: xyz_DelPress   (0:imax-1, 1:jmax, 1:kmax)

    integer  :: xy_KIndexRef   (0:imax-1, 1:jmax)
    real(DP) :: xy_URef        (0:imax-1, 1:jmax)
    real(DP) :: xy_VRef        (0:imax-1, 1:jmax)
    real(DP) :: xyz_WindSpeed  (0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xy_AbsWindSpeedRef(0:imax-1, 1:jmax)
    real(DP) :: xy_RhoRef         (0:imax-1, 1:jmax)
    real(DP) :: xy_NRef           (0:imax-1, 1:jmax)

    real(DP) :: xyz_Rho       (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_PotTemp   (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_N         (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_ZeroOne   (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xy_MomFluxRef (0:imax-1, 1:jmax)
    real(DP) :: xyz_MomFlux   (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: MomFluxSaturated
    real(DP) :: xyz_DMomFluxDU(0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xyz_DWindSpeedDt(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyr_MomFluxA (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_MomFluxXA(0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_MomFluxYA(0:imax-1, 1:jmax, 0:kmax)

    real(DP) :: WindSpeedTentative

!!$    real(DP) :: MomFluxTentative


    integer :: i               ! 経度方向に回る DO ループ用作業変数
                               ! Work variables for DO loop in longitude
    integer :: j               ! 緯度方向に回る DO ループ用作業変数
                               ! Work variables for DO loop in latitude
    integer :: k               ! 鉛直方向に回る DO ループ用作業変数
                               ! Work variables for DO loop in vertical direction
    integer :: kp
    integer :: kn

!!$    real(DP) :: xyz_WindSpeedTentative(0:imax-1, 1:jmax, 1:kmax)
!!$    integer :: itr


    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. gwd_m1987_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if

    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )


    ! Calculation of additional variables
    !
    xy_OrogEffHeight = 2.0_DP * xy_SurfHeightStd
    !
    do k = 1, kmax
      xyz_DelPress(:,:,k) = xyr_Press(:,:,k-1) - xyr_Press(:,:,k)
    end do
    !
    !   Determine reference level
    !
    if ( FlagDetermineRefLevByStd ) then
      xy_KIndexRef = 2
      do k = 1+1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( ( xyz_Height(i,j,k) - xy_SurfHeight(i,j) ) < xy_OrogEffHeight(i,j) ) then
              xy_KIndexRef(i,j) = k
            end if
          end do
        end do
      end do
    else
      xy_KIndexRef = 2
      do k = 1+1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( xyz_Press(i,j,k) / xyr_Press(i,j,0) > SigmaRef ) then
              xy_KIndexRef(i,j) = k
            end if
          end do
        end do
      end do
    end if

    !
    !   Set reference level wind velocity
    !
    do j = 1, jmax
      do i = 0, imax-1
        xy_URef = xyz_U(i,j,xy_KIndexRef(i,j))
        xy_VRef = xyz_V(i,j,xy_KIndexRef(i,j))
      end do
    end do
    xy_AbsWindSpeedRef = sqrt( xy_URef**2 + xy_VRef**2 )

    ! Calculation of additional variables
    !
    xyz_Rho = xyz_Press / ( GasRDry * xyz_Temp )
    do j = 1, jmax
      do i = 0, imax-1
        xy_RhoRef = xyz_Rho(i,j,xy_KIndexRef(i,j))
      end do
    end do
    !
    xyz_PotTemp = xyz_Temp / xyz_Exner
    !
    do k = 1, kmax
      kp = max( k - 1, 1    )
      kn = min( k + 1, kmax )
      xyz_N(:,:,k) =                                      &
        & Grav / xyz_PotTemp(:,:,k)                       &
        & * ( xyz_PotTemp(:,:,kn) - xyz_PotTemp(:,:,kp) ) &
        & / ( xyz_Height (:,:,kn) - xyz_Height (:,:,kp) )
    end do
!!$    xyz_N = max( xyz_N, 1.0e-6_DP )
    xyz_N = max( xyz_N, 0.0_DP )
    xyz_N = sqrt( xyz_N )
    do j = 1, jmax
      do i = 0, imax-1
        xy_NRef = xyz_N(i,j,xy_KIndexRef(i,j))
      end do
    end do
    !
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_AbsWindSpeedRef(i,j) /= 0.0_DP ) then
            xyz_WindSpeed(i,j,k) = &
              & ( xyz_U(i,j,k) * xy_URef(i,j) + xyz_V(i,j,k) * xy_VRef(i,j) ) &
              & / xy_AbsWindSpeedRef(i,j)
          else
            xyz_WindSpeed(i,j,k) = 0.0_DP
          end if
        end do
      end do
    end do
    ! Negative wind speed is inrelevant in the current formulation.
    xyz_WindSpeed = max( xyz_WindSpeed, 0.0_DP )


    ! Wave amplitude
    !
    ! Momentum flux parallel to the reference level flow at reference level
    !
    xy_MomFluxRef = - MomFluxFactor * xy_OrogEffHeight**2 &
      & * xy_RhoRef * xy_NRef * xy_AbsWindSpeedRef

    !
    ! Momentum flux parallel to the reference level flow
    !
    do j = 1, jmax
      do i = 0, imax-1
        ! Region at and below the reference level
        do k = 1, xy_KIndexRef(i,j)
          xyz_MomFlux(i,j,k) = xy_MomFluxRef(i,j)
          xyz_ZeroOne(i,j,k) = 0.0_DP
        end do
        ! Region above the reference level and below the highest model level
        do k = xy_KIndexRef(i,j)+1, kmax-1
          ! momentum flux is same as that in lower level tentatively
          xyz_MomFlux(i,j,k) = xyz_MomFlux(i,j,k-1)
          ! calculate momentum flux in the case of saturation
          if ( xyz_N(i,j,k) > 0.0_DP ) then
            MomFluxSaturated = &
              & - MomFluxFactor * CrtlFNumSq &
              &   * xyz_Rho(i,j,k) / xyz_N(i,j,k) * xyz_WindSpeed(i,j,k)**3
          else
            MomFluxSaturated = 0.0_DP
          end if
          ! comparison of momentum flux
          ! check whether saturationed or not
          ! It should be noted that momentum flux here is negative, since 
          ! a direction of the momentum flux is parallel to reference level 
          ! flow.
          if ( MomFluxSaturated > xyz_MomFlux(i,j,k) ) then
            ! saturation region
            xyz_MomFlux(i,j,k) = MomFluxSaturated
            xyz_ZeroOne(i,j,k) = 1.0_DP
          else
            ! non-saturation region
            xyz_ZeroOne(i,j,k) = 0.0_DP
          end if
        end do
        ! highest model level
        do k = kmax, kmax
          ! momentum flux is same as that in lower level
          xyz_MomFlux(i,j,k) = xyz_MomFlux(i,j,k-1)
          xyz_ZeroOne(i,j,k) = 0.0_DP
        end do
      end do
    end do
    !
    ! Momentum flux derivative with reference level flow
    !
!!$    xyz_DMomFluxDU =                                            &
!!$      & - 3.0_DP * MomFluxFactor * CrtlFNumSq * xyz_Rho / xyz_N &
!!$      &   * xyz_WindSpeed**2                                    &
!!$      &   * xyz_ZeroOne
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_N(i,j,k) > 0.0_DP ) then
            xyz_DMomFluxDU(i,j,k) =                   &
              & - 3.0_DP * MomFluxFactor * CrtlFNumSq &
              &   * xyz_Rho(i,j,k) / xyz_N(i,j,k)     &
              &   * xyz_WindSpeed(i,j,k)**2           &
              &   * xyz_ZeroOne(i,j,k)
          else
            xyz_DMomFluxDU(i,j,k) = 0.0_DP
          end if
        end do
      end do
    end do


    !
    ! calculation of wind velocity tendency
    !
    do j = 1, jmax
      do i = 0, imax-1
        ! No deceleration is assumed in the highest level
        k = kmax
        xyz_DWindSpeedDt(i,j,k) = 0.0_DP
        do k = kmax-1, 1+1, -1
          xyz_DWindSpeedDt(i,j,k) = &
            & Grav / xyz_DelPress(i,j,k) / 2.0_DP                         &
            & * (   xyz_MomFlux(i,j,k-1)                                  &
            &     - xyz_MomFlux(i,j,k+1)                                  &
            &     - xyz_DMomFluxDU(i,j,k+1) * xyz_DWindSpeedDt(i,j,k+1)   &
            &       * ( 2.0_DP * DelTime )                              ) &
            & / ( 1.0_DP - Grav / xyz_DelPress(i,j,k) / 2.0_DP &
            &              * xyz_DMomFluxDU(i,j,k) * ( 2.0_DP * DelTime ) )

          ! Wind speed tendency at k level and momentum flux at k-1 level 
          ! are estimated again.
          if ( k >= xy_KIndexRef(i,j) ) then
            ! Region above reference level
            WindSpeedTentative =       &
              &   xyz_WindSpeed(i,j,k) &
              & + xyz_DWindSpeedDt(i,j,k) * ( 2.0_DP * DelTime )
            if ( WindSpeedTentative < 0.0_DP ) then
              xyz_DWindSpeedDt(i,j,k) = &
                & ( 0.0_DP - xyz_WindSpeed(i,j,k) ) / ( 2.0_DP * DelTime )
              xyz_MomFlux(i,j,k-1) = &
                &   (   2.0_DP * xyz_DelPress(i,j,k) / Grav                   &
                &     - xyz_DMomFluxDU(i,j,k) * ( 2.0_DP * DelTime ) )        &
                &   * xyz_DWindSpeedDt(i,j,k)                                 &
                & + xyz_MomFlux(i,j,k+1)                                      &
                & + xyz_DMomFluxDU(i,j,k+1) * xyz_DWindSpeedDt(i,j,k+1)       &
                &       * ( 2.0_DP * DelTime )
            end if
          else
            ! below the reference level
            xyz_DWindSpeedDt(i,j,k) = 0.0_DP
            xyz_MomFlux(i,j,k-1) = &
              &   (   2.0_DP * xyz_DelPress(i,j,k) / Grav                   &
              &     - xyz_DMomFluxDU(i,j,k) * ( 2.0_DP * DelTime ) )        &
              &   * xyz_DWindSpeedDt(i,j,k)                                 &
              & + xyz_MomFlux(i,j,k+1)                                      &
              & + xyz_DMomFluxDU(i,j,k+1) * xyz_DWindSpeedDt(i,j,k+1)       &
              &       * ( 2.0_DP * DelTime )
          end if

        end do
        ! No deceleration is assumed in the lowest level
        k = 1
        xyz_DWindSpeedDt(i,j,k) = 0.0_DP
      end do
    end do



!!$    ! No deceleration is assumed in the lowest level
!!$    k = 1
!!$    xyz_DWindSpeedDt(:,:,k) = 0.0_DP
!!$    do k = 1+1, kmax-1
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          !
!!$          ! calculation of wind velocity tendency
!!$          !
!!$          xyz_DWindSpeedDt(i,j,k) = &
!!$            & Grav / xyz_DelPress(i,j,k) / 2.0_DP &
!!$            & * (   xyz_MomFlux(i,j,k-1) &
!!$            &     + xyz_DMomFluxDU(i,j,k-1) * xyz_DWindSpeedDt(i,j,k-1) &
!!$            &       * ( 2.0_DP * DelTime ) &
!!$            &     - xyz_MomFlux(i,j,k+1) ) &
!!$            & / ( 1.0_DP + Grav / xyz_DelPress(i,j,k) / 2.0_DP &
!!$            &              * xyz_DMomFluxDU(i,j,k) * ( 2.0_DP * DelTime ) )



!!$          xyz_DWindSpeedDt(i,j,k) = 0.0_DP
!!$
!!$          do itr = 1, 10
!!$
!!$            xyz_WindSpeedTentative(i,j,k) = &
!!$              & xyz_WindSpeed(i,j,k) + xyz_DWindSpeedDt(i,j,k) * ( 2.0_DP * DelTime )
!!$
!!$            if ( xyz_N(i,j,k) > 0.0_DP ) then
!!$              xyz_DMomFluxDU(i,j,k) =                   &
!!$                & - 3.0_DP * MomFluxFactor * CrtlFNumSq &
!!$                &   * xyz_Rho(i,j,k) / xyz_N(i,j,k)     &
!!$                &   * xyz_WindSpeedTentative(i,j,k)**2           &
!!$                &   * xyz_ZeroOne(i,j,k)
!!$            else
!!$              xyz_DMomFluxDU(i,j,k) = 0.0_DP
!!$            end if
!!$
!!$            xyz_DWindSpeedDt(i,j,k) = &
!!$              & Grav / xyz_DelPress(i,j,k) / 2.0_DP &
!!$              & * (   xyz_MomFlux(i,j,k-1) &
!!$              &     + xyz_DMomFluxDU(i,j,k-1) * xyz_DWindSpeedDt(i,j,k-1) &
!!$              &       * ( 2.0_DP * DelTime ) &
!!$              &     - xyz_MomFlux(i,j,k+1) ) &
!!$              & / ( 1.0_DP + Grav / xyz_DelPress(i,j,k) / 2.0_DP &
!!$              &              * xyz_DMomFluxDU(i,j,k) * ( 2.0_DP * DelTime ) )
!!$
!!$          end do


!!$          !
!!$          ! preparation for next level
!!$          !
!!$          !   estimated momentum flux at k and next time step
!!$          MomFluxTentative =                                    &
!!$            &   xyz_MomFlux(i,j,k)                              &
!!$            & + xyz_DMomFluxDU(i,j,k) * xyz_DWindSpeedDt(i,j,k) &
!!$            &   * ( 2.0_DP * DelTime )
!!$          xyz_MomFlux(i,j,k+1) = MomFluxTentative
!!$          !   calculate momentum flux at k+1 in the case of saturation
!!$          MomFluxSaturated = &
!!$            & - MomFluxFactor * CrtlFNumSq &
!!$            &   * xyz_Rho(i,j,k+1) / xyz_N(i,j,k+1) * xyz_WindSpeed(i,j,k+1)**3
!!$          !   comparison of momentum flux
!!$          !   check whether saturationed or not
!!$          if ( abs( MomFluxSaturated ) < abs( xyz_MomFlux(i,j,k+1) ) ) then
!!$            ! saturation region
!!$            !   set saturated momentum flux
!!$            xyz_MomFlux(i,j,k+1) = MomFluxSaturated
!!$            !   derivative of momentum flux with reference level flow
!!$            xyz_DMomFluxDU(i,j,k+1) =                 &
!!$              & - 3.0_DP * MomFluxFactor * CrtlFNumSq &
!!$              &   * xyz_Rho(i,j,k+1) / xyz_N(i,j,k+1) &
!!$              &   * xyz_WindSpeed(i,j,k+1)**2
!!$          else
!!$            ! non-saturation region
!!$            !   derivative of momentum flux with reference level flow
!!$            xyz_DMomFluxDU(i,j,k+1) = 0.0_DP
!!$          end if
!!$        end do
!!$      end do
!!$    end do
!!$    ! No deceleration is assumed in the highest level
!!$    k = kmax
!!$    xyz_DWindSpeedDt(:,:,k) = 0.0_DP

    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_AbsWindSpeedRef(i,j) /= 0.0_DP ) then
            xyz_DUDt(i,j,k) = &
              & xyz_DWindSpeedDt(i,j,k) * xy_URef(i,j) / xy_AbsWindSpeedRef(i,j)
            xyz_DVDt(i,j,k) = &
              & xyz_DWindSpeedDt(i,j,k) * xy_VRef(i,j) / xy_AbsWindSpeedRef(i,j)
          else
            xyz_DUDt(i,j,k) = 0.0_DP
            xyz_DVDt(i,j,k) = 0.0_DP
          end if
        end do
      end do
    end do

!!$    if ( FlagGWDDamp ) then
!!$
!!$      GWDDampCoef = 1.0_DP / DelTime * ( GWDDampPeriod - TimeN ) / GWDDampPeriod
!!$
!!$      if ( GWDDampCoef < 0.0_DP ) GWDDampCoef = 0.0_DP
!!$
!!$      xyz_DUDt = xyz_DUDt / ( 1.0_DP + 2.0_DP * DelTime * DivDampCoef )
!!$
!!$    end if


    ! estimated momentum flux at layer interface and at next time step
    do k = 0+1, kmax-1
      xyr_MomFluxA(:,:,k) =                                           &
        &   (   xyz_MomFlux(:,:,k)                                    &
        &     + xyz_MomFlux(:,:,k+1)                                  &
        &       + xyz_DMomFluxDU(:,:,k+1) * xyz_DWindSpeedDt(:,:,k+1) &
        &         * ( 2.0_DP * DelTime ) )                            &
        & / 2.0_DP
!!$      xyr_MomFluxA(:,:,k) =                                         &
!!$        &   (   xyz_MomFlux(:,:,k)                                  &
!!$        &       + xyz_DMomFluxDU(:,:,k) * xyz_DWindSpeedDt(:,:,k)   &
!!$        &         * ( 2.0_DP * DelTime )                            &
!!$        &     + xyz_MomFlux(:,:,k+1)                              ) &
!!$        & / 2.0_DP
    end do
    k = 0
    xyr_MomFluxA(:,:,k) = xyr_MomFluxA(:,:,k+1)
    k = kmax
    xyr_MomFluxA(:,:,k) = xyr_MomFluxA(:,:,k-1)


    do k = 0, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_AbsWindSpeedRef(i,j) /= 0.0_DP ) then
            xyr_MomFluxXA(i,j,k) = &
              & xyr_MomFluxA(i,j,k) * xy_URef(i,j) / xy_AbsWindSpeedRef(i,j)
            xyr_MomFluxYA(i,j,k) = &
              & xyr_MomFluxA(i,j,k) * xy_VRef(i,j) / xy_AbsWindSpeedRef(i,j)
          else
            xyr_MomFluxXA(i,j,k) = 0.0_DP
            xyr_MomFluxYA(i,j,k) = 0.0_DP
          end if
        end do
      end do
    end do


!!$    ! Set up of simultaneously linear equations
!!$    do k = 1, kmax
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          az_A(i+imax*(j-1),k) =                    &
!!$            &   Grav / xyz_DelPress(i,j,k) / 2.0_DP &
!!$            &   * xyz_DMomFluxDU(i,j,k-1)           &
!!$            &   * ( 2.0_DP * DelTime )
!!$          az_A(i+imax*(j-1),k) = - 1.0_DP
!!$          az_C(i+imax*(j-1),k) =                    &
!!$            & - Grav / xyz_DelPress(i,j,k) / 2.0_DP &
!!$            &   * xyz_DMomFluxDU(i,j,k+1)           &
!!$            &   * ( 2.0_DP * DelTime )
!!$          az_D(i+imax*(j-1),k) =                    &
!!$            & - Grav / xyz_DelPress(i,j,k) / 2.0_DP &
!!$            &  * ( xyz_MomFlux(i,j,k-1) - xyz_MomFlux(i,j,k+1) )
!!$        end do
!!$      end do
!!$    end do
!!$
!!$    mmax = imax*jmax
!!$    ms = 1
!!$    me = imax*jmax
!!$    call tridiag( mmax, kmax, az_A, az_B, az_C, az_D, ms, me )
!!$
!!$    do k = 1, kmax
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          xyz_DUDt = az_D(i+imax*(j-1),k) * xy_URef(i,j) / xy_WindSpeedRef(i,j)
!!$          xyz_DVDt = az_D(i+imax*(j-1),k) * xy_VRef(i,j) / xy_WindSpeedRef(i,j)
!!$        end do
!!$      end do
!!$    end do


    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'GWMomFlux' , xyr_MomFluxA  )
    call HistoryAutoPut( TimeN, 'GWMomFluxX', xyr_MomFluxXA )
    call HistoryAutoPut( TimeN, 'GWMomFluxY', xyr_MomFluxYA )
    call HistoryAutoPut( TimeN, 'DUDtGWD'   , xyz_DUDt      )
    call HistoryAutoPut( TimeN, 'DVDtGWD'   , xyz_DVDt      )
    call HistoryAutoPut( TimeN, 'DWSDtGWD'  , xyz_DWindSpeedDt )


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine GWDM1987

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

  !******************************************************************************
  !      subroutine tridiag
  !      tidiagonal solver
  !******************************************************************************
  !     a(j), b(j), and c(j) are, respectively, the subdiagonal, diagonal,
  !     and superdiagonal entries in row j.
  !     a(1) and c(jmx) need not be initialized.
  !     The output is in f; a, b, and c are unchanged.
  !******************************************************************************
  !     jmx    : dimensions of all the following arrays
  !     a      : sub (lower) diagonal
  !     b      : center diagonal
  !     c      : super (upper) diagonal
  !     f      : right hand side
  !******************************************************************************

  subroutine tridiag( mm, jmx, a, b, c, f, ms, me )

    integer , intent(in   ) :: mm, jmx
    real(DP), intent(in   ) :: a( mm, jmx ),b( mm, jmx )
    real(DP), intent(in   ) :: c( mm, jmx )
    real(DP), intent(inout) :: f( mm, jmx )
    integer , intent(in   ) :: ms, me


    ! Local variables
    !
    real(DP) :: q( mm, jmx ), p
    integer  :: j, m


    ! Forward elimination sweep
    !
    do m = ms, me
      q( m, 1 ) = - c( m, 1 ) / b( m, 1 )
      f( m, 1 ) =   f( m, 1 ) / b( m, 1 )
    end do

    do j = 2, jmx
      do m = ms, me
        p         = 1.0d0 / ( b( m, j ) + a( m, j ) * q( m, j-1 ) )
        q( m, j ) = - c( m, j ) * p
        f( m, j ) = ( f( m, j ) - a( m, j ) * f( m, j-1 ) ) * p
      end do
    end do

    ! Backward pass
    !
    do j = jmx - 1, 1, -1
      do m = ms, me
        f( m, j ) = f( m, j ) + q( m, j ) * f( m, j+1 )
      end do
    end do

  end subroutine tridiag

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

  subroutine GWDM1987Init
    !
    ! moist_conv_adjust モジュールの初期化を行います. 
    ! NAMELIST#moist_conv_adjust_nml の読み込みはこの手続きで行われます. 
    !
    ! "moist_conv_adjust" module is initialized. 
    ! "NAMELIST#moist_conv_adjust_nml" is loaded in this procedure. 
    !

    ! モジュール引用 ; USE statements
    !

    ! 日付および時刻の取り扱い
    ! Date and time handler
    !
    use dc_calendar, only: DCCalConvertByUnit

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! 文字列操作
    ! Character handling
    !
    use dc_string, only: StoA

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoAddVariable

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & AxnameX, &
      & AxnameY, &
      & AxnameZ, &
      & AxnameR, &
      & AxnameT

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only : &
      RPlanet
                              ! $ a $ [m].
                              ! 惑星半径.
                              ! Radius of planet


    ! 宣言文 ; Declaration statements
    !
    implicit none

    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                              ! IOSTAT of NAMELIST read

!!$    real(DP)         :: GWDDampPeriodValue
!!$    character(TOKEN) :: GWDDampPeriodUnit


    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /gwd_m1987_nml/ &
      & FlagDetermineRefLevByStd, &
      & SigmaRef, &
      & CrtlFNumSq, &
      & Efficiency, &
      & OrogEffWaveLength !, &
!!$      & FlagGWDDamp, GWDDampPeriodValue, GWDDampPeriodUnit
          ! デフォルト値については初期化手続 "moist_conv_adjust#CumAdjInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "moist_conv_adjust#MoistConvAdjustInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( gwd_m1987_inited ) return


    ! デフォルト値の設定
    ! Default values settings
    !
    FlagDetermineRefLevByStd = .false.

    SigmaRef           = 1.0_DP   ! lowest model level

    CrtlFNumSq         = 0.5_DP   ! value used by McFarlane (1987)
    Efficiency         = 1.0_DP   ! arbitrary
!!$    OrogEffWaveLength  = 2.0_DP * PI / ( 2.0_DP * 8.0e-6_DP / Efficiency )
                                  ! value used by McFarlane (1987)
    if ( imax /= 1 ) then
      OrogEffWaveLength  = 2.0_DP * PI * RPlanet / ( ( imax - 1 ) / 3.0_DP )
                                  ! wavelength of smallest resolved wave
    else
      OrogEffWaveLength  = 2.0_DP * PI / ( 2.0_DP * 8.0e-6_DP / Efficiency )
                                  ! value used by McFarlane (1987)
    end if

!!$    FlagGWDDamp        = .false.
!!$    GWDDampPeriodValue = 0.0_DP
!!$    GWDDampPeriodUnit  = 'day'


    ! 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 = gwd_m1987_nml,            &  ! (out)
        & iostat = iostat_nml )              ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if


    ! Calculate factor for momentum flux
    MomFluxFactor = Efficiency * 2.0_DP * PI / OrogEffWaveLength / 2.0_DP

    ! Check values
    !
    SigmaRef = max( min( SigmaRef, 1.0_DP ), 0.0_DP )

    ! Calculation of divergence damping period
    !
!!$    GWDDampPeriod = DCCalConvertByUnit( GWDDampPeriodValue, GWDDampPeriodUnit, 'sec' )

    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'GWMomFlux', &
      & (/ AxNameX, AxNameY, AxNameR, AxNameT /), &
      & 'gravity wave momentum flux', 'kg m-1 s-2' )
    call HistoryAutoAddVariable( 'GWMomFluxX', &
      & (/ AxNameX, AxNameY, AxNameR, AxNameT /), &
      & 'gravity wave zonal momentum flux', 'kg m-1 s-2' )
    call HistoryAutoAddVariable( 'GWMomFluxY', &
      & (/ AxNameX, AxNameY, AxNameR, AxNameT /), &
      & 'gravity wave meridional momentum flux', 'kg m-1 s-2' )
    call HistoryAutoAddVariable( 'DUDtGWD', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'zonal wind tendency by gravity wave drag', 'm-1 s-2' )
    call HistoryAutoAddVariable( 'DVDtGWD', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'meridional wind tendency by gravity wave drag', 'm-1 s-2' )
    call HistoryAutoAddVariable( 'DWSDtGWD', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'wind speed tendency by gravity wave drag', 'm-1 s-2' )


    ! Initialization of modules used in this module
    !


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '  FlagDetermineRefLevByStd = %b', l = (/ FlagDetermineRefLevByStd /) )
    call MessageNotify( 'M', module_name, '  SigmaRef            = %f', d = (/ SigmaRef /) )
    call MessageNotify( 'M', module_name, '  CrtlFNumSq          = %f', d = (/ CrtlFNumSq /) )
    call MessageNotify( 'M', module_name, '  Efficiency          = %f', d = (/ Efficiency /) )
    call MessageNotify( 'M', module_name, '  OrogEffWaveLength   = %f', d = (/ OrogEffWaveLength /) )
    call MessageNotify( 'M', module_name, '    MomFluxFactor     = %f', d = (/ MomFluxFactor /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    gwd_m1987_inited = .true.

  end subroutine GWDM1987Init

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

end module gwd_m1987