Class gwd_m1987
In: gwdrag/gwd_m1987.f90

Gravity wave drag by McFarlane (1987)

Gravity wave drag by McFarlane (1987)

Note that Japanese and English are described in parallel.

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

Methods

Included Modules

gridset dc_types namelist_util dc_message constants0 constants timeset gtool_historyauto dc_calendar dc_iounit dc_string axesset

Public Instance methods

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

Tendency of gravity wave drag based on McFarlane (1987)

[Source]

  subroutine GWDM1987( xyz_U, xyz_V, xyz_Temp, xyz_Press, xyr_Press, xyz_Exner, xyz_Height, 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_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_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 :

moist_conv_adjust モジュールの初期化を行います. NAMELIST#moist_conv_adjust_nml の読み込みはこの手続きで行われます.

"moist_conv_adjust" module is initialized. "NAMELIST#moist_conv_adjust_nml" is loaded in this procedure.

This procedure input/output NAMELIST#gwd_m1987_nml .

[Source]

  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

Private Instance methods

CrtlFNumSq
Variable :
CrtlFNumSq :real(DP), save
: Critical Froude number squared
Efficiency
Variable :
Efficiency :real(DP), save
: "Efficiency"
FlagDetermineRefLevByStd
Variable :
FlagDetermineRefLevByStd :logical , save
MomFluxFactor
Variable :
MomFluxFactor :real(DP), save
: Factor for momentum flux
OrogEffWaveLength
Variable :
OrogEffWaveLength :real(DP), save
: Orography effective wave length
SigmaRef
Variable :
SigmaRef :real(DP), save
: Sigma at reference level
gwd_m1987_inited
Variable :
gwd_m1987_inited = .false. :logical, save
: 初期設定フラグ. Initialization flag
module_name
Constant :
module_name = ‘gwd_m1987 :character(*), parameter
: モジュールの名称. Module name
Subroutine :
mm :integer , intent(in )
jmx :integer , intent(in )
a( mm, jmx ) :real(DP), intent(in )
b( mm, jmx ) :real(DP), intent(in )
c( mm, jmx ) :real(DP), intent(in )
f( mm, jmx ) :real(DP), intent(inout)
ms :integer , intent(in )
me :integer , intent(in )

[Source]

  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
version
Constant :
version = ’$Name: dcpam5-20150131 $’ // ’$Id: gwd_m1987.f90,v 1.1 2015/01/29 12:22:29 yot Exp $’ :character(*), parameter
: モジュールのバージョン Module version