subroutine physics_negq(xyz_Qvap, xyz_DNegQvapDt, xyr_Press, DelTimePhy)

    !==== Dependency
    use type_mod,    only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod, only: im, jm, km
    use spml_mod,    only: AvrLonLat_xy
    use dc_trace,    only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump

    implicit none

    !==== Input
    !
    real(DBKIND), intent(in) :: xyr_Press(im,jm,km+1) , DelTimePhy               ! intent(in): 2Δt

    !==== In/Out
    !
    real(DBKIND), intent(inout) :: xyz_Qvap(im,jm,km)     , xyz_DNegQvapDt(im,jm,km)     ! intent(inout): 比湿変化率

    !----- 作業用内部変数 -----
    character(STRING),  parameter:: subname = "physics_negq"
    real(DBKIND)        :: xyz_Qvap_b(im,jm,km)              , xyz_DPressDz(im,jm,km)            , xyz_Qvap_DPressDz(im,jm,km)       , xyz_DelQvap_DPressDz(im,jm,km)    , Qvap_DPressDz_AvrXYZ              , DelQvap_DPressDz_AvrXYZ               ! \int<ΔqΔp>dz

    real(DBKIND)        :: Work
    ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
    integer(INTKIND)    :: i,j, k


    continue

    !----------------------------------------------------------------
    !   開始処理
    !----------------------------------------------------------------
    call BeginSub(subname)

    !----------------------------------------------------------------
    !   負の水蒸気の補正
    !----------------------------------------------------------------

    !----- 調節前 Qvap の保存, ΔPress の計算 -----
    do k = 1, km
       do i = 1, im
          do j = 1, jm
             xyz_Qvap_b(i,j,k)   = xyz_Qvap(i,j,k)
             xyz_DPressDz(i,j,k) = xyr_Press(i,j,k) - xyr_Press(i,j,k+1)
          end do
       end do
    end do

    !----- 局所補正 ----- (最下層に負値を追いやる?)
    do k = km, 2, -1
       do i = 1, im
          do j = 1, jm
             if ( xyz_Qvap(i,j,k) < 0 ) then 
                Work = - xyz_Qvap(i,j,k) * xyz_DPressDz(i,j,k) / xyz_DPressDz(i,j,k-1) 
                xyz_Qvap(i,j,k) = 0.0d0
                xyz_Qvap(i,j,k-1) = xyz_Qvap(i,j,k-1) - Work
             end if
          end do
       end do
    end do

    !----- 全球補正 ----- 
    ! 各層における補正量の算出
    do k = 1, km
       do i = 1, im
          do j = 1, jm
             xyz_Qvap_DPressDz(i,j,k) = xyz_Qvap(i,j,k) * xyz_DPressDz(i,j,k) 
             if ( xyz_Qvap_DPressDz(i,j,k) < 0 ) then 
                xyz_DelQvap_DPressDz(i,j,k) = - xyz_Qvap_DPressDz(i,j,k)
             else
                xyz_DelQvap_DPressDz(i,j,k) = 0.0d0
             end if
          end do
       end do
    end do

    ! 補正量の東西南北鉛直平均
    Qvap_DPressDz_AvrXYZ = 0.0d0
    DelQvap_DPressDz_AvrXYZ = 0.0d0
    do k = 1, km
       Qvap_DPressDz_AvrXYZ    = Qvap_DPressDz_AvrXYZ + AvrLonLat_xy( xyz_Qvap_DPressDz(:,:,k) )
       DelQvap_DPressDz_AvrXYZ = DelQvap_DPressDz_AvrXYZ + AvrLonLat_xy( xyz_DelQvap_DPressDz(:,:,k) )
    end do


    ! 全体を引き下げる, 負値をゼロにする
    if ( Qvap_DPressDz_AvrXYZ .ne. 0.0d0 ) then 
       do k = 1, km
          do i = 1, im
             do j = 1, jm
                xyz_Qvap(i,j,k) = Qvap_DPressDz_AvrXYZ / (Qvap_DPressDz_AvrXYZ + DelQvap_DPressDz_AvrXYZ) * Max( xyz_Qvap(i,j,k) , 0.0d0 )             
             end do
          end do
       end do
    end if

    !----- 比湿変化の算出 ----- 
    do k = 1, km
       do i = 1, im
          do j = 1, jm
             xyz_DNegQvapDt(i,j,k) = xyz_DNegQvapDt(i,j,k) + ( xyz_Qvap(i,j,k) - xyz_Qvap_b(i,j,k) ) / DelTimePhy
          end do
       end do
    end do

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

  end subroutine physics_negq