物理過程 陰解法用モジュール

  subroutine physics_implicit_fluxcorrection ( xyr_VelLonFlux, xyr_VelLatFlux, xyr_TempFlux, xyr_QvapFlux, xyz_DVerdiffVelLonDt, xyz_DVerdiffVelLatDt, xyz_DVerdiffTempDt, xyz_DVerdiffSurfTempDt, xyz_DVerdiffQvapDt, xyzo_VelMatrix, xyzo_TempMatrix, xyzo_QvapMatrix, xy_SurfVelMatrix, xyoo_SurfTempMatrix, xyoo_SurfQvapMatrix, DelTimePhy )
    !
    ! フラックスの補正
    ! 
    ! 
    use type_mod,    only: DBKIND, INTKIND, STRING
    use grid_3d_mod, only: im, jm, km
    use constants_mod, only: EL, Cp 
    use dc_trace,  only: BeginSub, EndSub
    real(DBKIND), intent(inout) :: xyr_VelLonFlux(im,jm,km+1) ! Ｕのフラックス
    real(DBKIND), intent(inout) :: xyr_VelLatFlux(im,jm,km+1) ! Ｖのフラックス
    real(DBKIND), intent(inout) :: xyr_TempFlux(im,jm,km+1) ! Ｔのフラックス
    real(DBKIND), intent(inout) :: xyr_QvapFlux(im,jm,km+1) ! ｑのフラックス
    real(DBKIND), intent(in) :: xyz_DVerdiffVelLonDt(im,jm,km)
                                                      ! 東西運動量変化項ＵＡ
    real(DBKIND), intent(in) :: xyz_DVerdiffVelLatDt(im,jm,km) 
                                                      ! 南北運動量変化項ＶＡ
    real(DBKIND), intent(in) :: xyz_DVerdiffTempDt(im,jm,km) ! 温度時間変化項Ｈ
    real(DBKIND), intent(in) :: xyz_DVerdiffSurfTempDt(im,jm) ! 地表温度変化率
    real(DBKIND), intent(in) :: xyz_DVerdiffQvapDt(im,jm,km) ! 比湿時間変化項Ｒ
    real(DBKIND), intent(in) :: xyzo_VelMatrix(im,jm,km,-1:1) ! ｕ陰解行列
    real(DBKIND), intent(in) :: xyzo_TempMatrix(im,jm,0:km,-1:1) ! Ｔ陰解行列
    real(DBKIND), intent(in) :: xyzo_QvapMatrix(im,jm,km,-1:1) ! ｑ陰解行列
    real(DBKIND), intent(in) :: xy_SurfVelMatrix(im,jm) ! ｕ陰解行列：地表
    real(DBKIND), intent(in) :: xyoo_SurfTempMatrix(im,jm,0:1,-1:1) 
                                                       ! Ｔ陰解行列：地表
    real(DBKIND), intent(in) :: xyoo_SurfQvapMatrix(im,jm,0:1,-1:1)
                                                       ! ｑ陰解行列：地表
    real(DBKIND), intent(in) :: DelTimePhy ! 時間刻みΔt

    real(DBKIND) :: ELF
    INTEGER(INTKIND) :: k
    character(STRING), parameter:: subname = "physics_implicit_fluxcorrection"

    ! 開始処理
    call BeginSub(subname)

    ELF = EL/Cp

    DO k = 2, km
      xyr_VelLonFlux(:,:,k) = xyr_VelLonFlux(:,:,K) - (   xyzo_VelMatrix(:,:,k,-1)* xyz_DVerdiffVelLonDt(:,:,k-1) - xyzo_VelMatrix(:,:,k-1,1)* xyz_DVerdiffVelLonDt(:,:,k) ) * DelTimePhy

      xyr_VelLatFlux(:,:,k) = xyr_VelLatFlux(:,:,k) - (   xyzo_VelMatrix(:,:,k,-1) * xyz_DVerdiffVelLatDt(:,:,k-1) - xyzo_VelMatrix(:,:,k-1,1) * xyz_DVerdiffVelLatDt(:,:,k) ) * DelTimePhy

      xyr_TempFlux(:,:,k) = xyr_TempFlux(:,:,k) - (   xyzo_TempMatrix(:,:,k,-1) * xyz_DVerdiffTempDt(:,:,k-1) - xyzo_TempMatrix(:,:,k-1,1) * xyz_DVerdiffTempDt(:,:,k) ) * DelTimePhy

      xyr_QvapFlux(:,:,k) = xyr_QvapFlux(:,:,k) - (  xyzo_QvapMatrix(:,:,k,-1) * xyz_DVerdiffQvapDt(:,:,k-1) - xyzo_QvapMatrix(:,:,k-1,1) * xyz_DVerdiffQvapDt(:,:,k) ) * DelTimePhy * ELF
     end do

     xyr_VelLonFlux(:,:,1) = xyr_VelLonFlux(:,:,1) - xy_SurfVelMatrix(:,:) * xyz_DVerdiffVelLonDt(:,:,1) * DelTimePhy

     xyr_VelLatFlux(:,:,1) = xyr_VelLatFlux(:,:,1) - xy_SurfVelMatrix(:,:) * xyz_DVerdiffVelLatDt(:,:,1) * DelTimePhy

     xyr_TempFlux(:,:,1) = xyr_TempFlux(:,:,1) - (   xyoo_SurfTempMatrix(:,:,1,-1) * xyz_DVerdiffSurfTempDt(:,:) + xyoo_SurfTempMatrix(:,:,1,0) * xyz_DVerdiffTempDt(:,:,1) ) * DelTimePhy

     xyr_QvapFlux(:,:,1) = xyr_QvapFlux(:,:,1) - (  xyoo_SurfQvapMatrix(:,:,1,-1) * xyz_DVerdiffSurfTempDt(:,:) + xyoo_SurfQvapMatrix(:,:,1,0) * xyz_DVerdiffQvapDt(:,:,1) * ELF ) * DelTimePhy

    ! 終了処理
    call EndSub(subname)
     
  end subroutine physics_implicit_fluxcorrection

  subroutine physics_implicit_init( xyr_VelLonFlux       , xyr_VelLatFlux       , xyr_TempFlux         , xyr_QvapFlux         , xyzo_VelMatrix       , xyzo_TempMatrix      , xyzo_QvapMatrix      , xyr_Press            , DelTimePhy           , xy_SurfHeatCapacity  , xy_SurfCondition       ) ! (in) 地表状態
    !
    use type_mod,    only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod, only: im, jm, km
    use constants_mod, only: Cp, Grav
    use dc_trace,    only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump

    implicit none
    real(DBKIND), intent(out) :: xyr_VelLonFlux(im*jm,km+1)        , xyr_VelLatFlux(im*jm,km+1)        , xyr_TempFlux(im*jm,km+1)          , xyr_QvapFlux(im*jm,km+1)          , xyzo_VelMatrix(im*jm,km,-1:1)     , xyzo_TempMatrix(im*jm,0:km,-1:1)  , xyzo_QvapMatrix(im*jm,km,-1:1)        !比湿陰解行列
    real(DBKIND), intent(in) :: xyr_Press(im*jm,km+1)             , DelTimePhy                        , xy_SurfHeatCapacity(im*jm)            ! 地表熱容量
    integer(INTKIND), intent(in) :: xy_SurfCondition(im*jm)               ! 地表状態
    character(STRING),  parameter:: subname = "physics_implicit_init"
    integer(INTKIND)    :: ij, k
        ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)

    continue

    !   開始処理
    call BeginSub(subname)

    !----------------------------------------------------------------
    !   陰解行列初期化
    !----------------------------------------------------------------

    ! ---- 1. 質量, 熱容量の項 ----
    xyzo_VelMatrix  = 0.0 
    xyzo_TempMatrix = 0.0 
    xyzo_QvapMatrix = 0.0 

    do k = 1, km 
       xyzo_VelMatrix(:,k,0)  = ( xyr_Press(:,k) - xyr_Press(:,k+1) ) / Grav / DelTimePhy  
       xyzo_TempMatrix(:,k,0) = xyzo_VelMatrix(:,k,0) * Cp
       xyzo_QvapMatrix(:,k,0) = xyzo_VelMatrix(:,k,0) * Cp
    end do

    do ij = 1, im*jm
       if ( xy_SurfCondition(ij) .GE. 1 ) then 
          xyzo_TempMatrix(ij,0,0) = xy_SurfHeatCapacity(ij) / DelTimePhy 
       else
          xyzo_TempMatrix(ij,0,0) = 1.
       end if
    end do

    ! ---- 2. フラックスをリセット ----
    xyr_VelLonFlux = 0.0 
    xyr_VelLatFlux = 0.0 
    xyr_TempFlux   = 0.0 
    xyr_QvapFlux   = 0.0 

    ! 終了処理
    call EndSub(subname)

  end subroutine physics_implicit_init

  subroutine physics_implicit_integrate( xyz_DVerdiffVelLonDt    , xyz_DVerdiffVelLatDt    , xyz_DVerdiffTempDt      , xyz_DVerdiffSurfTempDt  , xyz_DVerdiffQvapDt      , xyr_VelLonFlux          , xyr_VelLatFlux          , xyr_TempFlux            , xyr_SurfRadSFlux        , xyr_SurfRadLFlux        , xy_GroundTempFlux       , xyr_QvapFlux            , xyzo_VelMatrix          , xyzo_TempMatrix         , xyzo_QvapMatrix         , xy_SurfVelMatrix        , xyoo_SurfTempMatrix     , xyoo_SurfQvapMatrix     , xyo_SurfRadLMatrix      , DelTimePhy              , xy_SurfCondition          ) ! (in) 地表状態
    !
    ! 時間変化率の計算 (implicit)
    use type_mod,    only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod, only: im, jm, km
    use constants_mod, only: EL, Cp 
    use dc_trace,    only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump
    implicit none
    real(DBKIND), intent(out) :: xyz_DVerdiffVelLonDt(im*jm,km), xyz_DVerdiffVelLatDt(im*jm,km), xyz_DVerdiffTempDt(im*jm,km), xyz_DVerdiffSurfTempDt(im*jm), xyz_DVerdiffQvapDt(im*jm,km) ! 鉛直拡散加湿率
    real(DBKIND), intent(in) :: xyr_VelLonFlux(im*jm,km+1), xyr_VelLatFlux(im*jm,km+1), xyr_TempFlux(im*jm,km+1), xyr_SurfRadSFlux(im*jm), xyr_SurfRadLFlux(im*jm), xy_GroundTempFlux(im*jm)            , xyr_QvapFlux(im*jm,km+1)            , xyzo_VelMatrix(im*jm,km,-1:1)       , xyzo_TempMatrix(im*jm,0:km,-1:1)    , xyzo_QvapMatrix(im*jm,km,-1:1)      , xy_SurfVelMatrix(im*jm)             , xyoo_SurfTempMatrix(im*jm,0:1,-1:1) , xyoo_SurfQvapMatrix(im*jm,0:1,-1:1) , xyo_SurfRadLMatrix(im*jm,-1:1)      , DelTimePhy                              ! ２Δt
    integer(INTKIND) :: xy_SurfCondition(im*jm) ! 地表状態

    !----- 作業用内部変数 -----
    character(STRING),  parameter:: subname = "physics_implicit_integrate"
    integer(INTKIND)    :: ij, k, l
      ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
    real(DBKIND) :: xyz_DelTempQvap(im*jm,-km:km)            , xyzo_TempQvapLUMatrix(im*jm,-km:km,-1:1) , xyzo_VelLUMatrix(im*jm,km,-1:1)              ! ＬＵ行列

    continue

    ! 開始処理
    call BeginSub(subname)

    !----------------------------------------------------------------
    !   陰解行列計算
    !----------------------------------------------------------------

    ! ---- 1. 速度 (Vlon, Vlat) の解 ----

    xyzo_VelLUMatrix  = xyzo_VelMatrix
    xyzo_VelLUMatrix(:,1,0)  = xyzo_VelLUMatrix(:,1,0) + xy_SurfVelMatrix(:)

    call lu_decomposition_tridiagonal( xyzo_VelLUMatrix, im*jm, km )

    do k = 1, km
       xyz_DVerdiffVelLonDt(:,k) = xyr_VelLonFlux(:,k) - xyr_VelLonFlux(:,k+1)
       xyz_DVerdiffVelLatDt(:,k) = xyr_VelLatFlux(:,k) - xyr_VelLatFlux(:,k+1)
    end do

    call lu_solve_tridiagonal( xyz_DVerdiffVelLonDt , xyzo_VelLUMatrix     , 1, im*jm, km )

    call lu_solve_tridiagonal( xyz_DVerdiffVelLatDt , xyzo_VelLUMatrix     , 1, im*jm, km )

    ! ---- 2. 温度と比湿の解 ----

    do l = -1, 1

       do k = 1, km
          xyzo_TempQvapLUMatrix(:,k,l)   = xyzo_TempMatrix(:,k,l)
          xyzo_TempQvapLUMatrix(:,-k,-l) = xyzo_QvapMatrix(:,k,l)
       end do
       
       xyzo_TempQvapLUMatrix(:,1,l)   = xyzo_TempMatrix(:,1,l) + xyoo_SurfTempMatrix(:,1,l)
       xyzo_TempQvapLUMatrix(:,-1,-l) = xyzo_QvapMatrix(:,1,l) + xyoo_SurfQvapMatrix(:,1,l)

    end do

    xyzo_TempQvapLUMatrix(:,0,0) = xyzo_TempMatrix(:,0,0) + xyoo_SurfTempMatrix(:,0,0) + xyoo_SurfQvapMatrix(:,0,0) + xyo_SurfRadLMatrix(:,0)

    xyzo_TempQvapLUMatrix(:,0,1) = + xyoo_SurfTempMatrix(:,0,1) + xyo_SurfRadLMatrix(:,1)

    xyzo_TempQvapLUMatrix(:,0,-1) =  xyoo_SurfQvapMatrix(:,0,1) 

    call lu_decomposition_tridiagonal( xyzo_TempQvapLUMatrix, im*jm, 2*km+1 )

    do k = 1, km
       xyz_DelTempQvap(:,k)  = xyr_TempFlux(:,k) - xyr_TempFlux(:,k+1)
       xyz_DelTempQvap(:,-k) = xyr_QvapFlux(:,k) - xyr_QvapFlux(:,k+1)
    end do

    xyz_DelTempQvap(:,0) = - xyr_SurfRadSFlux  - xyr_SurfRadLFlux - xyr_TempFlux(:,1) - xyr_QvapFlux(:,1) + xy_GroundTempFlux


    call lu_solve_tridiagonal( xyz_DelTempQvap       , xyzo_TempQvapLUMatrix , 1, im*jm, 2*km+1 )

    ! ---- 2. 時間変化率 ----
    do k = 1, km
       xyz_DVerdiffVelLonDt(:,k) = xyz_DVerdiffVelLonDt(:,k) / DeltimePhy 
       xyz_DVerdiffVelLatDt(:,k) = xyz_DVerdiffVelLatDt(:,k) / DeltimePhy 
       xyz_DVerdiffTempDt(:,k)   = xyz_DelTempQvap(:,k) / DeltimePhy 
       xyz_DVerdiffQvapDt(:,k)   = xyz_DelTempQvap(:,-k) / DeltimePhy / EL * Cp
    end do
    
    do ij = 1, im*jm
       if ( xy_SurfCondition(ij) .GE. 1 ) then 
          xyz_DVerdiffSurfTempDt(ij) = xyz_DelTempQvap(ij,0) / DeltimePhy 
       else
          xyz_DVerdiffSurfTempDt(ij) = 0.
       end if
    end do

    !----------------------------------------------------------------
    !   終了処理
    !----------------------------------------------------------------
    call EndSub(subname)

  end subroutine physics_implicit_integrate