! 繧ｵ繧ｦ繝ｳ繝�ぅ繝ｳ繧ｰ繝輔ぃ繧､繝ｫ縺九ｉ蛻晄悄蝣ｴ繧呈ｱｺ繧√ｋ縺溘ａ縺ｮ繧ｵ繝悶Ν繝ｼ繝√Φ
!
! Authors::   SUGIYAMA Koichiro, ODAKA Masatsugu
! Version::   $Id: initialdata_sounding.f90,v 1.8 2014/07/08 00:59:09 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!

module initialdata_sounding
  !
  ! 繧ｵ繧ｦ繝ｳ繝�ぅ繝ｳ繧ｰ繝輔ぃ繧､繝ｫ縺九ｉ蛻晄悄蝣ｴ繧呈ｱｺ繧√ｋ縺溘ａ縺ｮ繧ｵ繝悶Ν繝ｼ繝√Φ
  ! 繧ｵ繧ｦ繝ｳ繝�ぅ繝ｳ繧ｰ繝輔ぃ繧､繝ｫ縺ｯ繝�く繧ｹ繝医ヵ繧｡繧､繝ｫ縺ｧ譖ｸ縺九ｌ縺ｦ縺�ｋ.
  ! 蟆�擂逧�↓縺ｯ netCDF 縺ｫ螟画峩縺吶ｋ莠亥ｮ�.
  !
  ! 鬮伜ｺｦ譁ｹ蜷代�險育ｮ鈴�伜沺縺御ｸ弱∴繧峨ｌ縺溘ヵ繧｡繧､繝ｫ縺ｮ荳顔ｫｯ繧医ｊ繧るｫ倥＞蝣ｴ蜷医�,
  ! VIRA 繝��繧ｿ繧堤畑縺�※陬懷ｮ後☆繧�.
  !

  !繝｢繧ｸ繝･繝ｼ繝ｫ隱ｭ縺ｿ霎ｼ縺ｿ
  use dc_types, only: DP, STRING
  use gtool_history, only : GT_HISTORY

  !證鈴ｻ吶�蝙句ｮ｣險遖∵ｭ｢
  implicit none

  !蜈ｱ騾壼､画焚
  real(DP), save, private :: r_tmpAlt(10000)
  real(DP), save, private :: r_tmpTemp(10000)
  real(DP), save, private :: r_tmpPTemp(10000)
  real(DP), save, private :: r_tmpVelX(10000)
  real(DP), save, private :: r_tmpVelY(10000)
  real(DP), save, private :: TempTr = 0.0d0
  real(DP), save, private :: AltTr  = 0.0d0
  real(DP), save, private :: DelAlt = 4.0d3
  integer,  save, private :: NumRec = 0

  integer,  save, private :: AltCol   = 0
  integer,  save, private :: TempCol  = 0
  integer,  save, private :: PTempCol = 0
  integer,  save, private :: VelXCol  = 0
  integer,  save, private :: VelYCol  = 0

  integer , save, private :: NumVIRA = 23
  real(DP), save, private :: r_viraTemp(23)
  real(DP), save, private :: r_viraAlt(23)

  type(GT_HISTORY),  save, private :: rstat
  character(STRING), save, private :: OutputFile  = "BasicZ.Sounding.nc"

  !蜈ｬ髢句､画焚
  public  initialdata_sounding_init
  public  initialdata_sounding_basic
  public  initialdata_sounding_wind

contains

!!!------------------------------------------------------------------------------!!!
  subroutine initialdata_sounding_init
    !
    !繝輔ぃ繧､繝ｫ縺九ｉ繧ｵ繧ｦ繝ｳ繝�ぅ繝ｳ繧ｰ繝��繧ｿ繧定ｪｭ縺ｿ霎ｼ繧
    !

    !繝｢繧ｸ繝･繝ｼ繝ｫ隱ｭ縺ｿ霎ｼ縺ｿ
    use dc_types,      only: DP
    use dc_iounit,     only: FileOpen
    use namelist_util, only: namelist_filename

    !證鈴ｻ吶�蝙句ｮ｣險遖∵ｭ｢
    implicit none

    !蜀�Κ螟画焚
    integer             :: unit         !險ｭ螳壹ヵ繧｡繧､繝ｫ逕ｨ陬�ｽｮ逡ｪ蜿ｷ
    integer             :: io
    character(30)       :: SoundingFile

    integer, parameter  :: maxch=12
    character(len=100)  :: buf, eachcol(maxch)
    integer             :: num, MaxCol

    !險ｭ螳壹ヵ繧｡繧､繝ｫ縺九ｉ隱ｭ縺ｿ霎ｼ繧蜃ｺ蜉帙ヵ繧｡繧､繝ｫ諠��ｱ
    !
    NAMELIST /initialdata_sounding_nml/ SoundingFile, AltCol, TempCol, PTempCol, VelXCol, VelYCol, TempTr, AltTr, DelAlt

    !險ｭ螳壹ヵ繧｡繧､繝ｫ縺九ｉ隱ｭ縺ｿ霎ｼ縺ｿ
    !
    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=initialdata_sounding_nml)
    close(unit)

    ! 蛻晄悄蛹�
    !
    io = 0
    NumRec = 0
    MaxCol = max( AltCol, TempCol )
    r_tmpAlt  = 0.0d0; r_tmpTemp = 0.0d0
    r_tmpVelX = 0.0d0; r_tmpVelY = 0.0d0; r_tmpPTemp = 0.0d0;

    ! 繝輔ぃ繧､繝ｫ縺ｮ繧ｪ繝ｼ繝励Φ
    !
    open (17, file=SoundingFile, status='old')

    ! 繝輔ぃ繧､繝ｫ蜻ｼ縺ｳ蜃ｺ縺�
    !
    do while ( io == 0 )
      ! 1 陦悟�隱ｭ縺ｿ蜃ｺ縺�
      !
      read (17, '(a)', IOSTAT=io) buf

!!      write(*,*) buf

      ! 陦後ｒ繧ｫ繝ｳ繝槫玄蛻�ｊ縺ｧ蛻�牡
      !
      call devidecsv( buf, eachcol, maxch, num )

      ! num 縺ｮ蛟､縺悟ｰ上＆縺�ｂ縺ｮ縺ｯ繝倥ャ繝縺ｨ縺ｿ縺ｪ縺�.
      !
      if (num >= MaxCol) then
        ! 陦梧焚縺ｮ險育ｮ�
        !
        NumRec = NumRec + 1

        ! 蛟､縺ｮ莉｣蜈･
        !
        if (AltCol > 0)   read( eachcol(AltCol)(1:len_trim(eachcol(AltCol))), *)     r_tmpAlt(NumRec)
        if (TempCol > 0)  read( eachcol(TempCol)(1:len_trim(eachcol(TempCol))), *)   r_tmpTemp(NumRec)
        if (PTempCol > 0) read( eachcol(PTempCol)(1:len_trim(eachcol(PTempCol))), *) r_tmpPTemp(NumRec)
        if (VelXCol > 0)  read( eachcol(VelXCol)(1:len_trim(eachcol(VelXCol))), *)   r_tmpVelX(NumRec)
        if (VelYCol > 0)  read( eachcol(VelYCol)(1:len_trim(eachcol(VelYCol))), *)   r_tmpVelY(NumRec)

      end if

    end do

!!    write(*,*) r_tmpAlt(1:100)
!!    write(*,*) r_tmpTemp(1:100)

    ! 繝輔ぃ繧､繝ｫ縺ｮ繧ｯ繝ｭ繝ｼ繧ｺ
    !
    close (17)

    ! VIRA 繝��繧ｿ
    ! 譛ｬ譚･縺ｯ繝輔ぃ繧､繝ｫ縺九ｉ荳弱∴繧九∋縺阪□縺後ユ繧ｹ繝医→縺励※.
    !
    r_viraTemp = (/241.0d0, 235.4d0, 229.8d0, 224.1d0, 218.6d0, &
      &            212.1d0, 205.3d0, 197.1d0, 189.9d0, 183.8d0, &
      &            178.2d0, 173.6d0, 169.4d0, 167.2d0, 167.2d0, &
      &            169.2d0, 172.0d0, 175.4d0, 178.0d0, 189.0d0, &
      &            206.0d0, 225.0d0, 246.5d0 /)
    r_viraAlt = (/ 66000.0d0,  68000.0d0,  70000.0d0,  72000.0d0,  74000.0d0, &
      &            76000.0d0,  78000.0d0,  80000.0d0,  82000.0d0,  84000.0d0, &
      &            86000.0d0,  88000.0d0,  90000.0d0,  92000.0d0,  94000.0d0, &
      &            96000.0d0,  98000.0d0, 100000.0d0, 110000.0d0, 120000.0d0, &
      &           130000.0d0, 140000.0d0, 150000.0d0 /)

  end subroutine initialdata_sounding_init

!!!------------------------------------------------------------------------------!!!
  subroutine  initialdata_sounding_basic( z_Temp, z_Press )
    !
    ! 繧ｵ繧ｦ繝ｳ繝�ぅ繝ｳ繧ｰ繝��繧ｿ縺九ｉ蝓ｺ譛ｬ蝣ｴ繧剃ｽ懈�縺吶ｋ
    !

    !繝｢繧ｸ繝･繝ｼ繝ｫ隱ｭ縺ｿ霎ｼ縺ｿ
    use dc_types,  only: DP
    use gridset,   only: kmin, kmax,       &!驟榊�繧ｵ繧､繧ｺ (Z 譁ｹ蜷�)
      &                  nz                 !迚ｩ逅��伜沺縺ｮ螟ｧ縺阪＆ (Z 譁ｹ蜷�)
    use axesset,   only: r_Z, z_Z,         &!鬮伜ｺｦ
      &                  dz                 !譬ｼ蟄宣俣髫� (Z 譁ｹ蜷�)
    use constants, only: Grav,             &
      &                  GasRDry,          &
      &                  PressSfc, PressBasis, &
      &                  CpDry

    !證鈴ｻ吶�蝙句ｮ｣險遖∵ｭ｢
    implicit none

    !螟画焚螳夂ｾｩ
    real(DP), intent(out):: z_Press(kmin:kmax)           !蝨ｧ蜉�
    real(DP), intent(out):: z_Temp(kmin:kmax)            !貂ｩ蠎ｦ
    real(DP)             :: r_Press(kmin:kmax)           !蝨ｧ蜉�
    real(DP)             :: r_Temp(kmin:kmax)            !貂ｩ蠎ｦ
    real(DP)             :: z_PTemp(kmin:kmax)           !貂ｩ菴�
    real(DP)             :: r_PTemp(kmin:kmax)           !貂ｩ菴�
    real(DP)             :: z_Stab(kmin:kmax)            !螳牙ｮ壼ｺｦ
    real(DP)             :: Dens                         !蟇�ｺｦ
    real(DP)             :: z_DTempDZ(kmin:kmax)
    real(DP)             :: DTempDZ
    real(DP)             :: ratio
    integer              :: i, k, k1, k2, kmax2

    ! 蛻晄悄蛹�
    !
    z_PTemp = 0.0d0
    r_PTemp = 0.0d0
    z_Temp  = 0.0d0
    r_Temp  = 0.0d0
    z_Press = 0.0d0
    r_Press = 0.0d0
    z_Stab  = 0.0d0

    !!
    !! 邱壼ｽ｢陬懷ｮ�
    !!
    do k = kmin, kmax
      do i = 1, NumRec - 1
        if ( r_tmpAlt(i) == r_Z(k) ) then
          if (PTempCol > 0) r_PTemp(k) = r_tmpPTemp(i)
          if (TempCol  > 0) r_Temp(k)  = r_tmpTemp(i)
        elseif ( r_tmpAlt(i) < r_Z(k) .AND. r_Z(k) <= r_tmpAlt(i+1) ) then
          ratio = ( r_Z(k) - r_tmpAlt(i) ) / ( r_tmpAlt(i+1) - r_tmpAlt(i) )
          if (PTempCol > 0) r_PTemp(k) = r_tmpPTemp(i) + ( r_tmpPTemp(i+1) - r_tmpPTemp(i) ) * ratio
          if (TempCol  > 0) r_Temp(k)  = r_tmpTemp(i)  + ( r_tmpTemp(i+1)  - r_tmpTemp(i)  ) * ratio
        end if
      end do
    end do

    do k = kmin, kmax
      do i = 1, NumRec - 1
        if ( r_tmpAlt(i) == z_Z(k) ) then
          if (PTempCol > 0) z_PTemp(k) = r_tmpPTemp(i)
          if (TempCol  > 0) z_Temp(k)  = r_tmpTemp(i)
        elseif ( r_tmpAlt(i) < z_Z(k) .AND. z_Z(k) <= r_tmpAlt(i+1) ) then
          ratio = ( z_Z(k) - r_tmpAlt(i) ) / ( r_tmpAlt(i+1) - r_tmpAlt(i) )
          if (PTempCol > 0) z_PTemp(k) = r_tmpPTemp(i) + ( r_tmpPTemp(i+1) - r_tmpPTemp(i) ) * ratio
          if (TempCol  > 0) z_Temp(k)  = r_tmpTemp(i)  + ( r_tmpTemp(i+1)  - r_tmpTemp(i)  ) * ratio
        end if
      end do
    end do

!!    do k = kmin+1, kmax
!!      if (PTempCol > 0) z_PTemp(k) = ( r_PTemp(k-1) + r_PTemp(k) ) * 5.0d-1
!!      if (TempCol  > 0) z_Temp(k)  = ( r_Temp(k-1)  + r_Temp(k)  ) * 5.0d-1
!!    end do

    !! 繝��繧ｿ縺ｮ蜈･縺｣縺ｦ縺�ｋ驟榊�繧ｵ繧､繧ｺ繧貞ｾ励ｋ
    !!
    kmax2 = nz
    do k = 1, nz
      if ( max( r_Temp(k), r_PTemp(k) ) == 0.0d0 ) then
        kmax2 = k - 1
        exit
      end if
    end do

    !! 髱呎ｰｴ蝨ｧ蟷ｳ陦｡縺九ｉ蝓ｺ譛ｬ蝣ｴ縺ｮ蝨ｧ蜉帙ｒ菴懊ｊ逶ｴ縺�
    !! 鬮伜ｺｦ繝ｻ貂ｩ蠎ｦ => 豌怜悸
    !!
    if ( TempCol > 0 ) then
      r_Press(0) = PressSfc
      z_Press(1) = r_Press(0) * exp( - Grav * dz * 0.5d0 / GasRDry / r_Temp(0) )

      do k = 1, kmax2
        z_Press(k+1) = z_Press(k) * exp( - Grav * dz / GasRDry / r_Temp(k) )
      end do

      ! 蠢ｵ縺ｮ轤ｺ縺ｫ貂ｩ菴阪ｒ險育ｮ励＠縺ｦ縺翫￥
      !
      do k = 1, kmax2
        z_PTemp(k) = z_Temp(k) * (PressBasis / z_Press(k)) ** (GasRDry / CpDry)
      end do
    end if

    !! 髱呎ｰｴ蝨ｧ蟷ｳ陦｡縺九ｉ蝓ｺ譛ｬ蝣ｴ縺ｮ蝨ｧ蜉帙→貂ｩ蠎ｦ繧剃ｽ懊ｊ逶ｴ縺�
    !! 鬮伜ｺｦ繝ｻ貂ｩ菴� => 貂ｩ蠎ｦ繝ｻ豌怜悸
    !!
    if ( PTempCol > 0 ) then
      r_Press(0) = PressSfc
      r_Temp(0)  = ( PressBasis / r_Press(0) ) ** ( - GasRDry / CpDry ) * r_PTemp(0)

      do k = 1, kmax2
        Dens = r_Press(k-1) * ( (r_Press(k-1) / PressBasis )**(- GasRDry / CpDry)) / GasRDry / r_PTemp(k-1)
        z_Press(k) = r_Press(k-1) + ( - Dens * Grav ) * ( dz * 0.5d0 )
        z_Temp(k)  = ( PressBasis / z_Press(k) ) ** ( - GasRDry / CpDry ) * z_PTemp(k)

        Dens = z_Press(k) * ( (z_Press(k) / PressBasis )**(- GasRDry / CpDry)) / GasRDry / z_PTemp(k)
        r_Press(k) = z_Press(k) + ( - Dens * Grav ) * ( dz * 0.5d0 )
        r_Temp(k)  = ( PressBasis / r_Press(k) ) ** ( - GasRDry / CpDry ) * r_PTemp(k)
      end do
    end if

    !! 繧ｵ繧ｦ繝ｳ繝�ぅ繝ｳ繧ｰ繝輔ぃ繧､繝ｫ縺ｧ荳弱∴縺滄ｫ伜ｺｦ繧医ｊ荳顔ｩｺ縺ｮ蜃ｦ逅�
    !! VIRA 繝��繧ｿ繧剃ｽｿ縺��ｴ蜷�
    !!
    if ( kmax2 < nz ) then

      do k = kmax2+1, kmax

        if ( r_viraAlt(1) >= r_Z(k) ) then
          ratio = ( r_Z(k) - r_Z(kmax2) ) / ( r_viraAlt(1) - r_Z(kmax2) )
          r_Temp(k)  = r_Temp(kmax2)  + ( r_viraTemp(1)  - r_Temp(kmax2)  ) * ratio
        else
          do i = 1, NumVIRA-1
            if ( r_viraAlt(i) == r_Z(k) ) then
              r_Temp(k)  = r_viraTemp(i)
            elseif ( r_viraAlt(i) < r_Z(k) .AND. r_Z(k) <= r_viraAlt(i+1) ) then
              ratio = ( r_Z(k) - r_viraAlt(i) ) / ( r_viraAlt(i+1) - r_viraAlt(i) )
              r_Temp(k)  = r_viraTemp(i)  + ( r_viraTemp(i+1)  - r_viraTemp(i)  ) * ratio
            end if
          end do
        end if

        z_Temp(k)  = ( r_Temp(k-1) + r_Temp(k) ) * 5.0d-1
        z_Press(k) = z_Press(k-1) * exp( - Grav * dz / GasRDry / r_Temp(k-1) )
      end do
    end if

!!    do k = 1, nz
!!      write(*,*) k, z_Z(k), z_Temp(k)
!!      write(*,*) k, r_Z(k), r_Temp(k)
!!    end do

    !! 螳牙ｮ壼ｺｦ縺ｮ險育ｮ�
    !!
    do k = 1, nz
      z_Stab(k) = Grav / z_Temp(k)                  &
        &    * (                                    &
        &       + (r_Temp(k) - r_Temp(k-1)) / dz    &
        &       + Grav / CpDry                      &
        &      )
    end do

    !! 繧ｵ繧ｦ繝ｳ繝�ぅ繝ｳ繧ｰ繝��繧ｿ縺ｮ繝輔ぃ繧､繝ｫ蜃ｺ蜉�
    !!
    call CheckOutput

!!!
!!! 蝨冗阜髱｢縺悟ｭ伜惠縺吶ｋ蝣ｴ蜷医↓縺ｯ, 縺昴％縺ｧ縺ｮ貂ｩ蠎ｦ蜍ｾ驟阪ｒ邱ｩ繧�°縺ｫ縺吶ｋ.
!!!
    if ( AltTr > z_Z(1) .AND. AltTr < z_Z(nz) ) then

      ! 迴ｾ蝨ｨ縺ｮ貂ｩ蠎ｦ蛻�ｸ��蜍ｾ驟阪ｒ蜿悶ｋ
      !
      do k = 1, kmax
        z_DTempDZ(k) = (z_Temp(k) - z_Temp(k-1)) / dz
      end do

      ! 蟇ｾ豬∝恟逡碁擇繧医ｊ荳翫�謇ｱ縺�. 謖�ｮ壹＆繧後◆鬮伜ｺｦ縺ｮ貂ｩ蠎ｦ貂帷紫繧剃ｽｿ縺�ｶ壹￠繧�.
      !
      k1 = minloc( z_Z, 1, z_Z > AltTr )

      ! 貂ｩ蠎ｦ貂帷紫
      !
      DTempDZ = z_DTempDZ(k1-1)
      k2 = int( DelAlt / dz )    !謨ｴ謨ｰ蝙九∈螟画鋤

      do k = k1, kmax

        ! 貂ｩ蠎ｦ貅千紫繧偵ぞ繝ｭ縺ｫ霑代▼縺代ｋ.
        !
        DTempDZ = min( -1.0d-14, DTempDZ - DTempDZ / k2 * ( k - k1 ) )

        !蝓ｺ譛ｬ蝣ｴ縺ｮ貂ｩ蠎ｦ繧呈ｱｺ繧√ｋ
        z_Temp(k) = z_Temp(k-1) + DTempDZ * dz

        !蝨ｧ蜉帙ｒ髱呎ｰｴ蝨ｧ蟷ｳ陦｡縺九ｉ險育ｮ�
        z_Press(k) =                                      &
          &  z_Press(k-1) * ( ( z_Temp(k-1) / z_Temp(k) ) &
          &    ** (Grav / ( DTempDZ * GasRDry ) ) )

      end do

    end if

  contains

    function CalcCp(Temp)

      real(DP) :: Temp
      real(DP) :: CalcCp

      real(DP), parameter  :: CpA =  19.796d0
      real(DP), parameter  :: CpB =  7.344d-2
      real(DP), parameter  :: CpC = -5.602d-5
      real(DP), parameter  :: CpD =  1.715d-8

      CalcCp = ( CpA + CpB * Temp + CpC * Temp**2 + CpD * Temp**3 ) / 44.0d-3

    end function CalcCp

    subroutine CheckOutput

      use gtool_history,  only : HistoryCreate, HistoryClose, &
        &                        HistoryPut, HistoryAddVariable
      use fileset,        only : filetitle,      &!繝��繧ｿ縺ｮ陦ｨ鬘�
        &                        filesource,     &!繝��繧ｿ繧剃ｽ懈�縺吶ｋ謇矩��
        &                        FileInstitution  !譛邨ょ､画峩閠��邨�ｹ�
      use mpi_wrapper,    only : myrank

      real(DP) :: z_Stab2(kmin:kmax)            !螳牙ｮ壼ｺｦ
      real(DP) :: Cp

      ! CPU 縺ｮ繝ｩ繝ｳ繧ｯ縺� 0 縺ｮ蝣ｴ蜷医↓縺ｮ縺ｿ繝輔ぃ繧､繝ｫ蜃ｺ蜉�
      !
      if (myrank /= 0) return

      ! 蛻晄悄蛹�
      !
      z_Stab2 = 0.0d0

      !! 荳ｻ繝励Ο繧ｰ繝ｩ繝�縺ｮ險育ｮ励↓繧医▲縺ｦ, z_PTemp, z_Temp, z_Press 縺ｯ豎ゅ∪縺｣縺ｦ縺�ｋ.
      ! 繧ｵ繧ｦ繝ｳ繝�ぅ繝ｳ繧ｰ繝輔ぃ繧､繝ｫ縺ｧ貂ｩ蠎ｦ縺御ｸ弱∴繧峨ｌ縺溷�ｴ蜷医↓縺ｯ, 豈皮�縺ｮ貂ｩ蠎ｦ萓晏ｭ俶ｧ繧�
      ! 閠��縺励※螳牙ｮ壼ｺｦ繧定ｨ育ｮ励☆繧�.
      !
      if (TempCol > 0) then
        do k = 1, nz
          Cp = CalcCp( z_Temp(k) )

          z_Stab2(k) = Grav / z_Temp(k)                 &
            &    * (                                    &
            &       + (r_Temp(k) - r_Temp(k-1)) / dz    &
            &       + Grav / Cp                         &
            &      )
        end do
      end if

      !-------------------------------------------------------------
      ! 繝偵せ繝医Μ繝ｼ菴懈�
      !-------------------------------------------------------------
      call HistoryCreate(                              &
        & file = Outputfile,                           &
        & title = filetitle,                           &
        & source = filesource,                         &
        & institution = FileInstitution,               &
        & dims=(/'z'/),                                &
        & dimsizes=(/nz/),                             &
        & longnames=(/'Z-coordinate'/),                &
        & units=(/'m'/),                               &
        & xtypes=(/'float'/),                          &
        & history=rstat, quiet=.true. )

      call HistoryAddVariable(                         &
        & varname='Stab', dims=(/'z'/),                &
        & longname='Stability',                        &
        & units='1', xtype='double', history=rstat )

      call HistoryAddVariable(                         &
        & varname='Stab2', dims=(/'z'/),               &
        & longname='Stability with Cp(z)',             &
        & units='1', xtype='double', history=rstat )

      call HistoryAddVariable(                         &
        & varname='Temp', dims=(/'z'/),                &
        & longname='Temperature',                      &
        & units='1', xtype='double', history=rstat )

      call HistoryAddVariable(                         &
        & varname='PTemp', dims=(/'z'/),               &
        & longname='Temperature',                      &
        & units='1', xtype='double', history=rstat )

      call HistoryPut('z', z_Z(1:nz), rstat )

      call HistoryPut( 'Stab',   z_Stab(1:nz),  rstat )
      call HistoryPut( 'Stab2',  z_Stab2(1:nz), rstat )
      call HistoryPut( 'Temp',   z_Temp(1:nz),  rstat )
      call HistoryPut( 'PTemp' , z_PTemp(1:nz), rstat )

      call HistoryClose( rstat )

    end subroutine CheckOutput

  end subroutine Initialdata_sounding_basic

  subroutine initialdata_sounding_wind(pyz_VelX, xqz_VelY)

    !繝｢繧ｸ繝･繝ｼ繝ｫ隱ｭ縺ｿ霎ｼ縺ｿ
    use dc_types, only: DP
    use gridset,  only: imin, imax,       &!驟榊�繧ｵ繧､繧ｺ (X 譁ｹ蜷�)
      &                 jmin, jmax,       &!驟榊�繧ｵ繧､繧ｺ (Y 譁ｹ蜷�)
      &                 kmin, kmax         !驟榊�繧ｵ繧､繧ｺ (Z 譁ｹ蜷�)
    use axesset,  only: r_Z                !鬮伜ｺｦ

    implicit none

    real(DP), intent(out) :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xqz_VelY(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)              :: pyr_VelX(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)              :: xqr_VelY(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)              :: ratio
    integer               :: i, k

    !蛻晄悄蛹�
    pyz_VelX = 0.0d0
    pyr_VelX = 0.0d0
    xqz_VelY = 0.0d0
    xqr_VelY = 0.0d0

    ! 邱壼ｽ｢陬懷ｮ�
    !
    do k = kmin, kmax
      do i = 1, NumRec
        if ( r_tmpAlt(i) <= r_Z(k) .AND. r_Z(k) < r_tmpAlt(i+1) ) then
          ratio = ( r_Z(k) - r_tmpAlt(i) ) / ( r_tmpAlt(i+1) - r_tmpAlt(i) )
          pyr_VelX(:,:,k) = r_tmpVelX(i) + ( r_tmpVelX(i+1) - r_tmpVelX(i) ) * ratio
          xqr_VelY(:,:,k) = r_tmpVelY(i) + ( r_tmpVelY(i+1) - r_tmpVelY(i) ) * ratio
        end if
      end do
    end do

    ! r => z 縺ｮ螟画鋤
    !
    do k = kmin+1, kmax
      pyz_VelX(:,:,k) = ( pyr_VelX(:,:,k-1) + pyr_VelX(:,:,k) ) * 5.0d-1
      xqz_VelY(:,:,k) = ( xqr_VelY(:,:,k-1) + xqr_VelY(:,:,k) ) * 5.0d-1
    end do

  end subroutine initialdata_sounding_wind

  subroutine devidecsv( buf, eachcol, maxch, num )
    !
    ! csv 蠖｢蠑上�繝輔ぃ繧､繝ｫ縺九ｉ蠢�ｦ√↑諠��ｱ繧定ｪｭ縺ｿ蜃ｺ縺�
    !

    !證鈴ｻ吶�蝙句ｮ｣險遖∵ｭ｢
    implicit none

    !螟画焚螳夂ｾｩ
    integer          :: maxch, num
    character(len=*) :: buf, eachcol(maxch)
    integer          :: i, j, prev, now, now1, ll
    logical          :: quote

    quote = .FALSE.

    ll = len_trim(buf)

    prev=0; now=1; i=1
    do j=1, ll
      if( buf(j:j) == '"' ) then
        quote = .NOT. quote
      end if
      if( quote .eqv. .FALSE. ) then
        if(buf(j:j) == ',') then
          now = j
          now1 = now -1
          if( now1 < prev+1 ) then
            eachcol(i) = ''
          else
            eachcol(i) = buf(prev+1:now1)
          end if
          i=i+1
          if( i > maxch ) then
            write(0,*) 'maxch is too small'
            stop
          end if
          prev=now
        end if
      end if
    end do

    if( prev < ll ) then
      eachcol(i) =  buf(prev+1:ll)
      num = i
    else
      num=i-1
    end if

  end subroutine devidecsv

end module initialdata_sounding