Class sltt_dp
In: sltt/sltt_dp.f90

上流点探索

Finding Departure Point

Note that Japanese and English are described in parallel.

Procedures List

SLTTDPHor :水平方向の上流点探索
SLTTDPVer :鉛直方向の上流点探索
——————— :————
SLTTDPHor :Finding DP in Horizontal
SLTTDPVer :Finding DP in Vertical

NAMELIST

NAMELIST#

References

  • Williamson, D. L., and Rasch, P. J., 1989: Two-dimensional semi-Lagrangian transport with shape-preserving interpolation. Mon. Wea. Rev., 117, 102-129.

Methods

Included Modules

dc_types gridset sltt_const constants0 constants mpi_wrapper sltt_lagint axesset

Public Instance methods

Subroutine :
DelTime :real(DP), intent(in )
: $Delta t$
x_Lon(0:imax-1) :real(DP), intent(in )
: $lambda$ longitude
y_Lat(1:jmax/2) :real(DP), intent(in )
: $varphi$ latitude
y_SinLat(1:jmax/2) :real(DP), intent(in )
: $sin\varphi$
y_CosLat(1:jmax/2) :real(DP), intent(in )
: $cos\varphi$
iexmin :integer , intent(in )
iexmax :integer , intent(in )
jexmin :integer , intent(in )
jexmax :integer , intent(in )
x_ExtLon(iexmin:iexmax) :real(DP), intent(in )
: 経度の拡張配列 Extended array of Lon
y_ExtLat(jexmin:jexmax) :real(DP), intent(in )
: 緯度の拡張配列 Extended array of Lat
xyz_ExtU(iexmin:iexmax, jexmin:jexmax, 1:kmax) :real(DP), intent(in )
: 東西風速の拡張配列 Extended array of zonal velocity
xyz_ExtV(iexmin:iexmax, jexmin:jexmax, 1:kmax) :real(DP), intent(in )
: 南北風速の拡張配列 Extended array of meridional velocity
xyz_DPLon(0:imax-1, 1:jmax/2, 1:kmax) :real(DP), intent(out)
: 上流点の経度 Lon of departure point
xyz_DPLat(0:imax-1, 1:jmax/2, 1:kmax) :real(DP), intent(out)
: 上流点の緯度 Lat of departure point

水平方向の上流点探索 Finding DP in Horizontal

[Source]

  subroutine SLTTDPHor( DelTime, x_Lon, y_Lat, y_SinLat, y_CosLat, iexmin, iexmax, jexmin, jexmax, x_ExtLon, y_ExtLat, xyz_ExtU, xyz_ExtV, xyz_DPLon, xyz_DPLat )
    ! 水平方向の上流点探索
    ! Finding DP in Horizontal

    use sltt_const, only : PIx2, dtjw, nloop_dp_h


    real(DP), intent(in ) :: DelTime
                              ! $\Delta t$         
    real(DP), intent(in ) :: x_Lon(0:imax-1)
                              ! $\lambda$ longitude 
    real(DP), intent(in ) :: y_Lat(1:jmax/2)
                              ! $\varphi$ latitude
    real(DP), intent(in ) :: y_SinLat(1:jmax/2)
                              ! $\sin\varphi$
    real(DP), intent(in ) :: y_CosLat(1:jmax/2)
                              ! $\cos\varphi$
    integer , intent(in ) :: iexmin
    integer , intent(in ) :: iexmax
    integer , intent(in ) :: jexmin
    integer , intent(in ) :: jexmax
    real(DP), intent(in ) :: x_ExtLon(iexmin:iexmax)
                              ! 経度の拡張配列
                              ! Extended array of Lon
    real(DP), intent(in ) :: y_ExtLat(jexmin:jexmax)
                              ! 緯度の拡張配列
                              ! Extended array of Lat    
    real(DP), intent(in ) :: xyz_ExtU(iexmin:iexmax, jexmin:jexmax, 1:kmax)    
                              ! 東西風速の拡張配列
                              ! Extended array of zonal velocity
    real(DP), intent(in ) :: xyz_ExtV(iexmin:iexmax, jexmin:jexmax, 1:kmax)
                              ! 南北風速の拡張配列
                              ! Extended array of meridional velocity
    real(DP), intent(out) :: xyz_DPLon(0:imax-1, 1:jmax/2, 1:kmax)
                              ! 上流点の経度
                              ! Lon of departure point
    real(DP), intent(out) :: xyz_DPLat(0:imax-1, 1:jmax/2, 1:kmax)
                              ! 上流点の緯度
                              ! Lat of departure point

    !
    ! local variables
    !
    real(DP) :: xyz_MPLon (0:imax-1, 1:jmax/2, 1:kmax)
                              ! 中間点の経度
                              ! Lon of mid-point
    real(DP) :: xyz_MPLat (0:imax-1, 1:jmax/2, 1:kmax)
                              ! 中間点の緯度
                              ! Lat of mid-point
    real(DP) :: xyz_MPLonN(0:imax-1, 1:jmax/2, 1:kmax)
                              ! 中間点の経度(現在推定値)
                              ! Lon of mid-point (temporal estimation)
    real(DP) :: xyz_MPLatN(0:imax-1, 1:jmax/2, 1:kmax)
                              ! 中間点の緯度(現在推定値)
                              ! Lat of mid-point (now)
    logical  :: FlagWindInt   ! 初期化フラッグ
                              ! Flag for initialization
    logical  :: xyz_FlagEstDP(0:imax-1, 1:jmax/2, 1:kmax)

    integer:: i               ! 東西方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in zonal direction
    integer:: j               ! 南北方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in meridional direction
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer :: t              ! 推定回数の DO ループ用作業変数
                              ! Work variables for DO loop for DP estimation


    ! 初期化処理
    ! initialization
    !
    xyz_FlagEstDP = .true.
    do k = 1, kmax
      do j = 1, jmax/2
        do i = 0, imax-1
          xyz_MPLonN(i,j,k) = x_Lon(i)
          xyz_MPLatN(i,j,k) = y_Lat(j)
        end do
      end do
    end do


    FlagWindInt = .false.
    mp_search_loop : do t = 1, nloop_dp_h
      if( mod( t, 5 ) .eq. 0 ) write( 6, * ) "### Loop for searching departure points : ", t


      xyz_MPLon = xyz_MPLonN
      xyz_MPLat = xyz_MPLatN


      call SLTTDPHorCore( DelTime, x_Lon, y_Lat, y_SinLat, y_CosLat, iexmin, iexmax, jexmin, jexmax, x_ExtLon, y_ExtLat, xyz_ExtU, xyz_ExtV, xyz_MPLon, xyz_MPLat, xyz_MPLonN, xyz_MPLatN, xyz_FlagEstDP, FlagWindInt )
      FlagWindInt = .true.


    end do mp_search_loop

    xyz_FlagEstDP = .true.
    FlagWindInt   = .true.

    call SLTTDPHorCore( 2.0_DP * DelTime, x_Lon, y_Lat, y_SinLat, y_CosLat, iexmin, iexmax, jexmin, jexmax, x_ExtLon, y_ExtLat, xyz_ExtU, xyz_ExtV, xyz_MPLon, xyz_MPLat, xyz_DPLon, xyz_DPLat, xyz_FlagEstDP, FlagWindInt )


  end subroutine SLTTDPHor
Subroutine :
DelTime :real(DP), intent(in )
: $Delta t$
xyr_SigmaDot(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in )
: 鉛直流速(SigmaDot)
xyz_DPSigma(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(out)
: 上流点高度 Sigma of the departure point

鉛直方向の上流点探索 Finding DP in Vertical

[Source]

  subroutine SLTTDPVer( DelTime, xyr_SigmaDot, xyz_DPSigma )
    ! 鉛直方向の上流点探索
    ! Finding DP in Vertical

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only : r_Sigma, z_Sigma

    use sltt_const, only : nloop_dp_v


    real(DP), intent(in ) :: DelTime
                              ! $\Delta t$
    real(DP), intent(in ) :: xyr_SigmaDot(0:imax-1, 1:jmax, 0:kmax)
                              ! 鉛直流速(SigmaDot)    
    real(DP), intent(out) :: xyz_DPSigma (0:imax-1, 1:jmax, 1:kmax)
                              ! 上流点高度
                              ! Sigma of the departure point        


    !
    ! local variables
    !
    real(DP) :: xyz_MPSigma   (0:imax-1, 1:jmax, 1:kmax)
                              ! 中間点の高度
                              ! Sigma of mid-point
    real(DP) :: xyz_MPSigmaN  (0:imax-1, 1:jmax, 1:kmax)
                              ! 中間点の高度(現在推定値)
                              ! Lat of mid-point (temporal estimation)
    real(DP) :: xyz_MPSigmaDot(0:imax-1, 1:jmax, 1:kmax)
                              ! 中間点のSigmaDot
                              ! SigmaDot of mid-point
    real(DP) :: xyza_lcifz(0:imax-1, 1:jmax, 1:kmax, -1:2)
                              ! ラグランジュ補間のための作業変数
                              ! Work variable for Lagrange interpolation
    integer :: xyz_kk(0:imax-1, 1:jmax, 1:kmax)
                              ! ラグランジュ補間のための作業変数
                              ! Work variable for Lagrange interpolation
    

    integer:: i               ! 東西方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in zonal direction
    integer:: j               ! 南北方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in meridional direction
    integer:: k, kk, k2       ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer :: t              ! 推定回数の DO ループ用作業変数
                              ! Work variables for DO loop for DP estimation


    ! 初期化処理
    ! initialization
    !
    do k = 1, kmax
      xyz_MPSigmaN(:,:,k) = z_Sigma(k)
    end do

    ! 上流点推定のループ
    ! Loop for finding departure point
    mp_search_loop : do t = 1, nloop_dp_v


      xyz_MPSigma = xyz_MPSigmaN


      ! 中間点の鉛直流速をラグランジュ3次補間で求める 
      ! vertical wind velocity at mid-point is estimated by using Lagrange cubic 
      ! interpolation
      !
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1

            !
            ! Routine for dcpam
            !
            ! Departure points, xyz_MPSigma(:,:,k), must be located between 
            ! r_Sigma(kk) > xyz_MPSigma(k) > r_Sigma(kk+1).
            ! Further, 1 <= kk <= kmax-2.
            !

            !
            ! economical method
            !
            if( xyz_MPSigma(i,j,k) > r_Sigma(k) ) then
              k_search_1 : do k2 = k, 1, -1
                if( r_Sigma(k2) > xyz_MPSigma(i,j,k) ) exit k_search_1
              end do k_search_1
              xyz_kk(i,j,k) = k2
            else
              k_search_2 : do k2 = min( k+1, kmax ), kmax
                if( r_Sigma(k2) < xyz_MPSigma(i,j,k) ) exit k_search_2
              end do k_search_2
              xyz_kk(i,j,k) = k2 - 1
            end if
            xyz_kk(i,j,k) = min( max( xyz_kk(i,j,k), 1 ), kmax-2 )

          end do
        end do
      end do
      !
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            kk = xyz_kk(i,j,k)
            xyza_lcifz(i,j,k,-1) = ( xyz_MPSigma(i,j,k) - r_Sigma(kk  ) ) * ( xyz_MPSigma(i,j,k) - r_Sigma(kk+1) ) * ( xyz_MPSigma(i,j,k) - r_Sigma(kk+2) ) / ( ( r_Sigma(kk-1) - r_Sigma(kk  ) ) * ( r_Sigma(kk-1) - r_Sigma(kk+1) ) * ( r_Sigma(kk-1) - r_Sigma(kk+2) ) )
            xyza_lcifz(i,j,k, 0) = ( xyz_MPSigma(i,j,k) - r_Sigma(kk-1) ) * ( xyz_MPSigma(i,j,k) - r_Sigma(kk+1) ) * ( xyz_MPSigma(i,j,k) - r_Sigma(kk+2) ) / ( ( r_Sigma(kk  ) - r_Sigma(kk-1) ) * ( r_Sigma(kk  ) - r_Sigma(kk+1) ) * ( r_Sigma(kk  ) - r_Sigma(kk+2) ) )
            xyza_lcifz(i,j,k, 1) = ( xyz_MPSigma(i,j,k) - r_Sigma(kk-1) ) * ( xyz_MPSigma(i,j,k) - r_Sigma(kk  ) ) * ( xyz_MPSigma(i,j,k) - r_Sigma(kk+2) ) / ( ( r_Sigma(kk+1) - r_Sigma(kk-1) ) * ( r_Sigma(kk+1) - r_Sigma(kk  ) ) * ( r_Sigma(kk+1) - r_Sigma(kk+2) ) )
            xyza_lcifz(i,j,k, 2) = ( xyz_MPSigma(i,j,k) - r_Sigma(kk-1) ) * ( xyz_MPSigma(i,j,k) - r_Sigma(kk  ) ) * ( xyz_MPSigma(i,j,k) - r_Sigma(kk+1) ) / ( ( r_Sigma(kk+2) - r_Sigma(kk-1) ) * ( r_Sigma(kk+2) - r_Sigma(kk  ) ) * ( r_Sigma(kk+2) - r_Sigma(kk+1) ) )
          end do
        end do
      end do
      !
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            kk = xyz_kk(i,j,k)
            xyz_MPSigmaDot(i,j,k) = xyza_lcifz(i,j,k,-1) * xyr_SigmaDot(i,j,kk-1) + xyza_lcifz(i,j,k, 0) * xyr_SigmaDot(i,j,kk  ) + xyza_lcifz(i,j,k, 1) * xyr_SigmaDot(i,j,kk+1) + xyza_lcifz(i,j,k, 2) * xyr_SigmaDot(i,j,kk+2)
          end do
        end do
      end do



      do k = 1, kmax
        xyz_MPSigmaN(:,:,k) = z_Sigma(k) - xyz_MPSigmaDot(:,:,k) * DelTime
      end do
      xyz_MPSigmaN = min( max( xyz_MPSigmaN, r_Sigma(kmax) ), r_Sigma(0) )


    end do mp_search_loop

    ! 上流点高度の推定
    ! estimating departure point
    !
    do k = 1, kmax
      xyz_DPSigma(:,:,k) = z_Sigma(k) - xyz_MPSigmaDot(:,:,k) * 2.0_DP * DelTime
    end do


  end subroutine SLTTDPVer

Private Instance methods

Subroutine :
DelTime :real(DP), intent(in )
: $Delta t$
x_APLon(0:imax-1) :real(DP), intent(in )
: 到着点(グリッド上)の経度 Lon of arrival point (which is on grid)
y_APLat(1:jmax/2 ) :real(DP), intent(in )
: 到着点(グリッド上)の緯度 Lat of arrival point (which is on grid)
y_APSinLat(1:jmax/2) :real(DP), intent(in )
: 到着点の sin(緯度) sin(lat) of arrival point
y_APCosLat(1:jmax/2) :real(DP), intent(in )
: 到着点の cos(緯度) cos(lat) of arrival point
iexmin :integer , intent(in )
iexmax :integer , intent(in )
jexmin :integer , intent(in )
jexmax :integer , intent(in )
x_ExtLon(iexmin:iexmax) :real(DP), intent(in )
: 経度の拡張配列 Extended array of Lon
y_ExtLat(jexmin:jexmax) :real(DP), intent(in )
: 緯度の拡張配列 Extended array of Lat
xyz_ExtU(iexmin:iexmax, jexmin:jexmax, 1:kmax) :real(DP), intent(in )
: 東西風速の拡張配列 Extended array of zonal wind
xyz_ExtV(iexmin:iexmax, jexmin:jexmax, 1:kmax) :real(DP), intent(in )
: 南北風速の拡張配列 Extended array of meridional wind
xyz_MPLon(0:imax-1, 1:jmax/2, 1:kmax) :real(DP), intent(in )
: 中間点の経度 Lon of mid-point
xyz_MPLat(0:imax-1, 1:jmax/2, 1:kmax) :real(DP), intent(in )
: 中間点の緯度 Lat of mid-point
xyz_DPLon(0:imax-1, 1:jmax/2, 1:kmax) :real(DP), intent(inout)
: 上流点の経度 Lon of departure point
xyz_DPLat(0:imax-1, 1:jmax/2, 1:kmax) :real(DP), intent(inout)
: 上流点の緯度 Lat of departure point
xyz_FlagEstDP(0:imax-1, 1:jmax/2, 1:kmax) :logical , intent(in )
: 上流点探索のフラグ Flag for finding departure point
FlagWindInt :logical , intent(in )
: 初回フラグ Initial flag

水平方向の上流点探索のコア部分 Finding DP in Horizontal (Core)

[Source]

  subroutine SLTTDPHorCore( DelTime, x_APLon, y_APLat, y_APSinLat, y_APCosLat, iexmin, iexmax, jexmin, jexmax, x_ExtLon, y_ExtLat, xyz_ExtU, xyz_ExtV, xyz_MPLon, xyz_MPLat, xyz_DPLon, xyz_DPLat, xyz_FlagEstDP, FlagWindInt )
    ! 水平方向の上流点探索のコア部分
    ! Finding DP in Horizontal (Core)

    use constants0 , only : PI
    use constants  , only : RPlanet
    use mpi_wrapper, only : myrank
    use sltt_const , only : PIx2, PIH
    use sltt_lagint, only : SLTTLagIntCubCalcFactHor , SLTTLagIntCubIntHor


    real(DP), intent(in   ) :: DelTime
                              ! $\Delta t$    
    real(DP), intent(in   ) :: x_APLon(0:imax-1)
                              ! 到着点(グリッド上)の経度
                              ! Lon of arrival point (which is on grid)
    real(DP), intent(in   ) :: y_APLat(1:jmax/2 )
                              ! 到着点(グリッド上)の緯度
                              ! Lat of arrival point (which is on grid)
    real(DP), intent(in   ) :: y_APSinLat(1:jmax/2)
                              ! 到着点の sin(緯度)
                              ! sin(lat) of arrival point
    real(DP), intent(in   ) :: y_APCosLat(1:jmax/2)
                              ! 到着点の cos(緯度)
                              ! cos(lat) of arrival point    
    integer , intent(in   ) :: iexmin
    integer , intent(in   ) :: iexmax
    integer , intent(in   ) :: jexmin
    integer , intent(in   ) :: jexmax
    real(DP), intent(in   ) :: x_ExtLon(iexmin:iexmax)
                              ! 経度の拡張配列
                              ! Extended array of Lon
    real(DP), intent(in   ) :: y_ExtLat(jexmin:jexmax)
                              ! 緯度の拡張配列
                              ! Extended array of Lat
    real(DP), intent(in   ) :: xyz_ExtU(iexmin:iexmax, jexmin:jexmax, 1:kmax)
                              ! 東西風速の拡張配列
                              ! Extended array of zonal wind
    real(DP), intent(in   ) :: xyz_ExtV(iexmin:iexmax, jexmin:jexmax, 1:kmax)
                              ! 南北風速の拡張配列
                              ! Extended array of meridional wind
    real(DP), intent(in   ) :: xyz_MPLon(0:imax-1, 1:jmax/2, 1:kmax)
                              ! 中間点の経度
                              ! Lon of mid-point
    real(DP), intent(in   ) :: xyz_MPLat(0:imax-1, 1:jmax/2, 1:kmax)
                              ! 中間点の緯度
                              ! Lat of mid-point
    real(DP), intent(inout) :: xyz_DPLon(0:imax-1, 1:jmax/2, 1:kmax)
                              ! 上流点の経度
                              ! Lon of departure point    
    real(DP), intent(inout) :: xyz_DPLat(0:imax-1, 1:jmax/2, 1:kmax)
                              ! 上流点の緯度
                              ! Lat of departure point    
    logical , intent(in   ) :: xyz_FlagEstDP(0:imax-1, 1:jmax/2, 1:kmax)
                              ! 上流点探索のフラグ
                              ! Flag for finding departure point
    logical , intent(in   ) :: FlagWindInt
                              ! 初回フラグ
                              ! Initial flag

    !
    ! local variables
    !
    real(DP) :: xyz_MPU(0:imax-1, 1:jmax/2, 1:kmax)
                              ! 中間点の東西風速
                              ! Zonal wind at mid-point
    real(DP) :: xyz_MPV(0:imax-1, 1:jmax/2, 1:kmax)
                              ! 中間点の南北風速
                              ! Meridional wind at mid-point
    real(DP) :: MPCosLat      ! 中間点のcos(lat)
                              ! cos(lat) at mid-point

    real(DP) :: xyza_lcifx(0:imax-1, 1:jmax/2, 1:kmax, -1:2)
    real(DP) :: xyza_lcify(0:imax-1, 1:jmax/2, 1:kmax, -1:2)
    integer :: xyz_ii(0:imax-1, 1:jmax/2, 1:kmax)
    integer :: xyz_jj(0:imax-1, 1:jmax/2, 1:kmax)
                              ! ラグランジュ補間のための作業変数
                              ! Working variables for Lagrange interpolation


    integer:: i               ! 東西方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in zonal direction
    integer:: j               ! 南北方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in meridional direction
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数

    real(DP) :: InvRPlanet 


    InvRPlanet = 1.0_DP/RPlanet


    if( FlagWindInt ) then

      ! ラグランジュ3次補間
      ! Cubic interpolation
      call SLTTLagIntCubCalcFactHor( iexmin, iexmax, jexmin, jexmax, x_ExtLon, y_ExtLat, xyz_MPLon, xyz_MPLat, xyz_ii, xyz_jj, xyza_lcifx, xyza_lcify )
      ! 中間点の風速推定
      ! estimating wind velocities at mid-point
      !
      call SLTTLagIntCubIntHor( iexmin, iexmax, jexmin, jexmax, xyz_ii, xyz_jj, xyza_lcifx, xyza_lcify, xyz_ExtU, xyz_MPU )
      call SLTTLagIntCubIntHor( iexmin, iexmax, jexmin, jexmax, xyz_ii, xyz_jj, xyza_lcifx, xyza_lcify, xyz_ExtV, xyz_MPV )
    else

      !
      ! 推定風速の更新
      !
      xyz_MPU(0:imax-1,1:jmax/2,:) = xyz_ExtU(0:imax-1,1:jmax/2,:)
      xyz_MPV(0:imax-1,1:jmax/2,:) = xyz_ExtV(0:imax-1,1:jmax/2,:)

    end if


    do k = 1, kmax
      do j = 1, jmax/2
        do i = 0, imax-1

          if( xyz_FlagEstDP(i,j,k) ) then

!            MPCosLat = cos( xyz_MPLat(i,j,k) )
            xyz_DPLon(i,j,k) = x_APLon(i) - DelTime * xyz_MPU(i,j,k) / cos( xyz_MPLat(i,j,k) ) * InvRPlanet !/ RPlanet
            xyz_DPLat(i,j,k) = y_APLat(j) - DelTime * xyz_MPV(i,j,k)* InvRPlanet ! / RPlanet

          end if

        end do
      end do
    end do

    ! 普通の緯度経度の範囲に直す
    do k = 1, kmax
      do j = 1, jmax/2
        do i = 0, imax-1
          if( xyz_DPLat(i,j,k) < -PIH ) then
            xyz_DPLon(i,j,k) = xyz_DPLon(i,j,k) + PI
            xyz_DPLat(i,j,k) = -PIH + ( -PIH - xyz_DPLat(i,j,k) )
          else if( xyz_DPLat(i,j,k) > PIH ) then
            xyz_DPLon(i,j,k) = xyz_DPLon(i,j,k) + PI
            xyz_DPLat(i,j,k) = PIH - ( xyz_DPLat(i,j,k) - PIH )
          end if
          xyz_DPLon(i,j,k) = mod( xyz_DPLon(i,j,k) + PIx2, PIx2 )
        end do
      end do
    end do


  end subroutine SLTTDPHorCore