!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE INIBASE                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE inibase(nx,ny,nz,                                            & 1,6
           ubar,vbar,ptbar,pbar,ptbari,pbari,                           &
           rhostr,rhostri,qvbar,                                        &
           x,y,z,zp,j3, rhobar, zs,zuv,pibar)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Initialize the base-state variables.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/17/1991.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/03/92 (M. Xue)
!  Further documentation.
!
!  10/05/92 (M. Xue)
!  zero gradient in ptbar at the bottom and top boundaries are
!  no longer imposed.
!
!  2/10/93 (K. Droegemeier)
!  Cleaned up documentation.
!
!  9/10/94 (Weygandt & Y. Lu)
!  Cleaned up documentation.
!
!  10/31/94 (E. Adlerman)
!  Added option for dewpoint input on sounding and for
!  inibasopt = 5, i.e. Weisman&Klemp analytical thermodynamic
!  profile.
!
!  12/22/94 (Yuhe Liu)
!  Added option for base state wind input for inibasopt != 1 and
!  viniopt = 1 and 2.
!
!  5/10/95 (M. Xue)
!  Corrected a bug with the pibar calculation for inibasopt=4 case.
!
!  5/10/95 (M. Xue)
!  Changed lower boundary condition for ptbar and qvbar to constant
!  gradient extrapolation. rhobar, pbar and pibar at k=1 are obtained
!  from hydrostatic equation.
!
!  7/6/95 (M. Xue)
!  Zero gradient condition restored for ptbar and qvbar at the
!  bottom boundary. The extrapolated condition was introducing
!  undesirable turbulent heat and moisture fluxes at the ground.
!
!  2/11/1997 (M. Xue and D. Weber)
!  Corrected a problem with base-state initialization option 4,
!  the case of constant static stability.
!
!  8/4/1997 (M. Xue and D. Weber)
!  Fixed a problem with pibar calculation for option 4, introduced
!  on 2/11/1997.
!
!  5/1/1998 (G. Bassett and D. Weber)
!  Moved sounding output to subroutine writesnd.
!
!  9/15/1998 (D. Weber)
!  Added ptbari, pbari, and rhostri (inverse quantities).
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!  OUTPUT:
!
!    ubar     Base state zonal velocity component (m/s)
!    vbar     Base state meridional velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    ptbari   Inverse base state potential temperature (K)
!    pbari    Inverse base state pressure (Pascal)
!    rhostr   Base state density rhobar times j3 (kg/m**3)
!    rhostri  Inverse base state density rhobar times j3 (kg/m**3)
!    qvbar    Base state water vapor specific humidity (kg/kg)
!    pibar    Base state
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!    zp       Vertical coordinate of grid points in physical space(m)
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!
!  WORK ARRAYS:
!
!    rhobar   Temporary work array.
!    zs       Temporary work array.
!    zuv      Temporary work array.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'mp.inc'          ! Message passing parameters.


  INTEGER :: nx,ny,nz          ! The number of grid points in 3
                               ! directions

  REAL :: ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
  REAL :: ptbari(nx,ny,nz)     ! Inverse base state pot. temperature (K)
  REAL :: pbari (nx,ny,nz)     ! Inverse base state pressure (Pascal).
  REAL :: rhostr(nx,ny,nz)     ! Base state density rhobar times j3.
  REAL :: rhostri(nx,ny,nz)    ! Inverse base state density rhobar times j3.
  REAL :: qvbar (nx,ny,nz)     ! Base state water specific humidity
                               ! (kg/kg).
  REAL :: pibar (nx,ny,nz)     ! Base state Exner function
  REAL :: x     (nx)           ! x-coord. of the physical and compu-
                               ! tational grid. Defined at u-point.
  REAL :: y     (ny)           ! y-coord. of the physical and compu-
                               ! tational grid. Defined at v-point.
  REAL :: z     (nz)           ! z-coord. of the computational grid.
                               ! Defined at w-point on the staggered
                               ! grid.
  REAL :: zp    (nx,ny,nz)     ! Physical height coordinate defined at
                               ! w-point of staggered grid.
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian
                               ! as d(zp)/dz.

  REAL :: rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
  REAL :: zs    (nx,ny,nz)     ! Physical coordinate height at scalar
                               ! point (m)
  REAL :: zuv   (nx,ny,nz)     ! Physical coordinate height at u or v
                               ! point (m)

  INTEGER :: lvlprof
  REAL :: depthp

  PARAMETER( lvlprof=601, depthp = 3.0E4 )

!
!-----------------------------------------------------------------------
!
!  lvlprof:
!
!  The ARPS interpolates unevenly-spaced data from a vertical
!  sounding to evenly-spaced data for use in defining the base
!  state atmosphere.  In this process, an intermediate sounding is
!  generated at evenly-spaced altitudes, with the accuracy of the
!  associated interpolation controlled by the parameter lvlprof.
!  The larger lvlprof, the more accurate the interpolation
!  (we recommend using lvlprof=200 for a model run with about 50
!  points in the vertical direction).  Using the intermediate
!  sounding, the ARPS then generates a base state model sounding
!  for the particular vertical grid resolution chosen (i.e.,
!  the number of points in the vertical, nz, and the vertical grid
!  spacing, dz).
!
!  depthp:
!
!  The depth of atmosphere over which the interpolated profiles
!  will be defined.  depthp should be greater than or equal to
!  (nz-3)*dz, i.e., larger than the physical depth of the model
!  domain.  Otherwise, the code will extrapolate for gridpoints
!  outside the domain, leading to possible inconsistencies.  At all
!  costs, any such extrapolation should be avoided.
!
!-----------------------------------------------------------------------
!

  REAL :: uprof (lvlprof)      ! Temporary work array
  REAL :: vprof (lvlprof)      ! Temporary work array
  REAL :: ptprof(lvlprof)      ! Temporary work array
  REAL :: pprof (lvlprof)      ! Temporary work array
  REAL :: qvprof(lvlprof)      ! Temporary work array
  REAL :: zprof (lvlprof)      ! Temporary work array

  REAL :: temp1(lvlprof)       ! Temporary work array
  REAL :: temp2(lvlprof)       ! Temporary work array
  REAL :: temp3(lvlprof)       ! Temporary work array
  REAL :: temp4(lvlprof)       ! Temporary work array
  REAL :: temp5(lvlprof)       ! Temporary work array
  REAL :: temp6(lvlprof)       ! Temporary work array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  REAL :: factor,deltaz,pt0,t0
  INTEGER :: i,j,k,kdata, istat, nunit, onvf, itrnmin, jtrnmin
  REAL :: p0inv,cpdrd,ptbar0, nstatsq, ternmin, rho00, tem
  LOGICAL :: wk  ! Logical control for Weisman & Klemp sounding.
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  p0inv=1.0/p0
  cpdrd=1.0/rddcp

  IF( inibasopt == 5 ) THEN
    wk = .true.
  ELSE
    wk = .false.
  END IF

  onvf = 0
  CALL avgz(zp, onvf,                                                   &
            nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, zs)

  IF( inibasopt == 1 .OR. inibasopt == 5 ) THEN
!
!-----------------------------------------------------------------------
!
!  Initialize the base state from discrete sounding data.  Note that
!  The desired format of the sounding data is described inside
!  subroutine SOUNDG.
!
!-----------------------------------------------------------------------
!
    deltaz = depthp/(lvlprof-1)

    IF (myproc ==0) &
       WRITE(6,'(/'' THE INPUT SOUNDING DATA IS SAMPLED AT A '',           &
    &         f7.2,'' meter interval.''/)') deltaz
    
    DO k=1,lvlprof
      zprof(k) = (k-1)*deltaz
    END DO
!
!-----------------------------------------------------------------------
!
!  Obtain, from the input sounding, the base state vertical profiles
!  of the two horizontal wind components, potential temperature,
!  pressure, and water vapor specific humidity at vertical levels
!  defined by zprof.
!
!-----------------------------------------------------------------------
!
    CALL zprofil(zprof,lvlprof,                                         &
                 uprof,vprof,ptprof,pprof,qvprof,                       &
                 temp1,temp2,temp3,temp4,temp5,temp6,wk)
!
!
!-----------------------------------------------------------------------
!
!  Linearly interpolate the base state variables from the
!  intermediate evenly-spaced profile to the model grid.
!
!-----------------------------------------------------------------------
!

!-----------------------------------------------------------------------
!
!  Interpolation of scalars.
!
!-----------------------------------------------------------------------

    DO k=2,nz-2
      DO j=1,ny-1
        DO i=1,nx-1

          kdata = INT( zs(i,j,k)/deltaz ) + 1
          kdata = MIN( lvlprof-1, MAX(1,kdata) )
          factor = (zs(i,j,k)-zprof(kdata))/deltaz

          ptbar(i,j,k)=ptprof(kdata)+                                   &
                       factor*(ptprof(kdata+1)-ptprof(kdata))

          pbar (i,j,k)=pprof(kdata)+                                    &
                       factor*(pprof(kdata+1)-pprof(kdata))

          qvbar(i,j,k)=qvprof(kdata)+                                   &
                       factor*(qvprof(kdata+1)-qvprof(kdata))

          pibar(i,j,k) = (pbar(i,j,k)*p0inv)**rddcp
          rhobar(i,j,k)=pbar(i,j,k)/(rd*ptbar(i,j,k)*pibar(i,j,k))
        END DO
      END DO
    END DO

!-----------------------------------------------------------------------
!
!  Set top and bottom boundary conditions
!
!  Zero gradient is imposed on ptbar and qvbar at the top and bottom
!  boundaries. pt and qv above the top and below the bottom boundary
!  are used only in the calculations of turbulent fluxes through
!  the boundary. Zero gradient conditions ensure zero heat and
!  misture fluxes through the boundary. When such fluxes need to be
!  included, they are handled by the surface physics parameterization.
!
!  Hydrostatic relation is used to obtain pbar and rhobar above the
!  top and below the bottom boundary.
!
!-----------------------------------------------------------------------

    DO j=1,ny-1
      DO i=1,nx-1
        ptbar(i,j,1)=ptbar(i,j,2)
        qvbar(i,j,1)=qvbar(i,j,2)
      END DO
    END DO

    DO i=1,nx-1
      DO j=1,ny-1
        pibar(i,j,1)=pibar(i,j,2)+g*(zs(i,j,2)-zs(i,j,1))               &
                    /(0.5*cp*(ptbar(i,j,2)+ptbar(i,j,1)) )
        pbar (i,j,1) = p0 * pibar(i,j,1)**cpdrd
        rhobar(i,j,1)= pbar(i,j,1)/(rd*ptbar(i,j,1)*pibar(i,j,1))
      END DO
    END DO

    DO j=1,ny-1
      DO i=1,nx-1
        ptbar (i,j,nz-1)=ptbar(i,j,nz-2)
        qvbar (i,j,nz-1)=qvbar(i,j,nz-2)
      END DO
    END DO

    DO j=1,ny-1
      DO i=1,nx-1
        pibar(i,j,nz-1)=pibar(i,j,nz-2)-g*(zs(i,j,nz-1)-zs(i,j,nz-2))   &
                       /(0.5*cp*(ptbar(i,j,nz-2)+ptbar(i,j,nz-1)) )
        pbar (i,j,nz-1) = p0 * pibar(i,j,nz-1)**cpdrd
        rhobar(i,j,nz-1)= pbar(i,j,nz-1)/                               &
                         (rd*ptbar(i,j,nz-1)*pibar(i,j,nz-1))
      END DO
    END DO

!-----------------------------------------------------------------------
!
!  Interpolate u profile to model grid ubar
!  and subtract the domain translation speed when the option is on.
!
!-----------------------------------------------------------------------

    CALL avgsu(zs,nx,ny,nz, 1,ny-1, 1,nz-1, zuv, ubar) ! ubar used as temp

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx

          kdata = INT( zuv(i,j,k)/deltaz ) + 1
          kdata = MIN( lvlprof-1, MAX(1,kdata) )
          factor = (zuv(i,j,k)-zprof(kdata))/deltaz

          ubar(i,j,k)=uprof(kdata)+                                     &
                      factor*(uprof(kdata+1)-uprof(kdata))

        END DO
      END DO
    END DO

    IF(grdtrns == 1) THEN

      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx
            ubar(i,j,k) = ubar(i,j,k) - umove
          END DO
        END DO
      END DO

    END IF

!-----------------------------------------------------------------------
!
!  Interpolate v profile to model grid ubar
!  and subtract the domain translation speed when the option is on.
!
!-----------------------------------------------------------------------

!call test_dump (zs,"XXXbavgsv_zs",nx,ny,nz,0,1)
    CALL avgsv(zs,nx,ny,nz, 1,nx-1, 1,nz-1, zuv, vbar) ! vbar used as temp

    DO k=1,nz-1
      DO j=1,ny
        DO i=1,nx-1

          kdata = INT( zuv(i,j,k)/deltaz ) + 1
          kdata = MIN( lvlprof-1, MAX(1,kdata) )
          factor = (zuv(i,j,k)-zprof(kdata))/deltaz

          vbar(i,j,k)=vprof(kdata)+                                     &
                      factor*(vprof(kdata+1)-vprof(kdata))

        END DO
      END DO
    END DO

    IF(grdtrns == 1) THEN

      DO k=1,nz-1
        DO j=1,ny
          DO i=1,nx-1
            vbar(i,j,k) = vbar(i,j,k) - vmove
          END DO
        END DO
      END DO

    END IF


    DO j=1,ny-1
      DO i=1,nx
        ubar(i,j,1   )=ubar(i,j,2   )
        ubar(i,j,nz-1)=ubar(i,j,nz-2)
      END DO
    END DO

    DO j=1,ny
      DO i=1,nx-1
        vbar(i,j,   1)=vbar(i,j,   2)
        vbar(i,j,nz-1)=vbar(i,j,nz-2)
      END DO
    END DO

  ELSE IF( inibasopt == 2 ) THEN
!
!-----------------------------------------------------------------------
!
!  Specify a dry isentropic base state, using analytical formulations.
!  Hydrostatic relation is used.
!
!  This option can be switched on when needed.
!
!-----------------------------------------------------------------------
!
    pt0 = 300.0

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          pibar(i,j,k) = 1.0-g/(cp*pt0)*zs(i,j,k)
          pbar(i,j,k) = p0 * pibar(i,j,k)**cpdrd
          ptbar(i,j,k)=pt0
          rhobar(i,j,k)=pbar(i,j,k)/(rd*pt0*pibar(i,j,k))
          qvbar(i,j,k)=0.0
        END DO
      END DO
    END DO

  ELSE IF( inibasopt == 3) THEN
!
!-----------------------------------------------------------------------
!
!  Specify a dry isothermal base state, using analytical formulations.
!  Hydrostatic relation is used.
!
!  This option can be switched on when needed.
!
!-----------------------------------------------------------------------
!
    t0 = 250.0
    rho00 = p0/(rd*t0)

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          pbar(i,j,k) = p0 * EXP(-g*zs(i,j,k)/(rd*t0))
          pibar(i,j,k)=(pbar(i,j,k)*p0inv)**rddcp
          ptbar(i,j,k)=t0/pibar(i,j,k)
          rhobar(i,j,k)=pbar(i,j,k)/(rd*t0)
          IF(bsnesq == 1) rhobar(i,j,k)=rho00
          qvbar(i,j,k)=0.0
        END DO
      END DO
    END DO

  ELSE IF(  inibasopt == 4 ) THEN ! Case of constant static stability

    nstatsq = 0.0001             ! constant static stability (1/s**2)
    ptbar0 = 300.0               ! PTBAR at z=0.0

    rho00 = p0/(rd*ptbar0)

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          ptbar(i,j,k) = ptbar0*EXP(nstatsq*zs(i,j,k)/g )
        END DO
      END DO
    END DO

    tem = g*g/(nstatsq*cp*ptbar0)

    DO k=1,nz-1
      DO i=1,nx-1
        DO j=1,ny-1
          pibar(i,j,k)=1.0+tem*(EXP(-nstatsq/g*zs(i,j,k))-1.0)
        END DO
      END DO
    END DO

    DO k=1,nz-1
      DO i=1,nx-1
        DO j=1,ny-1
          pbar(i,j,k) =  p0*( pibar(i,j,k)**cpdrd )
          qvbar(i,j,k)=0.0
          rhobar(i,j,k)=pbar(i,j,k)/(rd*ptbar(i,j,k)*pibar(i,j,k))
          IF(bsnesq == 1) rhobar(i,j,k)=rho00
        END DO
      END DO
    END DO

  ELSE IF( inibasopt == 6 ) THEN

!-----------------------------------------------------------------------
!
!  Incompressible constant density and constant potential temperature
!  flow. Hydrostatic.
!
!-----------------------------------------------------------------------

    pt0 = 300.0

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          ptbar(i,j,k)=pt0
          rhobar(i,j,k)=p0/(rd*pt0*1.0)
          pbar(i,j,k)=p0 - rhobar(i,j,k)*g*zs(i,j,k)
          pibar(i,j,k) = (pbar(i,j,k)*p0inv)**rddcp
          qvbar(i,j,k)=0.0
        END DO
      END DO
    END DO

  ELSE

    IF (myproc ==0) WRITE(6,'(1x,a,i3,/1x,a)') &
        'Base state initialization option inibasopt=',inibasopt,        &
        'not implemented or available. Job stopped in INIBASE.'

    CALL arpsstop('arpsstop called from INIBASE with inibasopt',1)

  END IF

  IF ( inibasopt /= 1 ) THEN
!
!-----------------------------------------------------------------------
!
!  For all cases except when external sounding (inibasopt=1) is used,
!  set the base state winds according to the viniopt option.
!
!-----------------------------------------------------------------------
!

    IF ( viniopt == 1 ) THEN
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx
            ubar(i,j,k) = ubar0
          END DO
        END DO
      END DO

      DO k=1,nz-1
        DO j=1,ny
          DO i=1,nx-1
            vbar(i,j,k) = vbar0
          END DO
        END DO
      END DO

    ELSE IF ( viniopt == 2 ) THEN  ! User specifed wind profile
!
!-----------------------------------------------------------------------
!
!  When viniopt = 2, the profiles of ubar and vbar are supposed to be
!  specified by the user, by editing the following block of program.
!
!  As a default, ubar and vbar are set to zero.
!
!-----------------------------------------------------------------------
!
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx
            ubar(i,j,k) = 0.0
          END DO
        END DO
      END DO

      DO k=1,nz-1
        DO j=1,ny
          DO i=1,nx-1
            vbar(i,j,k) = 0.0
          END DO
        END DO
      END DO

    ELSE

     IF (myproc ==0) WRITE (6, '(1x,a,i4,a)')  &
            'Unexpected option for viniopt: ',viniopt,                  &
            ', Model did not stopped in subroutine INIBASE.'

    END IF

  END IF
!
!-----------------------------------------------------------------------
!
!  Calculate rhostr etc.
!
!-----------------------------------------------------------------------
!
  DO k= 1,nz-1
    DO j= 1,ny-1
      DO i= 1,nx-1
        ptbari(i,j,k)=1.0/ptbar(i,j,k)
        pbari(i,j,k) =1.0/pbar(i,j,k)
        rhostr(i,j,k)=ABS(j3(i,j,k))*rhobar(i,j,k)
        rhostri(i,j,k)=1.0/rhostr(i,j,k)
      END DO
    END DO
  END DO

  CALL writesnd(nx,ny,nz,ubar,vbar,ptbar,pbar,qvbar,zp, rhobar, zs)

  RETURN
END SUBROUTINE inibase

!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE ZPROFIL                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE zprofil(zprof,lvlprof,                                       & 1,2
           uprof,vprof,ptprof,pprof,qvprof,                             &
           zsnd,temp1,temp2,temp3,temp4,temp5,wk)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Compute the initial sounding profile at regularly-spaced grid
!  levels.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/17/1991.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/03/92 (M. Xue)
!  Further documentation.
!
!  2/10/93 (K. Droegemeier)
!  Cleaned up documentation.
!
!  9/10/94 (Weygandt & Y. Lu)
!  Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    zprof     1-D array defining the regularly-spaced vertical grid
!              levels at which base state profiles will be sampled(m)
!    lvlprof   The number of grid levels in zprof.
!    WK        Logical flag for Weisman&Klemp analytical sounding
!
!  OUTPUT:
!
!    uprof     Profile of u-velocity defined at levels given
!              in zprof (m)
!    vprof     Profile of v-velocity defined at levels given
!              in zprof (m)
!    ptprof    Profile of potential temperature defined at levels
!              given in zprof (K)
!    pprof     Profile of pressure defined at levels given in zprof
!              (Pascal)
!    qvprof    Profile of water vapor specific humidity defined at
!              levels given in zprof (kg/kg)
!
!  WORK ARRAYS:
!
!    zsnd     Temporary work array to be used to store the actual
!             height of the input sounding data.
!    temp1    Temporary work array
!    temp2    Temporary work array
!    temp3    Temporary work array
!    temp4    Temporary work array
!    temp5    Temporary work array
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: lvlprof          ! The number of grid levels in zprof

  REAL :: zprof (lvlprof)     ! 1-D array defining the regularly-spaced
                              ! vertical grid levels at which base
                              ! state profiles will be sampled (m)

  REAL :: uprof (lvlprof)     ! Profile of u-velocity defined at
                              ! levels given in zprof (m)
  REAL :: vprof (lvlprof)     ! Profile of v-velocity defined at
                              ! levels given in zprof (m)
  REAL :: ptprof(lvlprof)     ! Profile of potential temperature
                              ! defined at levels given in zprof (K)
  REAL :: pprof(lvlprof)      ! Profile of pressure defined at
                              ! levels given in zprof (Pascal)
  REAL :: qvprof(lvlprof)     ! Profile of water vapor specific
                              ! humidity defined at levels given in
                              ! zprof (kg/kg)

  REAL :: zsnd (lvlprof)      ! Temporary work array
  REAL :: temp1(lvlprof)      ! Temporary work array
  REAL :: temp2(lvlprof)      ! Temporary work array
  REAL :: temp3(lvlprof)      ! Temporary work array
  REAL :: temp4(lvlprof)      ! Temporary work array
  REAL :: temp5(lvlprof)      ! Temporary work array

  LOGICAL :: wk               ! Flag for W&K analytical sounding
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: lvlsnd
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Initialize the sounding profiles (analytical or from external
!  data).  All arguments are output.
!
!-----------------------------------------------------------------------
!
  CALL soundg(uprof,vprof,ptprof,pprof,qvprof,lvlprof,                  &
              zsnd,lvlsnd, temp1,temp2,temp3,temp4,temp5,wk )

!
!-----------------------------------------------------------------------
!
!  Interpolate the sounding profile at levels zsnd to zprof.
!  Outputs are also stored in arrays *prof, where * = u,v,w, etc.
!
!-----------------------------------------------------------------------
!
  CALL sndintrp(uprof,vprof,ptprof,pprof,qvprof,zsnd,lvlsnd,            &
                zprof,lvlprof, temp1)

  RETURN
END SUBROUTINE zprofil

!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE SOUNDG                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE soundg(usnd,vsnd,ptsnd,psnd,qvsnd,lvlprof,                   & 1,14
           zsnd,lvlsnd, tsnd,rhsnd,ptvsnd,pisnd,dpsnd,wk)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Specify the sounding profiles for u, v, pt, p and qv from analytical
!  functions or from an external sounding file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/17/1991.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/03/92 (M. Xue)
!  Further documentation.
!
!  2/10/93 (K. Droegemeier)
!  Cleaned up documentation.
!
!  9/10/94 (Weygandt & Y. Lu)
!  Cleaned up documentation.
!
!  10/31/94 (E. Adlerman)
!  Added option for dewpoint in input sounding
!  Added option for Weisman & Klemp sounding
!
!  9/22/1995 (M. xue)
!  Clarified the definition of sounding data recard 7.
!
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    lvlprof  The dimension of arrays usnd,vsnd,ptsnd,psnd,qvsnd,zsnd
!             This parameter is the number of levels in the
!             intermediate interpolated sounding.
!    WK       Logical flag for Weisman&Klemp thermodynamical
!             profile. Winds are set to zero.
!  OUTPUT:
!
!    usnd     u-velocity sounding data at original data levels (m/s)
!    vsnd     v-velocity sounding data at original data levels (m/s)
!    ptsnd    Potential temperature sounding data at original data
!             levels (K)
!    psnd     Pressure sounding data at original data levels (Pascal)
!    qvsnd    Water vapor specific humidity sounding data at the
!             original data levels (kg/kg)
!    zsnd     The height (above ground) of the sounding data levels
!             (m)
!    lvlsnd   The number of vertical levels in the input sounding
!
!  WORK ARRAYS:
!
!    tsnd     Temporary work array
!    ptvsnd   Temporary work array
!    pisnd    Temporary work array
!    rhsnd    Temporary work array
!    tdsnd    Temporary work array
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'mp.inc'          ! Message passing parameters.

  INTEGER :: lvlprof

  REAL :: usnd (lvlprof)     ! u-velocity from sounding data (m/s)
  REAL :: vsnd (lvlprof)     ! v-velocity from sounding data (m/s)
  REAL :: ptsnd(lvlprof)     ! Potential temperature sounding data (K)
  REAL :: psnd (lvlprof)     ! Pressure from sounding data (Pascal)
  REAL :: qvsnd(lvlprof)     ! Water vapor specific humidity
                             ! from sounding data (kg/kg)
  REAL :: zsnd (lvlprof)     ! The height of the sounding data levels
                             ! (m)

  REAL :: tsnd  (lvlprof)    ! Temporary work array
  REAL :: ptvsnd(lvlprof)    ! Temporary work array
  REAL :: pisnd (lvlprof)    ! Temporary work array
  REAL :: rhsnd (lvlprof)    ! Temporary work array
  REAL :: dpsnd (lvlprof)    ! Temporary work array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
!
!  The following three strings designate the type of sounding data.
!
!  if height(1:1)='h' or 'H', the sounding is given on height levels
!  if height(1:1)='p' or 'P', the sounding is given on pressure levels
!
!  if therm(1:1) ='t' or 'T', the sounding is specified in temperature
!  if therm(1:1) ='p' or 'P', the sounding is specified in potential
!                             temperature.
!
!  if humid(1:1) ='s' or 'S', the soundings uses specific humidity
!  if humid(1:1) ='r' or 'R', the sounding uses relative humidity,
!  if humid(1:1) ='d' or 'D', the soundings uses dewpoint temperature,
!
!  if wind(1:1) ='x' or 'Y', the sounding is specified in x and y
!                            component of velocity (u and v)
!  if wind(1:1) ='d' or 'D', the sounding is specified in direction
!                            and speed in m/s.
!  if wind(1:1) ='k' or 'K', the sounding is specified in direction
!                            and speed in knots.
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=30) :: height, therm, humid , wind
  LOGICAL :: in_is_p, in_is_t, in_is_rh, in_is_dp, wk
  CHARACTER (LEN=80) :: dummy
  CHARACTER (LEN=72) :: stime, sdate, sloc
  INTEGER :: k,iter,niter,lvlsnd,istat
  REAL :: psnd1,zsnd1
  INTEGER :: lenstr
  REAL :: p0inv,cpdrd
  REAL :: pt0, pttrop, ttrop, htrop, qvmix, mixtop
  INTEGER :: item1, item2
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_qvsat

!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!*$*  inline routine (f_qvsat)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  p0inv=1.0/p0
  cpdrd=1.0/rddcp
!
!-----------------------------------------------------------------------
!
!  Create Weisman & Klemp thermodynamic profile over 20km. depth.
!  Winds are set to zero. See W&K, JAS 1982 for details.
!
!-----------------------------------------------------------------------
!

  IF ( wk ) THEN

    lvlsnd = 201
    zsnd1  = 0.0
    psnd1  = 100000.0
    height = 'height'
    humid  = 'rh'
    therm  = 'pt'
    wind   = 'uv'

    pttrop = 343.0           ! Tropopause pot.temp.
    ttrop  = 213.0           ! Tropopause temp.
    pt0    = 300.0           ! Surface pot.temp.
    htrop  = 12000.0         ! Tropopause height
    qvmix  = 0.015           ! Mixed layer mixing ratio
    mixtop = 1200.0          ! Mixed layer height

    DO k=1, lvlsnd
      usnd(k) = 0.0
      vsnd(k) = 0.0
    END DO

    zsnd(1) = zsnd1
    DO k=2, lvlsnd
      zsnd(k)=(20000.0/REAL(lvlsnd-1)) * REAL(k-1)
    END DO

    DO k=1, lvlsnd
      IF ( zsnd(k) <= htrop ) THEN
        ptsnd(k) =                                                      &
            pt0 + (pttrop-pt0)*( (zsnd(k)/htrop)**(5.0/4.0) )
        qvsnd(k) =                                                      &
            1.0 - 0.75*( (zsnd(k)/htrop)**(5.0/4.0) )
      ELSE
        ptsnd(k) = pttrop*EXP( (zsnd(k)-htrop)*g/(ttrop*cp) )
        qvsnd(k) = 0.25
      END IF
    END DO


  ELSE

!
!-----------------------------------------------------------------------
!
!              Sounding File Format for ARPS 3.0
!
!  Record 1: a header line, such as "1-D Sounding Input for ARPS"
!            (skipped)
!  Record 2: miscellaneous description of sounding (skipped)
!  Record 3: time of sounding (character*72)
!  Record 4: date of sounding (character*72)
!  Record 5: location of sounding (character*72)
!  Record 6: three character strings designate the sounding
!            data type, e.g.
!            'pressure' 'potential_temperature' 'relative_humidity'.
!            Only the first character of the strings
!            is decoded, and thus the first character should not be
!            left blank.  Note that either upper or lower case may be
!            used.  A more detailed explanation is provided in the
!            portion of the code where these strings are declared.
!  Record 7: Ground-level height (m) and the correspoding pressure
!            (Pascal) when sounding is specified at height levels.
!            When it is given on pressure levels, this record is
!            not used. The last sounding data is assumed to be at
!            the ground level (z=0).
!  Record 8: number of data levels in sounding (variable lvlsnd )
!  Record 9: a line of data column labels (not read in)
!  Record 10 to the end:
!            sounding data in the order 'z/p, pt/t, qv/rh, u, v'
!
!  For records 10  to the end, there is one line of data
!            corresponding to each sounding level.
!
!  Important convention:
!            Line 10 corresponds to the level k = lvlsnd
!            Line 11 corresponds to the level k = (lvlsnd -1)
!            etc.
!
!  Units of the data:
!            pressure: Pascals, height: meters, temperature, dewpoint,
!            and potential temperature: degrees Kelvin, specific
!            humidity kg/kg, relative humidity: 0 to 1.
!
!  CAUTION: The sounding data have to be listed in the order of
!           decreasing height or increasing pressure.
!
!-----------------------------------------------------------------------
!

    istat = 0

!
!-----------------------------------------------------------------------
!
!  Open the sounding file
!
!-----------------------------------------------------------------------
!
    CALL getunit( sfunit )

    OPEN(UNIT=sfunit, FILE=trim(sndfile), STATUS='old',                 &
                IOSTAT=istat)

    IF(istat /= 0) THEN
      IF (myproc == 0 ) THEN
        WRITE (6,1000) 'Error in opening sounding file'
        WRITE (6,'(1x,a,a)') 'Check file ',sndfile
      END IF
      CALL arpsstop('arpsstop called from INIBASE with opening snd file',1)
    ELSE
      lenstr = 80
      CALL strlnth( sndfile, lenstr )
      IF (myproc == 0 ) &
         WRITE(6,'(1x,/2a,a,i3/)') 'Sounding file ', sndfile(1:lenstr),    &
              ' opened using FORTRAN unit',sfunit
    END IF

    READ (sfunit,1000,IOSTAT=istat) dummy ! header line
    IF (myproc ==0) WRITE(6,1000)dummy

    READ (sfunit,1000,IOSTAT=istat) dummy ! miscellaneous descriptions
    IF (myproc ==0) WRITE(6,1000)dummy

    READ (sfunit,1000,IOSTAT=istat) stime ! time of sounding
    IF (myproc ==0) WRITE(6,1000)stime


    READ (sfunit,1000,IOSTAT=istat) sdate ! date of sounding
    IF (myproc ==0) WRITE(6,1000)sdate

    READ (sfunit,1000,IOSTAT=istat) sloc  ! location of sounding
    IF (myproc ==0) WRITE(6,1000)sloc
!
!-----------------------------------------------------------------------
!
!  The following is for maintaining backwards compatibility with
!  ARPS 3.* sounding files which did not have a wind data type
!  indicator.
!
!-----------------------------------------------------------------------
!
    READ (sfunit,1111,IOSTAT=istat) dummy ! tem string for data types
    1111  FORMAT(a)

    item1 = 0
    item1 = INDEX( dummy(item1+1:80), '''' ) + item1
    item2 = INDEX( dummy(item1+1:80), '''' ) + item1
    height = dummy((item1+1):item2-1)

    item1 = item2
    item1 = INDEX( dummy(item1+1:80), '''' ) + item1
    item2 = INDEX( dummy(item1+1:80), '''' ) + item1
    therm = dummy(item1+1:item2-1)

    item1 = item2
    item1 = INDEX( dummy(item1+1:80), '''' ) + item1
    item2 = INDEX( dummy(item1+1:80), '''' ) + item1
    humid = dummy((item1+1):item2-1)

    item1 = item2
    item1 = INDEX( dummy(item1+1:80), '''' ) + item1
    item2 = INDEX( dummy(item1+1:80), '''' ) + item1
    IF ( item1 == item2 ) THEN
      wind = 'uv'
    ELSE
      wind  = dummy(item1+1:item2-1)
    END IF

    IF (myproc == 0 ) WRITE(6,'(1x,4a)') height, therm, humid, wind

    READ (sfunit,*,IOSTAT=istat) zsnd1,psnd1
    IF (myproc ==0)  &
       WRITE(6,'(1x,''SNDL='',F13.2,''PSNDL='',F13.2)') zsnd1,psnd1

    READ (sfunit,*,IOSTAT=istat) lvlsnd   ! actual number of data
                                          ! levels in the sounding
                                          ! file
    WRITE(6,'(1x,i6)')

    IF(lvlsnd > lvlprof) THEN

      IF (myproc == 0 ) THEN
        WRITE (6,1000) 'Error in sounding file'
        WRITE (6,'(1x,a,i6)') 'lvlsnd cannot be greater than ',lvlprof
        WRITE (6,1000) 'Check subroutine SOUNDG'
      END IF

      CALL arpsstop('arpsstop called from INIBASE with number of levels &
        & in the sounding ',1)

    END IF


    READ (sfunit, 1000, IOSTAT=istat) dummy ! a separator line,
    IF (myproc ==0) WRITE(6,1000)dummy

    IF(istat /= 0) THEN

      IF (myproc == 0 ) THEN
        WRITE (6,1000) 'Error in sounding file'
        WRITE (6,1000) 'Error in sounding file header'
        WRITE (6,1000) 'Check sounding file', sndfile
      END IF

      CALL arpsstop('arpsstop called from INIBASE with opening snd file',1)

    END IF

!-----------------------------------------------------------------------
!
!  Read in the sounding data:
!
!-----------------------------------------------------------------------

    DO k = lvlsnd ,1,-1

      READ (sfunit, *, IOSTAT=istat)                                    &
            zsnd(k), ptsnd(k), qvsnd(k), usnd(k), vsnd(k)

      IF (myproc == 0) WRITE(6,'(1x,i4,5f13.6)')  &
          k,zsnd(k), ptsnd(k), qvsnd(k), usnd(k), vsnd(k)

      IF(istat /= 0) THEN

        IF (myproc == 0) THEN
          WRITE (6,1000) 'Error in sounding file'
          WRITE (6,1000) 'Error in sounding file data'
          WRITE (6,1000) 'Check ',sndfile
        END IF
        CALL arpsstop('arpsstop called from INIBASE opening snd file',1)

      END IF

    END DO

    CLOSE (UNIT=sfunit)      ! close the sounding file
    CALL retunit( sfunit )

    IF (myproc == 0) WRITE(6,'(/1x,a)')'Sounding data reading complete.'
!
!-----------------------------------------------------------------------
!
!  End of sounding data input.
!
!-----------------------------------------------------------------------
!

  END IF
!
!-----------------------------------------------------------------------
!
!  If wind input is direction and speed, then find u,v.
!  For now there is no rotation of winds due to the sounding
!  not coinciding with the map projection longitude.
!  We tell the conversion program the sounding is at trulat1
!  and trulon.  In the future this can be the actual sounding
!  lat, lon or the center of the ARPS grid.
!  Here tsnd, pisnd and ptvsnd are used as dummy work arrays.
!
!-----------------------------------------------------------------------
!
  IF( wind(1:1) == 'd' .OR. wind(1:1) == 'D' .OR.                       &
        wind(1:1) == 'k' .OR. wind(1:1) == 'K' ) THEN

    DO k=1,lvlsnd
      pisnd(k)=trulon
    END DO

    IF( wind(1:1) == 'k' .OR. wind(1:1) == 'K' ) THEN
      DO k=1,lvlsnd
        vsnd(k) = vsnd(k)*0.514444
      END DO
    END IF

    CALL ddrotuv(lvlsnd,pisnd,usnd,vsnd,ptvsnd,                         &
                 usnd,vsnd)

  END IF
!
!-----------------------------------------------------------------------
!
!  If input data are not in height units, copy into the pressure
!  array and set the input data type flag.
!
!-----------------------------------------------------------------------
!

  IF( height(1:1) /= 'h' .AND. height(1:1) /= 'H' ) THEN

    in_is_p = .true.
    DO k=1,lvlsnd
      psnd(k) =zsnd(k)
    END DO

  ELSE

    in_is_p = .false.

  END IF
!
!-----------------------------------------------------------------------
!
!  If input is temperature, copy it from the potential temperature
!  array and set the input data type flag.
!
!-----------------------------------------------------------------------
!
  IF( therm(1:1) == 't' .OR. therm(1:1) == 'T' ) THEN

    in_is_t = .true.
    DO k=1,lvlsnd
      tsnd(k) =ptsnd(k)
    END DO

  ELSE

    in_is_t = .false.

  END IF

!
!-----------------------------------------------------------------------
!
!  If input is relative humidity, copy it from the specific humidity
!  array and set the input data type flag. If input is dewpoint, copy
!  it from the specific humidity array and set the input data type
!  flag. Otherwise, humidity is already mixing ratio/specific humid.
!  and the data flag types are set appropriately.
!
!-----------------------------------------------------------------------
!
  IF( humid(1:1) == 'r' .OR. humid(1:1) == 'R' ) THEN


    in_is_rh = .true.
    in_is_dp = .false.
    DO k=1,lvlsnd
      rhsnd(k) =qvsnd(k)
    END DO

  ELSE IF ( humid(1:1) == 'd' .OR. humid(1:1) == 'D' ) THEN

    in_is_dp = .true.
    in_is_rh = .false.
    DO k=1, lvlsnd
      dpsnd(k) = qvsnd(k)
    END DO

  ELSE

    in_is_dp = .false.
    in_is_rh = .false.
    DO k=1,lvlsnd
      rhsnd(k) = 0.0
      dpsnd(k) = 0.0
    END DO


  END IF

!
!-----------------------------------------------------------------------
!
!  Calculate the derived variables from the input sounding
!
!-----------------------------------------------------------------------
!

  IF( in_is_p ) THEN
!
!-----------------------------------------------------------------------
!
!  If the sounding data are specified on pressure levels...
!
!-----------------------------------------------------------------------
!
    DO k=1,lvlsnd

      IF( psnd(k) <= 0.0) THEN

        IF (myproc == 0) WRITE(6,'(2(/1x,a))')  &
             'Negative pressure found in the sounding profile.',        &
             'Job stopped. Check subroutine INIBASE.'
        CALL arpsstop('arpsstop called from INIBASE due to neg. pressure ',1)

      END IF

    END DO

    IF( in_is_t ) THEN

      DO k=1,lvlsnd
        ptsnd(k) = tsnd(k)*(p0/psnd(k))**rddcp
      END DO

    ELSE

      DO k=1,lvlsnd
        tsnd(k) = ptsnd(k)*(psnd(k)*p0inv)**rddcp
      END DO

    END IF

    IF( in_is_rh ) THEN

      CALL getqvs(1,1,lvlsnd, 1,1,1,1,1,lvlsnd, psnd,tsnd,qvsnd)

      DO k=1,lvlsnd
        qvsnd(k) = qvsnd(k) * rhsnd(k)
      END DO

    ELSE IF ( in_is_dp ) THEN

      CALL getqvs(1,1,lvlsnd, 1,1,1,1,1,lvlsnd, psnd,dpsnd,qvsnd)

    END IF
!
!-----------------------------------------------------------------------
!
!  The virtual temperature Tv:
!
!-----------------------------------------------------------------------
!
    DO k = 1,lvlsnd
      ptvsnd(k) = tsnd(k)*(1.0+rvdrd*qvsnd(k))/(1.0+qvsnd(k))
    END DO

!
!-----------------------------------------------------------------------
!
!  Integrate hydrostatic relation dz/dp = - RTv/(pg) for z
!
!-----------------------------------------------------------------------
!
    zsnd(1) = zsnd1

    DO k = 2,lvlsnd
      zsnd(k)=zsnd(k-1) - rd/g                                          &
            *(ptvsnd(k)+ptvsnd(k-1))/(psnd(k)+psnd(k-1))                &
            *(psnd(k)-psnd(k-1))
    END DO
!

  ELSE IF( .NOT. in_is_p ) THEN

!
!-----------------------------------------------------------------------
!
!  If the sounding data are specified at height levels
!
!-----------------------------------------------------------------------
!
    IF( in_is_rh .OR. in_is_dp ) THEN
!
!-----------------------------------------------------------------------
!
!  If the input fields are height and relative humidity/dewpoint, the
!  pressure, temperature and specific humidity profiles are coupled;
!  thus, one must iterate to compute the pressure and temperature.
!  We use a predetermined number (10) of iterations, which should
!  ensure convergence in even the most difficult situations.
!
!-----------------------------------------------------------------------
!
      niter = 10

      DO k=1,lvlsnd
        qvsnd(k) =0.0
      END DO
    ELSE
      niter = 1
    END IF


    DO iter = 1, niter

      IF( .NOT. in_is_t) THEN
!
!-----------------------------------------------------------------------
!
!  Thermodynamic variable is potential temperature...
!
!-----------------------------------------------------------------------
!

!-----------------------------------------------------------------------
!
!  Calculate ptvsnd, the virtual potential temperature sounding array
!
!-----------------------------------------------------------------------
!
        DO k = 1,lvlsnd
          ptvsnd(k) = ptsnd(k)*(1.0+rvdrd*qvsnd(k))/(1.0+qvsnd(k))
        END DO

!
!-----------------------------------------------------------------------
!
!  Integrate the hydrostatic relation to get PISND, the
!  nondimensional pressure: d(pi)/dz = -g/(cp * ptv)
!
!-----------------------------------------------------------------------
!
!

        pisnd(1) = (psnd1*p0inv)**rddcp

        DO k = 2,lvlsnd
          pisnd(k) = pisnd(k-1) - g/cp                                  &
                     *(zsnd(k)-zsnd(k-1))                               &
                     /(0.5*(ptvsnd(k)+ptvsnd(k-1)))
        END DO

        DO k = 1,lvlsnd
          psnd(k) = p0 * (pisnd(k)**cpdrd)
          tsnd(k) = ptsnd(k) * pisnd(k)
        END DO

      ELSE
!
!-----------------------------------------------------------------------
!
!  Thermodynamic variable is air temperature...
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  Calculate the virtual temperature and store it in ptvsnd
!
!-----------------------------------------------------------------------
!
        DO k = 1,lvlsnd
          ptvsnd(k) = tsnd(k)*(1.0+rvdrd*qvsnd(k))/(1.0+qvsnd(k))
        END DO

!
!-----------------------------------------------------------------------
!
!  Calculate log(p) according to  d(log(p))/dz = -g/(r * tv)
!
!-----------------------------------------------------------------------
!
        pisnd(1) = LOG( psnd1 )

        DO k = 2,lvlsnd
          pisnd(k) = pisnd(k-1) - g/rd                                  &
                     *(zsnd(k)-zsnd(k-1))                               &
                     /(0.5*(ptvsnd(k)+ptvsnd(k-1)))
        END DO

        DO k = 1,lvlsnd
          psnd(k) = EXP( pisnd(k) )
          ptsnd(k) = tsnd(k) * (p0/psnd(k))**rddcp
        END DO

      END IF

      IF ( in_is_rh ) THEN

        CALL getqvs(1,1,lvlsnd, 1,1,1,1,1,lvlsnd, psnd,tsnd,qvsnd)

        DO k=1,lvlsnd
          qvsnd(k) = qvsnd(k) * rhsnd(k)
        END DO

      ELSE IF ( in_is_dp ) THEN

        CALL getqvs(1,1,lvlsnd, 1,1,1,1,1,lvlsnd, psnd,dpsnd,qvsnd)

      END IF

    END DO

  END IF

  IF( .NOT. in_is_rh ) THEN

    CALL getqvs(1,1,lvlsnd, 1,1,1,1,1,lvlsnd, psnd,tsnd,rhsnd)

    DO k=1,lvlsnd
      rhsnd(k) = qvsnd(k)/rhsnd(k)
    END DO

  END IF

!
!-----------------------------------------------------------------------
!
!  Add fixed mixed layer to WK profile and recalculate RH values
!
!-----------------------------------------------------------------------
!
  IF ( wk ) THEN

    DO k=1, lvlsnd

      IF (zsnd(k) <= mixtop) THEN
        qvsnd(k) = qvmix
        rhsnd(k) = qvmix / f_qvsat( psnd(k), tsnd(k) )
      END IF

    END DO

  END IF


!
!-----------------------------------------------------------------------
!
!  Write out the uninterpolated sounding profile
!
!-----------------------------------------------------------------------
!
  IF (myproc == 0) THEN
    WRITE(6,'(/1x,a/)')                                                   &
      'Base state profiles at original sounding data levels:'

    WRITE(6,'(1x,2a)')                                                    &
       '   n    z(m)   p(Pascal)  PT(K)   T(K)',                        &
       '    Qv(kg/kg)   RH   ubar(m/s) vbar(m/s)'

    DO k=lvlsnd,1,-1
      WRITE(6,'(1x,i4,f9.2,f9.1,2f9.3,f9.6,3f9.3)')                       &
          k,zsnd(k),psnd(k),ptsnd(k),                                   &
          tsnd(k), qvsnd(k), rhsnd(k), usnd(k), vsnd(k)
    END DO
  END IF

  RETURN

  1000  FORMAT(1X,a)

END SUBROUTINE soundg

!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE SNDINTRP                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE sndintrp(uprof,vprof,ptprof,pprof,qvprof,zsnd,lvlsnd,        & 1,5
           zprof,lvlprof, temp1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Interpolate the input sounding profile on the grid "zsnd" to a
!  profile on the model grid "zprof".
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/17/1991.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/03/92 (M. Xue)
!  Further documentation.
!
!  2/10/93 (K. Droegemeier)
!  Cleaned up documentation.
!
!  9/10/94 (Weygandt & Y. Lu)
!  Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    uprof    u-velocity sounding data at original data levels (m/s)
!    vprof    v-velocity sounding data at original data levels (m/s)
!    ptprof   Potential temperature sounding data at original data
!             levels (K)
!    pprof    Pressure sounding data at original data levels (Pascal)
!    qvprof   Water vapor specific humidity sounding data at the
!             original data levels (kg/kg)
!    zsnd     The height of the input sounding data levels (m)
!    lvlsnd   The number of input sounding data levels
!
!    zprof    The height of equally spaced grid levels to which
!             the sounding data are interpolated (m).
!    lvlprof  The number of such grid levels.
!
!  OUTPUT:
!
!    uprof    u-velocity data defined on grid levels given in zprof
!    vprof    v-velocity data defined on grid levels given in zprof
!    ptprof   Potential temperature data defined on grid levels
!             given in zprof
!    pprof    Pressure data defined on grid levels given in zprof
!    qvprof   Water vapor specific humidity data defined on grid
!             levels given in zprof
!
!  WORK ARRAY:
!
!    temp1     Temporary work array
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: lvlprof, lvlsnd

  REAL :: uprof (lvlprof)     ! Profile of u-velocity defined at
                              ! levels given in zprof (m)
  REAL :: vprof (lvlprof)     ! Profile of v-velocity defined at
                              ! levels given in zprof (m)
  REAL :: ptprof(lvlprof)     ! Profile of potential temperature
                              ! defined at levels given in zprof (K)
  REAL :: pprof(lvlprof)      ! Profile of pressure defined at
                              ! levels given in zprof (Pascal)
  REAL :: qvprof(lvlprof)     ! Profile of water vapor specific
                              ! humidity defined at levels given in
                              ! zprof (kg/kg)
  REAL :: zsnd (lvlprof)      ! Temporary work array
  REAL :: zprof (lvlprof)     ! 1-D array defining the regularly-spaced
                              ! vertical grid levels at which base
                              ! state
                              ! profiles will be sampled (m)
  REAL :: temp1 (lvlsnd)      ! Temporary work array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: k
!
!-----------------------------------------------------------------------
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=1,lvlsnd
    temp1(k) = uprof(k)
  END DO

  CALL inte1d(temp1,zsnd,lvlsnd, uprof,zprof,lvlprof)

  DO k=1,lvlsnd
    temp1(k) = vprof(k)
  END DO

  CALL inte1d(temp1,zsnd,lvlsnd, vprof,zprof,lvlprof)

  DO k=1,lvlsnd
    temp1(k) = ptprof(k)
  END DO

  CALL inte1d(temp1,zsnd,lvlsnd, ptprof,zprof,lvlprof)

  DO k=1,lvlsnd
    temp1(k) = pprof(k)
  END DO

  CALL inte1d(temp1,zsnd,lvlsnd, pprof,zprof,lvlprof)

  DO k=1,lvlsnd
    temp1(k) = qvprof(k)
  END DO

  CALL inte1d(temp1,zsnd,lvlsnd, qvprof,zprof,lvlprof)

  RETURN
END SUBROUTINE sndintrp
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE INTE1D                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE inte1d(fdat,zdat,ndat,fpro,zpro,npro) 10
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Interpolate data from sounding to a model grid-column with uniform
!  grid spacing.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/17/1991.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/03/92 (M. Xue)
!  Further documentation.
!
!  2/10/93 (K. Droegemeier)
!  Cleaned up documentation.
!
!  9/10/94 (Weygandt & Y. Lu)
!  Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
!  Input :
!
!    fdat     1-D array defined at levels in zdat to be
!             interpolated to levels defined by zpro.
!    zdat     The height of the input data given by fdat.
!    ndat     The number of data in array fdat.
!
!    zpro     The grid level height to which data are interpolated
!    npro     The number of interpolated data levels
!
!  Output:
!
!    fpro     Output data interpolated to grid levels defined by zpro.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: ndat             ! The number of data in array fdat
  INTEGER :: npro             ! The number of interpolated data levels

  REAL :: fdat(ndat)          ! 1-D array defined at levels in zdat to
                              ! be interpolated to levels defined by
                              ! zpro
  REAL :: zdat(ndat)          ! The height of the input data given by
                              ! fdat
  REAL :: zpro(npro)          ! The grid level height to which data are
                              ! interpolated
  REAL :: fpro(npro)          ! Output data interpolated to grid levels
                              ! defined by zpro
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i0,i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  i0=1
  DO j=1,npro

    IF(zpro(j) <= zdat(1))THEN
      fpro(j)=fdat(1)
      CYCLE
    END IF

    IF(zpro(j) >= zdat(ndat))THEN
      fpro(j)=fdat(ndat)
      CYCLE
    END IF

!
!-----------------------------------------------------------------------
!
!  Find index for interpolation
!
!-----------------------------------------------------------------------
!
    DO i=i0,ndat-1
      IF(zpro(j) >= zdat(i).AND.zpro(j) < zdat(i+1))GO TO 15
    END DO

    CYCLE
    15      i0=i

    fpro(j)=fdat(i)+(fdat(i+1)-fdat(i))*(zpro(j)-zdat(i))               &
            /(zdat(i+1)-zdat(i))

  END DO

  RETURN
END SUBROUTINE inte1d