! ã‚µã‚¦ãƒ³ãƒ‡ã‚£ãƒ³ã‚°ãƒ•ã‚¡ã‚¤ãƒ«ã‹ã‚‰åˆæœŸå ´ã‚’æ±ºã‚ã‚‹ãŸã‚ã®ã‚µãƒ–ãƒ«ãƒ¼ãƒãƒ³
!
! 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_sounding2
  !
  ! ã‚µã‚¦ãƒ³ãƒ‡ã‚£ãƒ³ã‚°ãƒ•ã‚¡ã‚¤ãƒ«ã‹ã‚‰åˆæœŸå ´ã‚’æ±ºã‚ã‚‹ãŸã‚ã®ã‚µãƒ–ãƒ«ãƒ¼ãƒãƒ³
  ! ã‚µã‚¦ãƒ³ãƒ‡ã‚£ãƒ³ã‚°ãƒ•ã‚¡ã‚¤ãƒ«ã¯ãƒ†ã‚ã‚¹ãƒˆãƒ•ã‚¡ã‚¤ãƒ«ã§æ›¸ã‹ã‚Œã¦ã„ã‚‹.
  ! å°†æ¥çš„ã«ã¯ netCDF ã«å¤‰æ›´ã™ã‚‹äºˆå®š.
  !

  !ãƒ¢ã‚¸ãƒ¥ãƒ¼ãƒ«èªã¿è¾¼ã¿
  use dc_types, only: DP

  !æš—é»™ã®åž‹å®£è¨€ç¦æ¢
  implicit none

  !å…±é€šå¤‰æ•°
  real(DP), save, private :: r_tmpAlt(10000)
  real(DP), save, private :: r_tmpTemp(10000)
  real(DP), save, private :: r_tmpPress(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 :: r_tmpQrad(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

  !å…¬é–‹å¤‰æ•°
  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             :: AltCol = 0
    integer             :: TempCol = 0
    integer             :: PTempCol = 0
    integer             :: PressCol = 0
    integer             :: VelXCol = 0
    integer             :: VelYCol = 0
    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, PressCol, 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, max( TempCol, PressCol ) )
    r_tmpAlt  = 0.0d0; r_tmpTemp = 0.0d0; r_tmpPress = 0.0d0;
    r_tmpVelX = 0.0d0; r_tmpVelY = 0.0d0

    ! ãƒ•ã‚¡ã‚¤ãƒ«ã®ã‚ªãƒ¼ãƒ—ãƒ³
    !
    open (17, file=SoundingFile, status='old')

    ! ãƒ•ã‚¡ã‚¤ãƒ«å‘¼ã³å‡ºã—
    !
    do while ( io == 0 )
      ! 1 è¡Œåˆ†èªã¿å‡ºã—
      !
      read (17, '(a)', IOSTAT=io) 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(TempCol)(1:len_trim(eachcol(PTempCol))), *)  r_tmpPTemp(NumRec)
        if (PressCol > 0) read( eachcol(PressCol)(1:len_trim(eachcol(PressCol))), *) r_tmpPress(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

    ! ãƒ•ã‚¡ã‚¤ãƒ«ã®ã‚¯ãƒãƒ¼ã‚º
    !
    close (17)

  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

    !æš—é»™ã®åž‹å®£è¨€ç¦æ¢
    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_DTempDZ(kmin:kmax)
    real(DP)             :: DTempDZ
    real(DP)             :: ratio
    integer              :: i, k, k1, k2

    real(DP)             :: r_PTemp(kmin:kmax)
    real(DP)             :: z_PTemp(kmin:kmax)

    ! åˆæœŸåŒ–
    !
    z_Temp  = 0.0d0
    r_Temp  = 0.0d0
    z_PTemp = 0.0d0
    r_PTemp = 0.0d0
    z_Press = 0.0d0
    r_Press = 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) )
          r_PTemp(k)  = r_tmpPTemp(i)  + ( r_tmpPTemp(i+1) - r_tmpPTemp(i) ) * ratio
        end if
      end do
    end do

    ! r => z ã®å¤‰æ›
    !
    do k = kmin+1, kmax
      z_PTemp(k)  = ( r_PTemp(k-1)  + r_PTemp(k)  ) / 2.0d0
    end do

    ! é™æ°´åœ§ã‹ã‚‰åœ§åŠ›ã‚’æ±ºã‚ã‚‹
    !
    do k = 1, nz
      z_Press(k+1)
    end do

    ! æ¸©ä½ã‹ã‚‰æ¸©åº¦ã‚’ä½œã‚‹
    !
    z_Temp = xyz_PTemp / ( (PressBasis / xyz_Press) ** (GasRDry / CpDry) )

    !

!!!
!!! åœç•Œé¢ãŒå˜åœ¨ã™ã‚‹å ´åˆã«ã¯, ãã“ã§ã®æ¸©åº¦å‹¾é…ã‚’ç·©ã‚„ã‹ã«ã™ã‚‹.
!!!
    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

  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, z_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_sounding2