!= 火星のダスト巻き上げ
!
!= Mars Dust Lifting!

! Authors::   OGIHARA Hirotaka
! Version::   $Id: Mars_dust_lifting.f90,v 0.1 2012-07-17 16:18:31 yot Exp $ 
! Tag Name::  $Name: dcpam5-20120813-2$
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2009. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module mars_dust_dbug
  !
  != 火星のダスト巻き上げ
  !
  != Surface flux
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! Kahre et al (2006) に基づいて火星のダスト巻き上げ率を計算する. 
  !
  ! 
  !
  !== References
  !
  ! Melinda A. Kahre, James R. Murphy, and Robert M. Haberle
  ! Modeling the Martian dust cycle and surface dust
  ! reservoirs with the NASA Ames general circulation
  ! model, Journal of Geophysical Research, vol 111, E06008, 2006 
  !
  !== Procedures List
  !
  ! MarsDustLifting   :: 火星のダスト巻き上げ率の計算
  ! ------------      :: ------------
  !
  !== NAMELIST
  !
  ! NAMELIST#mars_dust_lifting_nml
  !

  ! モジュール引用 ; USE statements
  !
!  ! 格子点設定
!  ! Grid points settings
!  !
!  use gridset, only: imax, & ! 経度格子点数. 
!                             ! Number of grid points in longitude
!    &                jmax, & ! 緯度格子点数. 
!                             ! Number of grid points in latitude
!    &                kmax    ! 鉛直層数. 
!                             ! Number of vertical level
!
!
!  ! 種別型パラメタ
!  ! Kind type parameter
!  !
!  use dc_types, only: 8, &      ! 倍精度実数型. Double precision. 
!    &                 STRING     ! 文字列.       Strings. 
!
!  ! メッセージ出力
!  ! Message output
!  !
!  use dc_message, only: MessageNotify
  ! 宣言文 ; Declaration statements
  !
!  implicit none
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public :: DustLifting
  public :: MarsDustLiftingDD
  public :: DustLiftingDDtest

  ! 公開変数
  ! Public variables
  !

  ! 非公開変数
  ! Private variables
  !
  logical, save :: mars_dust_lifting_inited = .false.
                              ! 初期化設定フラグ
                              ! Initialization flag
  logical, save :: dbug = .true.
  real(8), save :: StressThreshold
  real(8), save :: KMHEfficiencyFactor
  real(8), save :: NewmanEfficiencyFactor
  real(8), save :: DDAEfficiencyFactor
  real(8), save :: DTHEfficiencyFactor
  real(8), save :: DustDevilEfficiencyFactor
  real(8), save :: VtangThreshold
  real(8), save :: FrictionRate
  integer, save :: IDMechanismDD
  integer, parameter :: IDMechanismDDA = 1
  integer, parameter :: IDMechanismDTH = 2
  integer, parameter :: IDMechanismDDALL = 3

contains

  !-------------------------------------------------------------------------------------

  subroutine DustLifting( &
      & bunkatu, Stress, & ! (in)
      & KMHMarsDustLifting, NewmanMarsDustLifting          & ! (out)
      & )
    !
    ! バルク法で求めた応力, 熱フラックスを使って火星のダスト巻き上げ率を
    ! 計算します.

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

    implicit none
    integer,intent(in):: bunkatu !分割数
    real(8), intent(in):: Stress(0:bunkatu) !応力
    real(8), intent(out):: KMHMarsDustLifting(0:bunkatu)
    real(8), intent(out):: NewmanMarsDustLifting(0:bunkatu)
    real(8),parameter:: Rho = 0.01d0
    real(8), parameter:: Grav = 3.72d0
    ! 作業変数
    ! Work variables
    !

    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 dimension of constituents

    ! 実行文 ; Executable statement
    ! 初期化確認
    ! Initialization check
    !
    !    if ( .not. mars_dust_lifting_inited ) call MarsDustLiftingInit
    StressThreshold = 0.0225d0
    KMHEfficiencyFactor = 0.1d0
    NewmanEfficiencyFactor = 9.0d0 * 0.000001d0
    !    StressThreshold = 22.5d0

    do i = 0, bunkatu
      if(Stress(i) < StressThreshold) then
        KMHMarsDustLifting(i) = 0.0d0
      else
        KMHMarsDustLifting(i) =                                                      & 
          & ((Stress(i))**2) * ((Stress(i) - StressThreshold)    &
          & / StressThreshold) ! KMH スキームの応力の項


        KMHMarsDustLifting(i) = (2.3d0 * 0.001d0) * KMHMarsDustLifting(i) ! 地球観測から求めた係数 2.3 x 10^-3 を掛ける

        !          KMHMarsDustLifting(i) = ((10.0d0)**6) * KMHMarsDustLifting(i)      ! 単位を kg m-1 s-1 にそろえる

        KMHMarsDustLifting(i) = KMHEfficiencyFactor * KMHMarsDustLifting(i)    ! 効率因子を掛けた

      end if
    end do

    do i = 0, bunkatu
      if(Stress(i) < StressThreshold) then
        NewmanMarsDustLifting(i) = 0.0d0
      else
        NewmanMarsDustLifting(i) =                                                      &
          & ((Stress(i))**1.5)/(Grav*sqrt(Rho)) * ( 1.0d0 - sqrt(StressThreshold/Stress(i))) &
          & * ( 1.0d0 + sqrt(StressThreshold/Stress(i)))**2 ! Newman スキームの応力の項

        NewmanMarsDustLifting(i) = 2.61d0 * NewmanMarsDustLifting(i) ! 風洞実験から求めた値 2.61を掛ける

        NewmanMarsDustLifting(i) = NewmanEfficiencyFactor * NewmanMarsDustLifting(i)    ! 効率因子を掛けた

      end if
    end do
  end subroutine DustLifting

  subroutine DustLiftingDDtest( &
     & bunkatu, Press_top, Heatflux, Temp, Ps, &
     & xy_B,                                              &
     & xy_MarsDustLiftingDD,                               & ! (out)
     & xy_AnotherMarsDustLiftingDD                         & ! (out) オプション
     & )

   ! モジュール引用 ; USE statements
   !
   ! 宣言文 ; Declaration statements
   !
   implicit none

   integer,intent(in):: bunkatu !分割数
   real(8), intent(in):: Ps
   real(8), intent(in):: Press_top(0:bunkatu) !応力
   real(8), intent(in):: Heatflux(0:bunkatu)
   real(8), intent(in):: Temp(0:bunkatu)
   real(8), intent(out) :: xy_B(0:bunkatu)
   real(8), intent(out):: xy_MarsDustLiftingDD(0:bunkatu,0:bunkatu)
   !火星のダスト巻き上げフラックス (kgm-2s-1)
   real(8), intent(out), optional:: xy_AnotherMarsDustLiftingDD(0:bunkatu,0:bunkatu)
   ! デバッグ用の別のスキームによる火星のダスト巻き上げフラックス (kgm-2s-1)
   real(8), parameter:: Grav = 3.72d0
   real(8), parameter:: CpDry = 843.9d0
   real(8), parameter:: MolWtDry = 43.5d-3
   real(8), parameter:: GasRUniv = 8.3144621d0
   real(8), parameter:: Surftemp = 224.0d0
   real(8), parameter:: Rho = 0.014d0
   real(8) GasRDry
   real(8) :: xy_Eta(0:bunkatu)
   ! ダストデビルの熱効率 (1 - b)
   real(8) :: xy_EtaH(0:bunkatu)
   ! ダストデビルの水平熱効率
   !  
   ! 地表面大気の平均温度に対する熱消滅のエントロピーで平均した温度の比
   real(8) :: xy_Vtang(0:bunkatu,0:bunkatu)
   ! ダストデビル渦の接線風速の二乗
   real(8) :: xy_DeltaP(0:bunkatu,0:bunkatu)
   ! ダストデビルの潜在的な強さ ( 渦壁と周りの圧力勾配)
   real(8) :: Chi
   ! R/c_p


   ! 作業変数
   ! Work variables
   !

   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 dimension of constituents

   ! 実行文 ; Executable statement

   GasRDry = GasRUniv/MolWtDry
   ! chi の計算
   Chi = GasRDry/CpDry
   !   対流境界層の上端の気圧を測定

   DDAEfficiencyFactor = 1.0d0 * 0.0000000001d0
   DTHEfficiencyFactor = 4.0d0 * 0.0000000001d0

   VtangThreshold = 15.0d0
   FrictionRate = 0.5d0

   !    do k = 0, kmax
   !      do i = 0, imax-1
   !        do j = 1, jmax
   !          if (k == 0) then
   !            if(xyr_HeatFlux(i,j,0) <= 0.0) then
   !              ktop(i,j) = k
   !            end if
   !          else 
   !            if (xyr_HeatFlux(i,j,k) == 0.0) then
   !              ktop(i,j) = k
   !            else if (xyr_HeatFlux(i,j,k) < 0.0) then
   !              ktop(i,j) = k - 1
   !              ! k での値がちょうど 0 でない場合は符号が入れ替わる前の k
   !              ! を最上端とする
   !            end if
   !          end if
   !        end do
   !      end do
   !    end do

   ! DDA スキームによるダスト巻き上げ

   do i = 0, bunkatu
     do j = 0, bunkatu
       if (j == 0) then
         xy_B(j) = 1.0d0
         xy_Eta(j) = 1.0d0 - xy_B(j)
       else
         xy_B(j) = (Ps**(Chi + 1.0d0) - Press_top(j)**(Chi + 1.0d0)) &
           & / ((Ps - Press_top(j)) * (Chi + 1.d0) * Ps**(Chi))
         xy_Eta(j) = 1.0d0 - xy_B(j)
       end if
       xy_MarsDustLiftingDD(i,j) = xy_Eta(j) * Heatflux(i)
       xy_MarsDustLiftingDD(i,j) = DDAEfficiencyFactor * xy_MarsDustLiftingDD(i,j)
     end do
   end do
   ! DTH によるダスト巻き上げスキーム(Newman スキーム)
   do i = 0, bunkatu
     do j = 0, bunkatu
       xy_EtaH(i) = (Temp(i) - Surftemp)/Surftemp
       ! これは Newman et al 2002 でのもの. Renno et al 1998 では $\frac{(T_0 - \bat{T_s})}{\bar{T_s}}$. $T_0$は大気最下層の温度. $\bar{T_s}$は地表面での平均温度. 
       xy_DeltaP(i,j) = Ps & 
         & * ( 1.0d0 - exp((FrictionRate * xy_Eta(i)) &
         & /(FrictionRate * xy_Eta(j) - 1.0d0) * ( xy_EtaH(i) / Chi)))

       xy_Vtang(i,j) = xy_DeltaP(i,j)/Rho

       xy_AnotherMarsDustLiftingDD(i,j) = max(0.0d0,((Rho * xy_Vtang(i,j) - VtangThreshold)/Grav))
       xy_AnotherMarsDustLiftingDD(i,j) = DTHEfficiencyFactor * xy_AnotherMarsDustLiftingDD(i,j)
     end do
   end do
  end subroutine DustLiftingDDtest

  subroutine MarsDustLiftingDD( &
      & imax, jmax, kmax,                                   & ! (in)
      & xy_SurfMajCompIce,                                  & ! (in)
      & xyr_Press, xyr_VirTemp,                             & ! (in)
      & xy_SurfTemp, xyr_HeatFlux,                          & ! (in)
      & xy_MarsDustLiftingDD,                               & ! (out)
      & xy_AnotherMarsDustLiftingDD                         & ! (out) オプション
      & )

    ! モジュール引用 ; USE statements
    !
    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(8), parameter:: Grav = 3.72d0
    real(8), parameter:: CpDry = 843.9d0
    real(8), parameter:: MolWtDry = 43.5d-3
    real(8), parameter:: GasRUniv = 8.3144621d0
    real(8), parameter:: CO2IceThreshold = 0.20d0
    real(8) GasRDry
    real(8), intent(in):: xy_SurfMajCompIce(0:imax-1, 1:jmax)
    ! $ M_mcs $ . 地表面主成分氷量 (kg m-2)
    ! Surface major component ice amount (kg m-2)
    real(8), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
    ! $ T_v $ . 仮温度 (半整数レベル). 
    ! Virtual temperature (half level)
    real(8), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
    ! $ p $ . 気圧 (半整数レベル). 
    ! pressure (half level)
    real(8), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
    ! $F_h$ . 熱フラックス(半整数レベル)
    ! Heat flux (half level)
    real(8), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
    ! $T_s$ . 地表面温度
    ! surface temperature
    real(8), intent(out):: xy_MarsDustLiftingDD(0:imax-1, 1:jmax)
    !火星のダスト巻き上げフラックス (kgm-2s-1)
    real(8), intent(out), optional:: xy_AnotherMarsDustLiftingDD(0:imax-1, 1:jmax)
    ! デバッグ用の別のスキームによる火星のダスト巻き上げフラックス (kgm-2s-1)
    real(8) :: MajCompIceThreshold
    ! 地表面で主成分が氷結しているかの閾値
    real(8) :: xy_SurfRho (0:imax-1, 1:jmax)
    ! 地表面密度
    real(8) :: xy_Eta(0:imax-1, 1:jmax)
    ! ダストデビルの熱効率 (1 - b)
    real(8) :: xy_EtaH(0:imax-1, 1:jmax)
    ! ダストデビルの水平熱効率
    real(8) :: xy_B(0:imax-1, 1:jmax)
    !  
    ! 地表面大気の平均温度に対する熱消滅のエントロピーで平均した温度の比
    real(8) :: xy_PressPBLTop(0:imax-1, 1:jmax)
    ! 対流境界層の上端での気圧
    real(8) :: xy_Vtang(0:imax-1, 1:jmax)
    ! ダストデビル渦の接線風速の二乗
    real(8) :: xy_DeltaP (0:imax-1,1:jmax)
    ! ダストデビルの潜在的な強さ ( 渦壁と周りの圧力勾配)
    real(8) :: Chi
    ! R/c_p


    ! 作業変数
    ! Work variables
    !

    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 dimension of constituents
    real(8):: ktop(0:imax-1,1:jmax)
    ! 境界層上端での層番号
    integer:: kk
    ! 境界層上端の層番号を入れる仮変数
    integer:: imax, jmax, kmax

    ! 実行文 ; Executable statement

    GasRDry = GasRUniv/MolWtDry
    !   地表面に氷が張っていると判断する閾値の設定
    MajCompIceThreshold = CO2IceThreshold

    ! chi の計算
    Chi = GasRDry/CpDry
    !   対流境界層の上端の気圧を測定

    IDMechanismDD = IDMechanismDDALL
    DDAEfficiencyFactor = 1.0d0 * 0.0000000001d0
    DTHEfficiencyFactor = 4.0d0 * 0.0000000001d0
    if (IDMechanismDD == IDMechanismDDA)then
    DustDevilEfficiencyFactor = DDAEfficiencyFactor
elseif(IDMechanismDD == IDMechanismDTH) then
    DustDevilEfficiencyFactor = DTHEfficiencyFactor
  else
    DustDevilEfficiencyFactor = 0.0d0
  end if

  VtangThreshold = 15.0d0
  FrictionRate = 0.5d0

    do k = 0, kmax
      do i = 0, imax-1
        do j = 1, jmax
          if (k == 0) then
            if(xyr_HeatFlux(i,j,0) <= 0.0) then
              ktop(i,j) = k
            end if
          else 
            if (xyr_HeatFlux(i,j,k) == 0.0) then
              ktop(i,j) = k
            else if (xyr_HeatFlux(i,j,k) < 0.0) then
              ktop(i,j) = k - 1
              ! k での値がちょうど 0 でない場合は符号が入れ替わる前の k
              ! を最上端とする
            end if
          end if
        end do
      end do
    end do

    select case ( IDMechanismDD )
    case ( IDMechanismDDA )
      ! DDA スキームによるダスト巻き上げ
      do i = 0, imax-1
        do j = 1, jmax
          if (xy_SurfMajCompIce(i,j) > MajCompIceThreshold ) then
            xy_MarsDustLiftingDD(i,j) = 0.0d0
          else 
            kk = ktop(i,j)
            if (kk == 0) then
              xy_B(i,j) = 0.0d0
            else
              write(*,*) "surface press",xyr_Press(i,j,0)
              write(*,*) "top press",xyr_Press(i,j,kk)
              write(*,*) "chi=",Chi
              xy_B(i,j) = (xyr_Press(i,j,0)**(Chi + 1.0d0) - xyr_Press(i,j,kk)**(Chi + 1)) &
                & / ((xyr_Press(i,j,0) - xyr_Press(i,j,kk)) * (Chi + 1) * xyr_Press(i,j,0)**(Chi))
              write(*,*) "b=",xy_B(i,j)
              xy_Eta(i,j) = 1.0d0 - xy_B(i,j)
              xy_MarsDustLiftingDD(i,j) = xy_Eta(i,j) * xyr_HeatFlux(i,j,0)
              xy_MarsDustLiftingDD(i,j) = DustDevilEfficiencyFactor * xy_MarsDustLiftingDD(i,j)
            end if
          end if
        end do
      end do
    case ( IDMechanismDTH )

      !地表面密度の計算
      do i = 0, imax-1
        do j = 1, jmax
          xy_SurfRho(i,j) = xyr_Press(i,j,0) /(GasRDry * xyr_VirTemp(i,j,0))
        end do
      end do
      ! DTH によるダスト巻き上げスキーム(Newman スキーム)

      do i = 0, imax-1
        do j = 1, jmax
          if (xy_SurfMajCompIce(i,j) > MajCompIceThreshold ) then
            xy_MarsDustLiftingDD(i,j) = 0.0d0
          else 
            kk = ktop(i,j)
            xy_B(i,j) = (xyr_Press(i,j,0)**(Chi + 1.0d0) - xyr_Press(i,j,kk)**(Chi + 1)) &
              & / ((xyr_Press(i,j,0) - xyr_Press(i,j,kk)) * (Chi + 1) * xyr_Press(i,j,0)**(Chi))
            xy_Eta(i,j) = 1.0d0 - xy_B(i,j)
            xy_EtaH(i,j) = (xyr_VirTemp(i,j,0) - xy_SurfTemp(i,j))/xy_SurfTemp(i,j)
            ! これは Newman et al 2002 でのもの. Renno et al 1998 では $\frac{(T_0 - \bat{T_s})}{\bar{T_s}}$. $T_0$は渦の中心の温度. $\bar{T_s}$は地表面層大気での平均温度. 
            xy_DeltaP(i,j) = xyr_Press(i,j,0) * &
              & ( 1.0d0 - exp((FrictionRate * xy_Eta(i,j))/& 
              & (FrictionRate * xy_Eta(i,j) - 1.0d0) * ( xy_EtaH(i,j) / Chi)))

            xy_Vtang(i,j) = xy_DeltaP(i,j)/xy_SurfRho(i,j)

            xy_MarsDustLiftingDD(i,j) = max(0.0d0,((xy_SurfRho(i,j) * xy_Vtang(i,j) - VtangThreshold)/Grav))
            xy_MarsDustLiftingDD(i,j) = DustDevilEfficiencyFactor * xy_MarsDustLiftingDD(i,j)
          end if
        end do
      end do

    case ( IDMechanismDDALL )
      ! デバッグ用両方のスキームの出力. DDA は MarsDustLifting へ DTH は AnotherMarsDustLifting へ
      ! DDA スキームによるダスト巻き上げ
      do i = 0, imax-1
        do j = 1, jmax
          if (xy_SurfMajCompIce(i,j) > MajCompIceThreshold ) then
            xy_MarsDustLiftingDD(i,j) = 0.0d0
          else 
            kk = ktop(i,j)
            xy_B(i,j) = (xyr_Press(i,j,0)**(Chi + 1.0d0) - xyr_Press(i,j,kk)**(Chi + 1)) &
              & / ((xyr_Press(i,j,0) - xyr_Press(i,j,kk)) * (Chi + 1) * xyr_Press(i,j,0)**(Chi))
            xy_Eta(i,j) = 1.0d0 - xy_B(i,j)
            xy_MarsDustLiftingDD(i,j) = xy_Eta(i,j) * xyr_HeatFlux(i,j,0)
            xy_MarsDustLiftingDD(i,j) = DDAEfficiencyFactor * xy_MarsDustLiftingDD(i,j)
          end if
        end do
      end do
      !地表面密度の計算
      do i = 0, imax-1
        do j = 1, jmax
          xy_SurfRho(i,j) = xyr_Press(i,j,0) /(GasRDry * xyr_VirTemp(i,j,0))
        end do
      end do
      ! DTH によるダスト巻き上げスキーム(Newman スキーム)
      do i = 0, imax-1
        do j = 1, jmax
          if (xy_SurfMajCompIce(i,j) > MajCompIceThreshold ) then
            xy_AnotherMarsDustLiftingDD(i,j) = 0.0d0
          else 
            kk = ktop(i,j)
            xy_EtaH(i,j) = (xyr_VirTemp(i,j,0) - xy_SurfTemp(i,j))/xy_SurfTemp(i,j)
            ! これは Newman et al 2002 でのもの. Renno et al 1998 では $\frac{(T_0 - \bat{T_s})}{\bar{T_s}}$. $T_0$は大気最下層の温度. $\bar{T_s}$は地表面での平均温度. 
            xy_DeltaP(i,j) = xyr_Press(i,j,0) & 
              & * ( 1.0d0 - exp((FrictionRate * xy_Eta(i,j)) &
              & /(FrictionRate * xy_Eta(i,j) - 1.0d0) * ( xy_EtaH(i,j) / Chi)))

            xy_Vtang(i,j) = xy_DeltaP(i,j)/xy_SurfRho(i,j)

            xy_AnotherMarsDustLiftingDD(i,j) = max(0.0d0,((xy_SurfRho(i,j) * xy_Vtang(i,j) - VtangThreshold)/Grav))
            xy_AnotherMarsDustLiftingDD(i,j) = DTHEfficiencyFactor * xy_AnotherMarsDustLiftingDD(i,j)
          end if
        end do
      end do
    end select

  end subroutine MarsDustLiftingDD
  !-------------------------------------------------------------------

end module mars_dust_dbug


program duststressdepend

  use mars_dust_dbug, only: DustLifting, DustLiftingDDtest, MarsDustLiftingDD
  implicit none

  real(8), allocatable:: xy_KMHMarsDustLifting(:)
  real(8), allocatable:: xy_NewmanMarsDustLifting(:)
  real(8), allocatable:: xy_B(:)
  real(8), allocatable:: Stress(:)
  real(8), allocatable:: MarsDustLiftingDDA(:,:)
  real(8), allocatable:: MarsDustLiftingDTH(:,:)
  real(8), allocatable:: Press_top(:)
  real(8), allocatable:: Temp(:)
  real(8), allocatable:: Heat(:)
!  real(8), allocatable:: xy_MarsDustLiftingDDA(:,:)
!  real(8), allocatable:: xy_MarsDustLiftingDTH(:,:)
!  real(8), allocatable:: HeatFlux(:,:,:)
!  real(8), allocatable:: Press(:,:,:)
!  real(8), allocatable:: Ps(:,:)
!  real(8), allocatable:: SurfTemp(:,:)
!  real(8), allocatable:: VirTemp(:,:,:)
!  real(8), allocatable:: Ice(:,:)
!  real(8) :: Sens
  real, allocatable:: KMHMarsDustLifting(:)
  real, allocatable:: NewmanMarsDustLifting(:)
  real, allocatable:: B(:)
  real, allocatable:: DDADustLifting(:,:)
  real, allocatable:: DTHDustLifting(:,:)
  real, allocatable:: Press(:)
  real, allocatable:: SurfStress(:)
  real, allocatable:: PressPBLtop(:)
  real, allocatable:: VirTemp(:)
  real, allocatable:: SurfSens(:)
  integer,parameter::bunkatu=100
!  real(8),parameter::Ps = 600.0d0
  real(8),parameter::Ps = 729.0d0
  ! lon 326 lat19 での大体の年間平均地表面圧力
  integer:: i,j,k
!  integer:: imax,jmax,kmax
!  integer:: num
!  character(136) :: text

  allocate(xy_KMHMarsDustLifting(0:bunkatu))
  allocate(xy_NewmanMarsDustLifting(0:bunkatu))
  allocate(KMHMarsDustLifting(0:bunkatu))
  allocate(NewmanMarsDustLifting(0:bunkatu))
  allocate(MarsDustLiftingDDA(0:bunkatu,0:bunkatu))
  allocate(MarsDustLiftingDTH(0:bunkatu,0:bunkatu))
  allocate(DDADustLifting(0:bunkatu,0:bunkatu))
  allocate(DTHDustLifting(0:bunkatu,0:bunkatu))
  allocate(Stress(0:bunkatu))
  allocate(SurfStress(0:bunkatu))
  allocate(Press_top(0:bunkatu))
  allocate(Temp(0:bunkatu))
  allocate(Heat(0:bunkatu))
  allocate(B(0:bunkatu))
  allocate(xy_B(0:bunkatu))
  allocate(Press(0:bunkatu))
  
!  open(20, file="./dat/heat.dat")
!  open(21, file="./dat/press.dat")
!  open(22, file="./dat/ps.dat")
!  open(23, file="./dat/surftemp.dat")
!  open(24, file="./dat/temp.dat")
!  open(25, file="./dat/ice.dat")
!  open(26, file="./dat/sens.dat")
!
!  read(24,'(a)') text
!  imax = count(transfer(trim(text),' ', size=len_trim(text)) ==' ')
!  if (imax == 0) then
!    imax = 1
!  endif
!  jmax = imax
!  rewind(24)
!
!  kmax = 0
!  rewind(24)
!  do 
!    read(24,*,end=458)
!    kmax = kmax + 1
!  end do
!  458 write(*,*) kmax
!  rewind(24)
!  allocate(xy_MarsDustLiftingDDA(0:imax-1,1:jmax))
!  allocate(xy_MarsDustLiftingDTH(0:imax-1,1:jmax))
!  allocate(HeatFlux(0:imax-1,1:jmax,0:kmax))
!  allocate(Press(0:imax-1,1:jmax,0:kmax))
!  allocate(Ps(0:imax-1,1:jmax))
!  allocate(SurfTemp(0:imax-1,1:jmax))
!  allocate(VirTemp(0:imax-1,1:jmax,0:kmax))
!  allocate(Ice(0:imax-1,1:jmax))
!  do k = 0,kmax-1
!    do i = 0,imax-1
!      do j = 1,jmax
!        read(20,*) HeatFlux(i,j,k)
!        read(21,*) Press(i,j,k)
!        write(*,*) k
!        read(24,*) VirTemp(i,j,k)
!      end do
!    end do
!  enddo
!  do i = 0,imax-1
!    do j = 1, jmax
!      read(22,*) Ps(i,j)
!      read(23,*) SurfTemp(i,j)
!      read(25,*) ice(i,j)
!      read(26,*) sens
!    end do
!  end do
!
!  HeatFlux(0,1,0) = sens

!  !実行文----------------------------------
!   ! 風応力
!  do i = 0, bunkatu
!    Stress(i) = Stress(i) + 0.0008d0 * i
!    !Stress(i) = Stress(i) + 0.8d0 * dble(i)
!  end do
!  call DustLifting(bunkatu,Stress,xy_KMHMarsDustLifting, xy_NewmanMarsDustLifting)
!  write(*,*) "1"
!  KMHMarsDustLifting= real(xy_KMHMarsDustLifting)
!  NewmanMarsDustLifting= real(xy_NewmanMarsDustLifting)
!  SurfStress = real(Stress)
!!  write(*,*) Stress(bunkatu)
!!  write(*,*) xy_KMHMarsDustLifting(28)
!!  write(*,*) KMHMarsDustLifting(29)
!  open(10,file='KMH.dat')
!  open(11,file='Newman.dat')
!  open(12,file='KMHandNewman.dat')
!  do i = 0, bunkatu
!    write(10,*) Stress(i), xy_KMHMarsDustLifting(i)
!    write(11,*) Stress(i), xy_NewmanMarsDustLifting(i)
!    write(12,*) Stress(i), xy_KMHMarsDustLifting(i), xy_NewmanMarsDustLifting(i)
!  end do 

  ! ダストデビル---------------------------
!  open(30, file="DDA.dat")
!  open(31, file="DTH.dat")
!  call MarsDustLiftingDD( & 
!      & imax, jmax, kmax,  & ! (in)
!      & Ice, Press, VirTemp, SurfTemp, HeatFlux, &
!      & xy_MarsDustLiftingDDA,                   & ! (out)
!      & xy_MarsDustLiftingDTH)
!  write(30,*) xy_MarsDustLiftingDDA
!  write(31,*) xy_MarsDustLiftingDTH
   Press_top(0) = Ps
   Temp(0) = 140.0d0 
   Heat(0) = 0.0d0
  do i = 1, bunkatu
    Press_top(i) = Press_top(i-1) - Ps/bunkatu 
    Temp(i) = Temp(i-1) + 1.6
    Heat(i) = Heat(i-1) + 8.0d0/bunkatu
  end do
  call DustLiftingDDtest( bunkatu, Press_top, Heat, Temp, Ps, xy_B, MarsDustLiftingDDA,MarsDustLiftingDTH)
  open(13,file='DDA.dat')
  open(14,file='DTH.dat')
  open(15,file='DDAandDTH.dat')
  open(16,file='B.dat')
  do i = 0, bunkatu
    do j = 0, bunkatu
      write(13,*) Press_top(i), Heat(j), MarsDustLiftingDDA(i,j)
      write(14,*) Press_top(i), Temp(j),  MarsDustLiftingDTH(i,j)
      write(15,*) Press_top(i), Heat(j), Temp(j), MarsDustLiftingDDA(i,j), MarsDustLiftingDTH(i,j)
    end do 
    write(16,*) Press_top(i), xy_B(i)
  end do

  DDADustLifting= real(MarsDustLiftingDDA)
  DTHDustLifting= real(MarsDustLiftingDTH)
  Press = real(Press_top)
  SurfSens = real(Heat)
  VirTemp = real(Temp)
  B = real(xy_B)


  !ファイルのクローズ
  close(10)
  close(11)
  close(12)
  close(13)
  close(14)
  close(15)
  close(16)
!  close(21)
!  close(22)
!  close(23)
!  close(24)
!  close(25)
!  close(26)
!  close(30)
!  close(31)

  !DCL-------------------------------------------------------
  !DCL of wind stress
  !call sgpwsn
  !!call gropn(2)
  !call gropn(4)
  !call grfrm
  !call ussttl("SurfStress","N m-2", "dust lifting rate","kg m-2 s-1")
  !call grswnd( 0.0, 0.08, 1.0e-9, 1.0e-5)
  !call grsvpt(0.2,0.8,0.2,0.8)
  !!call usspnt(72,SurfStress(29:bunkatu),KMHMarsDustLifting(29:bunkatu))
  !!call usspnt(72,SurfStress(29:bunkatu),NewmanMarsDustLifting(29:bunkatu))
  !call grstrn(2)
  !!call uspfit
  !call grstrf
  !call usdaxs
  !call uulin(bunkatu+1,SurfStress(0:bunkatu),KMHMarsDustLifting(0:bunkatu))
  !!call uuslnt(2)
  !!call uulin(bunkatu+1,SurfStress(0:bunkatu),NewmanMarsDustLifting(0:bunkatu))
  !!call uulin(72,SurfStress(29:bunkatu),KMHMarsDustLifting(29:bunkatu))
  !!call uuslnt(2)
  !!call uulin(72,SurfStress(29:bunkatu),NewmanMarsDustLifting(29:bunkatu))
  !call grcls
  !-----------------------------
  !DCL of dust devil
  call sgpwsn
  call gropn(2)
!  call gropn(4)
  call grfrm
  call ussttl("Sens","W m-2", "Press-top","m-1 s-2")
  call grswnd(0.0, 8.0, 729.0, 0.0)
!  call ussttl('Press-top','m-1 s-2','B',' ')
!  call ussttl('Press-top','m-1 s-2','eta',' ')
!  call grswnd( 729.0, 0.0, 0.0, 0.2)
  call grsvpt(0.2,0.8,0.2,0.8)
!  !call uspfit
  call grstrf
  call usdaxs
!  call uulin(bunkatu+1, Press(0:bunkatu),  B(0:bunkatu))
   ! B の 対流層上端の圧力に対する依存性
!  call uulin(bunkatu+1, Press(0:bunkatu), 1.0e0 - B(0:bunkatu))
   ! 鉛直熱効率 の 対流層上端の圧力に対する依存性
  call udcntr(DTHDustLifting,bunkatu+1,bunkatu+1,bunkatu+1)
  call grcls
end program duststressdepend
