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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, GasRDry
                              ! $ R $ [J kg-1 K-1]. 
                              ! 乾燥大気の気体定数. 
                              ! Gas constant of air

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, 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_Amp       (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 :: mmax
    integer :: ms
    integer :: me
    real(DP) :: az_A(0:imax*jmax-1, 1:kmax)
    real(DP) :: az_B(0:imax*jmax-1, 1:kmax)
    real(DP) :: az_C(0:imax*jmax-1, 1:kmax)
    real(DP) :: az_D(0:imax*jmax-1, 1:kmax)



    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 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

    integer:: k


    ! 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, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = gwd_m1987_nml, 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

  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