!=================================================================
!  鉛直 1 次元火星大気放射対流平衡モデル:
!  鉛直構造計算
!
!  1999/10/17 小高正嗣
!=================================================================

      SUBROUTINE VPROF
     I   ( ZGRID, PSFC,
     O     PRES0, DENS0, TEMP0, TPOT0, PPAI0,
     O     QVAP0, RVAP0,
     C     GRAV , RDRY , CPDRY, PSD,  RMVAP,
     G     FDZ  , FDZM , DX   ,
     D     NZ   , NZB                        )

!-----------------------------------------------------------------
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------

      DIMENSION
     \   PRES0(-2:NZB), TEMP0(-2:NZB), TPOT0(-2:NZB),
     \   PPAI0(-2:NZB), DENS0(-2:NZB), ZGRID(-2:NZB),
     \   RVAP0(-2:NZB), QVAP0(-2:NZB), FDZ(-2:NZB), FDZM(-2:NZB)

      dimension PPLOG(-2:NZB)

!-----------------------------------------------------------------
! 静水圧の式の積分

      PPLOG(0) = log(PSFC) -
     \     GRAV / RDRY / TEMP0(0) * FDZ(0) * DX * 0.5D0
      PRES0(0) = exp( PPLOG(0) )
      DENS0(0) = PRES0(0) / RDRY / TEMP0(0)

      do 101 IZ = 1, NZ
         tmid = ( TEMP0( IZ ) + TEMP0( IZ-1 ) ) * 0.5D0

         PPLOG( IZ ) = PPLOG( IZ - 1 )
     \        - GRAV / RDRY / tmid * FDZM( IZ ) * DX
         PRES0( IZ ) = exp( PPLOG(IZ) )
         DENS0( IZ ) = PRES0( IZ ) / ( RDRY * TEMP0(IZ) )
 101  continue

c      write(*,*) 'IZ   ','TEMP0     ','TPOT0      ','PRES0      ',
c     \        'DENS0      ','PPAI0'
c      do iz = 0, 30
c      write(*,*) iz, TEMP0(iz), TPOT0(iz), PRES0(iz),
c     \        DENS0(iz),PPAI0(iz)
c      end do

!-----------------------------------------------------------------
! 温位, 湿潤変数の計算

      do 201 IZ = 0, NZ

         PPAI0( IZ ) = ( PRES0(IZ) / PSD ) ** ( RDRY / CPDRY )
         TPOT0( IZ ) = TEMP0( IZ ) / PPAI0( IZ )

         RVAP0( IZ ) = 0.0D0
         QVAP0(IZ) = exp ( 0.007D0 * ( 1.0D0 - psfc / PRES0(iz) ) )
c         QVAP0(IZ) = RVAP0(IZ) * EVSAT( TEMP0(IZ) ) / PRES0(IZ)
c     \                         * RMVAP

 201  continue

!-----------------------------------------------------------------
! 境界条件の計算

c      ZGRID(   -1 ) = ZGRID( 0 )
c      ZGRID(   -2 ) = ZGRID( 1 )
c      ZGRID( NZ   ) = ZGRID( NZ-1 )
c      ZGRID( NZ+1 ) = ZGRID( NZ-2 )
c      ZGRID( NZ+2 ) = ZGRID( NZ-3 )
C
      PRES0(   -1 ) = PRES0(    0 )
      PRES0(   -2 ) = PRES0(    1 )
c      PRES0( NZ   ) = PRES0( NZ-1 )
c      PRES0( NZ+1 ) = PRES0( NZ-2 )
c      PRES0( NZ+2 ) = PRES0( NZ-3 )
C
      PPAI0(   -1 ) = PPAI0(    0 )
      PPAI0(   -2 ) = PPAI0(    1 )
c      PPAI0( NZ   ) = PPAI0( NZ-1 )
c      PPAI0( NZ+1 ) = PPAI0( NZ-2 )
c      PPAI0( NZ+2 ) = PPAI0( NZ-3 )
C
      DENS0(   -1 ) = DENS0(    0 )
      DENS0(   -2 ) = DENS0(    1 )
c      DENS0( NZ   ) = DENS0( NZ-1 )
c      DENS0( NZ+1 ) = DENS0( NZ-2 )
c      DENS0( NZ+2 ) = DENS0( NZ-3 )
C
      TEMP0(   -1 ) = TEMP0(    0 )
      TEMP0(   -2 ) = TEMP0(    1 )
c      TEMP0( NZ   ) = TEMP0( NZ-1 )
c      TEMP0( NZ+1 ) = TEMP0( NZ-2 )
c      TEMP0( NZ+2 ) = TEMP0( NZ-3 )
C
      TPOT0(   -1 ) = TPOT0(    0 )
      TPOT0(   -2 ) = TPOT0(    1 )
c      TPOT0( NZ   ) = TPOT0( NZ-1 )
c      TPOT0( NZ+1 ) = TPOT0( NZ-2 )
c      TPOT0( NZ+2 ) = TPOT0( NZ-3 )
C
      RVAP0(   -1 ) = RVAP0(    0 )
      RVAP0(   -2 ) = RVAP0(    1 )
c      RVAP0( NZ   ) = RVAP0( NZ-1 )
c      RVAP0( NZ+1 ) = RVAP0( NZ-2 )
c      RVAP0( NZ+2 ) = RVAP0( NZ-3 )
C
      QVAP0(   -1 ) = QVAP0(    0 )
      QVAP0(   -2 ) = QVAP0(    1 )
c      QVAP0( NZ   ) = QVAP0( NZ-1 )
c      QVAP0( NZ+1 ) = QVAP0( NZ-2 )
c      QVAP0( NZ+2 ) = QVAP0( NZ-3 )
C
C---------------------------------------------------------------------
c      WRITE(6,*) '**DENS0 OBTAINED IN SUBROUTINE VBASIC'
c      WRITE(6,*) DENS0
C
C
      RETURN
      E N D