c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE INIBASE                    ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE INIBASE(nx,ny,nz, ubar,vbar,ptbar,pbar,rhobar,qvbar)

c
c#######################################################################
c
c     PURPOSE:
c
C     INIitialize the BASE state variables.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     3/17/1991.
c
c     MODIFICATION HISTORY:
c
c     5/02/92 (M. Xue)
c     Added full documentation.
C
c     5/03/92 (M. Xue)
c     Further documentation. 
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c     OUTPUT:
c
c       ubar     Base state zonal velocity component (m/s)
c       vbar     Base state meridional velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state density (kg/m**3)
c       qvbar    Base state water vapor specific humidity (kg/kg)
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      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 rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
      real qvbar (nx,ny,nz)     ! Base state water specific humidity
                                ! (kg/kg).
c
      integer lvlprof
      real    depthp

      parameter( lvlprof=601, depthp = 3.0e4 )

c
c#######################################################################
c
c     lvlprof: 
c
c     The number of levels of the equally spaced profiles. 
c     lvlprof is chosen by the user but must be no smaller than the 
c     number of levels in a input sounding
c 
c     The greater the number of levels , the better these profiles 
c     represent the original sounding. Typically 200 levels is good 
c     enough for a model run with about 50 grid levels.
c
c     depthp:  
c
c     The depth of atmosphere these profiles will cover. 
c     depth >= (nz-3)*dz, i.e. depth should be larger than model depth. 
c     Otherwise, extrapolations will occur for grid points higher 
c     than this depth. 
c
c#######################################################################
c 

      real uprof (lvlprof)      ! Temporary array
      real vprof (lvlprof)      ! Temporary array
      real ptprof(lvlprof)      ! Temporary array
      real pprof (lvlprof)      ! Temporary array
      real qvprof(lvlprof)      ! Temporary array
      real zprof (lvlprof)      ! Temporary array
c
      real temp1(lvlprof)       ! Temporary array
      real temp2(lvlprof)       ! Temporary array
      real temp3(lvlprof)       ! Temporary array
      real temp4(lvlprof)       ! Temporary array
      real temp5(lvlprof)       ! Temporary array
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      real factor,zsgrid,pibar,zgrid,zugrid,zvgrid,deltaz,pt0,pi,z
      integer i,j,k,kdata,isentrop
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
      include 'phycst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external zprofil
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      isentrop = 0
c
c#######################################################################
c
c     If the atmosphere is isentropic, i.e. ptbar = constant
c     then specify the base state directly from analytical functions
c
c#######################################################################
c
      IF( isentrop.eq.1) THEN 
c
c#######################################################################
c
c     To obtain sounding profiles specified on equally spaced
c     grid levels.
c
c#######################################################################
c
        pt0 = 300.0

        DO 100 k=1,nz-1
        DO 100 j=1,ny-1
        DO 100 i=1,nx-1
          z = (k-1.5)*dz
          pi = 1.0-g/(cp*pt0)*z
          pbar(i,j,k) = p0 * pi**(cp/rd)
          ptbar(i,j,k)=pt0
          rhobar(i,j,k)=pbar(i,j,k)/(rd*pt0*pi)
          qvbar(i,j,k)=0.0
100     CONTINUE

        DO 200 k=1,nz-1
        DO 200 j=1,ny-1
        DO 200 i=1,nx
          ubar(i,j,k)=0.0
200     CONTINUE

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

      ELSE
c
c#######################################################################
c
c     Specify the base state from discrete profiles
c
c#######################################################################
c

        deltaz = depthp/(lvlprof-1)

        write(6,'(/'' The input sounding data is sampled at a '',
     :            f7.2,'' meter interval.''/)') deltaz

        DO 400 k=1,lvlprof
          zprof(k) = (k-1)*deltaz
400     CONTINUE 
c
c#######################################################################
c
c     To obtain the vertical profiles of u, v, potential temperature,
c     pressure and water vapor specific humidity of the base state
c     atmospere sampled as vertical grid levels defined by zprof.
c
c#######################################################################
c
        CALL zprofil(zprof,lvlprof,
     :               uprof,vprof,ptprof,pprof,qvprof, 
     :               temp1,temp2,temp3,temp4,temp5)
c
c
c#######################################################################
c
c     Linearly interpolating the profiles to model grid levels
c     for base sate varaibles ubar, vbar, ptbar, and pbar.
c
c#######################################################################
c

c#######################################################################
c
c     Interpolation for scalars. 
c
c#######################################################################

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

          zsgrid = (k-1.5)*dz
          kdata = int( zsgrid/deltaz ) + 1
          kdata = min( lvlprof-1, max(1,kdata) )
          factor = (zsgrid-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 = (pbar(i,j,k)/p0)**(rddcp)
          rhobar(i,j,k)=pbar(i,j,k)/(rd*ptbar(i,j,k)*pibar)
10      CONTINUE

c#######################################################################
c
c     Interpolation for ubar, and subtract the x component of the
c     Domain translation speed from the base state x-velocity.
c
c#######################################################################

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

          zugrid = (k-1.5)*dz
          kdata = int( zugrid/deltaz ) + 1
          kdata = min( lvlprof-1, max(1,kdata) )
          factor = (zugrid-zprof(kdata))/deltaz

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

          ubar(i,j,k) = ubar(i,j,k) - umove
20    CONTINUE

c#######################################################################
c
c     Interpolation for vbar, and subtract the y component of the 
c     Domain translation speed from the base state y-velocity. 
c
c#######################################################################

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

          zvgrid = (k-1.5)*dz
          kdata = int( zvgrid/deltaz ) + 1
          kdata = min( lvlprof-1, max(1,kdata) )
          factor = (zvgrid-zprof(kdata))/deltaz

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

          vbar(i,j,k) = vbar(i,j,k) - vmove
30      CONTINUE
   
      ENDIF
c
c#######################################################################
c
c     Set top and bottom boundary conditions for the base state variables.  
c
c#######################################################################
c
      DO 50 j=1,ny-1
      DO 50 i=1,nx-1
        qvbar(i,j,1)=qvbar(i,j,2)
        qvbar(i,j,nz-1)=qvbar(i,j,nz-2)
50    CONTINUE

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

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

c
c#######################################################################
c
c     Print out a list of sounding. 
c
c#######################################################################
c
      i = (nx-1)/2+1
      j = (ny-1)/2+1

      write(6,'(/a,i4,a,i4/)')' Base state z-profiles at I=',i,', J=',j

      write(6,'(1x,2a)') 
     :     '   k      z(m)    pbar(Pascal)  ptbar(K)',
     :     '   Qvbar(kg/kg)   U(m/s)     V(m/s)'

      DO 1100 k=nz-1,1,-1

        zgrid = (k-1.5)*dz
        write(6,'(1x,i4,2f12.3,f12.5,f12.8,2f12.5)') 
     :        k,zgrid, pbar(i,j,k),
     :        ptbar(i,j,k),qvbar(i,j,k),ubar(i,j,k), vbar(i,j,k)

1100  CONTINUE

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE ZPROFIL                    ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE ZPROFIL(zprof,lvlprof, 
     :           uprof,vprof,ptprof,pprof,qvprof,
     :           zsnd,temp1,temp2,temp3,temp4) 

c
c#######################################################################
c
c     PURPOSE:
c
c     Obtain the initial sounding profiles at regularly spaced grid 
c     levels. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     3/17/1991.
c
c     MODIFICATION HISTORY:
c
c     5/02/92 (M. Xue)
c     Added full documentation.
C
c     5/03/92 (M. Xue)
c     Further documentation. 
C
c#######################################################################
c
c     INPUT :
c
c       zprof     1-d array defining regularly spaced vertical grid 
c                 levels at which base state profiles will be sampled (m)
c       lvlprof   The number of such grid levels.
c
c     OUTPUT:
c
c       uprof     Profile of u-velocity defined at levels given 
c                 in zprof (m)
c       vprof     Profile of v-velocity defined at levels given 
c                 in zprof (m)
c       ptprof    Profile of potential temperature defined at levels 
c                 given in zprof (K)
c       pprof     Profile of pressure defined at levels given in zprof 
c                 (Pascal)
c       qvprof    Profile of water vapor specific humidity defined at 
c                 levels given in zprof (kg/kg)
c
c     WORK ARRAYS:
c
c       zsnd    A temporary array to be used to store the height 
c               actual input sounding data.
c       temp1   Locally defined array
c       temp2   Locally defined array
c       temp3   Locally defined array
c       temp4   Locally defined array
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer lvlprof

      real zprof (lvlprof)

      real uprof (lvlprof)
      real vprof (lvlprof)
      real ptprof(lvlprof)
      real pprof (lvlprof)
      real qvprof(lvlprof) 
c
      real zsnd (lvlprof)      ! 1-D temporary array
      real temp1(lvlprof)      ! 1-D temporary array
      real temp2(lvlprof)      ! 1-D temporary array
      real temp3(lvlprof)      ! 1-D temporary array
      real temp4(lvlprof)      ! 1-D temporary array
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer lvlsnd
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external soundg
      external sndintrp
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Initialize sounding profiles (analytical or from external data)
c     all arguments are output.
c
c#######################################################################
c
      CALL soundg(uprof,vprof,ptprof,pprof,qvprof,lvlprof,
     :            zsnd,lvlsnd, temp1,temp2,temp3,temp4 )

c
c#######################################################################
c
c     Interpolate sounding profiles at levels ZSND to levels zproF.
c     Output are also stored in arrays *PROF. 
c
c#######################################################################
c
      CALL sndintrp(uprof,vprof,ptprof,pprof,qvprof,zsnd,lvlsnd,
     :              zprof,lvlprof, temp1)

      RETURN
      END

c 
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE SOUNDG                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c

      SUBROUTINE SOUNDG(usnd,vsnd,ptsnd,psnd,qvsnd,lvlprof,
     :                  zsnd,lvlsnd, tsnd,rhsnd,ptvsnd,pisnd)

c
c#######################################################################
c
c     PURPOSE:
c
c     Specify the sounding profiles for u, v, pt, p and qv from analytical
c     functions or by reading from a sounding file.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     3/17/1991.
c
c     MODIFICATION HISTORY:
c
c     5/02/92 (M. Xue)
c     Added full documentation.
C
c     5/03/92 (M. Xue)
c     Further documentation. 
C
c#######################################################################
c
c     INPUT :
c
c       lvlprof  The dimension of arrays usnd,vsnd,ptsnd,psnd,qvsnd,zsnd
c
c     OUTPUT:
c
c       usnd     u-velocity sounding data at original data levels (m/s)
c       vsnd     v-velocity sounding data at original data levels (m/s)
c       ptsnd    Potential temperature sounding data at original data 
c                levels (K)
c       psnd     Pressure sounding data at original data levels (Pascal)
c       qvsnd    Water vapor specific humidity sounding data at the 
c                original data levels (kg/kg)
c       zsnd     The height of the sounding data levels (m)
c       lvlsnd   The number of sounding data
c
c     WORK ARRAYS:
c  
c       tsnd     Temporary array
c       ptvsnd   Temporary array
c       pisnd    Temporary array
c       rhsnd    Temporary array
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      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)
c 
      real tsnd  (lvlprof)    ! Temporary array
      real ptvsnd(lvlprof)    ! Temporary array
      real pisnd (lvlprof)    ! Temporary array
      real rhsnd (lvlprof)    ! Temporary array
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
c
c     The following three strings designate the type of sounding data.
c
c     if height(1:1)='h' or 'H', the sounding is given in geometric
c        height, otherwise in pressure.
c
c     if therm(1:1) ='t' or 'T', the sounding uses temperature instead
c        of potential temperature.
c
c     if humid(1:1) ='r' or 'R', the sounding uses relative humidity,
c        otherwise uses sepecific humidity or mixing ratio.
c
c#######################################################################
c
      character*30 height, therm, humid 
      logical in_is_p, in_is_t, in_is_rh
      character*80 dummy
      CHARACTER*72 STIME, SDATE, SLOC
      integer k,iter,niter,lvlsnd,istat
      real psnd1,zsnd1
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
      include 'phycst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external getqvs
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

c
c#######################################################################
c
c                 Sounding File Format for ARPS 3.0
c
c     Record 1: a header line, such as "1-D Sounding Input for ARPS" 
c               (skipped)
c     Record 2: miscellaneous description of sounding (skipped)
c     Record 3: time of sounding (character*72)
c     Record 4: date of sounding (character*72)
c     Record 5: location of sounding (character*72)
c     Record 6: three character strings designating the sounding 
c               data type, e.g. 
c               'pressure' 'potential_temperature' 'relative_humity'. 
c               Here, only the first character of the strings 
c               is significant and in case insensitive. 
c               For this reason, the first character should not be 
c               left blank. More detailed explanations are given 
c               in the declarations  of these strings. 
c     Record 7: height (m) and pressure (Pascal) at the first (lowest) 
c               sounding data level 
c     Record 8: number of data levels in sounding (variable lvlsnd ) 
c     Record 9: a line of data column labels. skipped when reading. 
c     Record 10 to the end:
c               sounding data in the order 'z/p, pt/t, qv/rh, u, v'     
c
c     For records 10  to the end, there is one line of data corresponding to 
c               each level of the sounding.
c
c     Important convention:
c               Line 10 corresponds to the level K = lvlsnd 
c               Line 11 corresponds to the level K = (lvlsnd -1)
c               etc.
c
c     Units of the data: 
c               pressure: Pascal, height: meter, temperature and 
c               potential temperature: Kelvin, specific humidity kg/kg, 
c               relative humidity: 0 to 1 scale. 
c
c#######################################################################
c

      istat = 0
c
c#######################################################################
c
c     Open the sounding file
c
c#######################################################################
c
      open (unit=sfunit, file=sndfile,status='old',
     :            iostat=istat)

      IF(istat.ne.0) THEN
            write (6,1000) 'Error in opening sounding file'
            write (6,'(1x,a,a)') 'Check file ',sndfile
            STOP
      ELSE 
        write(6,'(1x,/2a,a,i3/)') 'Sounding file ', sndfile,
     :          ' opened using FORTRAN unit',sfunit
      ENDIF

      read (sfunit,1000,iostat=istat) dummy ! header line
      write(6,1000)dummy

      read (sfunit,1000,iostat=istat) dummy ! miscellaneous descriptions
      write(6,1000)dummy

      read (sfunit,1000,iostat=istat) stime ! time of sounding
      write(6,1000)stime

      read (sfunit,1000,iostat=istat) sdate ! date of sounding
      write(6,1000)sdate

      read (sfunit,1000,iostat=istat) sloc  ! location of sounding
      write(6,1000)sloc

      read (sfunit,*,iostat=istat) height,therm,humid ! data type flags
      write(6,'(3(1x,a))') height, therm, humid

      read (sfunit,*,iostat=istat) zsnd1,psnd1  
      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 .gt.lvlprof) 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'

        STOP

      ENDIF


      read (sfunit, 1000, iostat=istat) dummy ! a seperator line,
           ! possibly descriptive of the file contents for human readers
      write(6,1000)dummy

      IF(istat.ne.0) THEN

        write (6,1000) 'Error in sounding file'
        write (6,1000) 'Error in sounding file header'
        write (6,1000) 'Check sounding file', sndfile

        STOP

      ENDIF

c#######################################################################
c
c     Read in the sounding data:
c
c#######################################################################

      DO 10 K = lvlsnd ,1,-1 

        read (sfunit, *, iostat=istat)
     1        zsnd(k), ptsnd(k), qvsnd(k), usnd(k), vsnd(k)

        write(6,'(1x,i4,5f13.6)')
     :  k,zsnd(k), ptsnd(k), qvsnd(k), usnd(k), vsnd(k)

        IF(istat.ne.0) THEN

          write (6,1000) 'Error in sounding file'
          write (6,1000) 'Error in sounding file data'
          write (6,1000) 'Check ',sndfile
          STOP

        ENDIF

10    CONTINUE

      close (unit=sfunit)      ! close the sounding file

      write(6,'(/1x,a)')'Sounding data reading complete.'
      write(6,'(1x,a,i6,a/)')'FORTRAN unit',sfunit,' closed.'
c
c#######################################################################
c
c     End of sounding data reading.
c
c#######################################################################
c

c
c#######################################################################
c
c     If input is not in height, copy it into pressure array and
c     set the input data type flag.
c
c#######################################################################
c
      IF( height(1:1).ne.'h' .and. height(1:1).ne.'H' ) THEN

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

      ELSE

        in_is_p = .false.

      ENDIF
c
c#######################################################################
c
c     If input is in temperature, copy it from potential temperature 
c     array and set the input data type flag.
c
c#######################################################################
c
      IF( therm(1:1).eq.'t' .or. therm(1:1).eq.'T' ) THEN

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

      ELSE

        in_is_t = .false.

      ENDIF
c
c#######################################################################
c
c     If input is in relative humidity, copy it from specific humidity  
c     array and set the input data type flag.
c
c#######################################################################
c
      IF( humid(1:1).eq.'r' .or. humid(1:1).eq.'R' ) THEN

        in_is_rh = .true.
        DO 3001 k=1,lvlsnd 
          rhsnd(k) =qvsnd(k)
3001    CONTINUE

      ELSE

        in_is_rh = .false.
        DO 3002 k=1,lvlsnd 
          rhsnd(k) = 0.0
3002    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Now calculate the derived variables from the input sounding
c
c#######################################################################
c

      IF( in_is_p ) THEN 
c
c#######################################################################
c
c     If sounding data are specified on pressure levels
c
c#######################################################################
c
        DO 3045 k=1,lvlsnd

          IF( psnd(k).le.0.0) THEN

            write(6,'(2(/1x,a))')
     :           'Negative pressure found in the sounding profile.',
     :           'Job stopped. Check subroutine INIBASE.'
            STOP

          ENDIF

3045    CONTINUE

        IF( in_is_t ) THEN 

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

        ELSE

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

        ENDIF

        IF( in_is_rh ) THEN 

          CALL getqvs( psnd, tsnd, lvlsnd, qvsnd )

          DO 3052 k=1,lvlsnd
            qvsnd(k) = qvsnd(k) * rhsnd(k)
3052      CONTINUE

        ENDIF
c
c#######################################################################
c
c     The virtual temperature Tv:
c
c#######################################################################
c
        DO 3055 k = 1,lvlsnd
          ptvsnd(k) = tsnd(k)*(1.0+rvdrd*qvsnd(k))/(1.0+qvsnd(k))
3055    CONTINUE
c 
c#######################################################################
c
c     Integrate hydrostatic relation dz/dp = - RTv/(pg) for z
c
c#######################################################################
c
        zsnd(1) = zsnd1

        DO 3060 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))
3060    CONTINUE
c

      ELSEIF( .not. in_is_p ) THEN 

c
c#######################################################################
c
c     If sounding data are specified on height levels
c
c#######################################################################
c
        IF( in_is_rh ) THEN  
c
c#######################################################################
c
c     If inputs are height and relative humidity, the pressure,
c     temperature and specific humidity profiles are coupled, so
c     interation is needed. We use a predetermined number (10) 
c     of iterations, which should be enough for convergence. 
c
c#######################################################################
c
          niter = 10

          DO 4001 k=1,lvlsnd 
            qvsnd(k) =0.0
4001      CONTINUE
        ELSE
          niter = 1
        ENDIF


        DO 1010 iter = 1, niter

          IF( .not. in_is_t) THEN 
c 
c#######################################################################
c
c     Temperature data are given as potential temperature 
c
c#######################################################################
c

c#######################################################################
c
c     Calculate ptvsnd, the virtual potential temperature sounding array
c
c#######################################################################
c
            DO 20 k = 1,lvlsnd
              ptvsnd(k) = ptsnd(k)*(1.0+rvdrd*qvsnd(k))/(1.0+qvsnd(k))
20          CONTINUE
c
c#######################################################################
c
c     Integrate the hydrostatic relation to get PISND, the 
c     nondimensional pressure: d(pi)/dz = -g/(cp * ptv)
c
c#######################################################################
c
c

            pisnd(1) = (psnd1/p0)**(rddcp) 

            DO 30 k = 2,lvlsnd
              pisnd(k) = pisnd(k-1) - g/cp  
     :                   *(zsnd(k)-zsnd(k-1)) 
     :                   /(0.5*(ptvsnd(k)+ptvsnd(k-1)))
30          CONTINUE

            DO 40 k = 1,lvlsnd
              psnd(k) = p0 * (pisnd(k)**(cp/rd))
              tsnd(k) = ptsnd(k) * pisnd(k)
40          CONTINUE

          ELSE
c
c#######################################################################
c
c     Temperature data are given as temperature 
c
c#######################################################################

c#######################################################################
c
c     Calculate the virtual temperature and store it in ptvsnd
c
c#######################################################################
c
            DO 60 k = 1,lvlsnd
              ptvsnd(k) = tsnd(k)*(1.0+rvdrd*qvsnd(k))/(1.0+qvsnd(k))
60          CONTINUE
c
c#######################################################################
c
c     Calculate log(p) according to  d(log(p))/dz = -g/(r * tv)
c
c#######################################################################
c
            pisnd(1) = log( psnd1 )

            DO 70 k = 2,lvlsnd
              pisnd(k) = pisnd(k-1) - g/rd  
     :                   *(zsnd(k)-zsnd(k-1))
     :                   /(0.5*(ptvsnd(k)+ptvsnd(k-1)))
70          CONTINUE

            DO 80 k = 1,lvlsnd
              psnd(k) = exp( pisnd(k) ) 
              ptsnd(k) = tsnd(k) * (p0/psnd(k))** (rddcp)
80          CONTINUE

          ENDIF

          IF( in_is_rh ) THEN

            CALL getqvs( psnd, tsnd, lvlsnd, qvsnd )

            DO 1003 k=1,lvlsnd
              qvsnd(k) = qvsnd(k) * rhsnd(k)
1003        CONTINUE

          ENDIF

1010    CONTINUE

      ENDIF

      IF( .not. in_is_rh ) THEN

        CALL getqvs( psnd, tsnd, lvlsnd, rhsnd )

        DO 1005 k=1,lvlsnd
          rhsnd(k) = qvsnd(k)/rhsnd(k)
1005    CONTINUE

      ENDIF
c
c#######################################################################
c
c     Write out the uninterpolated sounding profiles
c
c#######################################################################
c
      write(6,'(/1x,a/)') 'Uninterpolated sounding profiles:'
      write(6,'(1x,2a)') 
     :     '   n    z(m)   p(Pascal)  PT(K)   T(K)',
     :     '    Qv(kg/kg)   RH      U(m/s)   V(m/s)'

      DO 1100 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)
1100  CONTINUE

      RETURN

1000  format(1x,a)

      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE SNDINTRP                   ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE SNDINTRP(uprof,vprof,ptprof,pprof,qvprof,zsnd,lvlsnd,
     :              zprof,lvlprof, temp1)

c
c#######################################################################
c
c     PURPOSE:
c
c     Interpolate the input sounding profiles on the grid "zsnd" to profiles 
c     on the regular grid "zprof". 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     3/17/1991.
c
c     MODIFICATION HISTORY:
c
c     5/02/92 (M. Xue)
c     Added full documentation.
C
c     5/03/92 (M. Xue)
c     Further documentation. 
C
c#######################################################################
c
c     INPUT :
c
c       uprof    u-velocity sounding data at original data levels (m/s)
c       vprof    v-velocity sounding data at original data levels (m/s)
c       ptprof   Potential temperature sounding data at original data
c                levels (K)
c       pprof    Pressure sounding data at original data levels (Pascal)
c       qvprof   Water vapor specific humidity sounding data at the
c                original data levels (kg/kg)
c       zsnd     The height of the input sounding data levels (m)
c       lvlsnd   The number of the input sounding data
c
c       zprof    The height of equally spaced grid levels to which 
c                the sounding data are interpolated (m).
c       lvlprof  The number of such grid levels.
c
c     OUTPUT: 
c
c       uprof    u-velocity data defined on grid levels given in zprof
c       vprof    v-velocity data defined on grid levels given in zprof
c       ptprof   Potential temperature data defined on grid levels 
c                given in zprof
c       pprof    Pressure data defined on grid levels given in zprof
c       qvprof   Water vapor specific humidity data defined on grid  
c                levels given in zprof
c
c     WORK ARRAY:
c
c       temp1    A work array
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer lvlprof, lvlsnd

      real uprof (lvlprof)
      real vprof (lvlprof)
      real ptprof(lvlprof)
      real pprof (lvlprof)
      real qvprof(lvlprof)
      real zprof (lvlprof)
      real zsnd  (lvlsnd)
c
      real temp1 (lvlsnd)      ! Temporary array
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer k
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external inte1d
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c


      DO 10 k=1,lvlsnd 
        temp1(k) = uprof(k)
10    CONTINUE

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

      DO 20 k=1,lvlsnd 
        temp1(k) = vprof(k)
20    CONTINUE

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

      DO 30 k=1,lvlsnd 
        temp1(k) = ptprof(k)
30    CONTINUE

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

      DO 40 k=1,lvlsnd 
        temp1(k) = pprof(k)
40    CONTINUE

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

      DO 50 k=1,lvlsnd 
        temp1(k) = qvprof(k)
50    CONTINUE

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

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE GETQVS                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE GETQVS( p, t, n, qvs )
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the saturation specific humidity using Teten's formula.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     3/17/1991.
c
c     MODIFICATION HISTORY:
c
c     5/02/92 (M. Xue)
c     Added full documentation.
C
c     5/03/92 (M. Xue)
c     Further documentation. 
C
c#######################################################################
c
c     INPUT :
c
c       p        Pressure (Pascal)
c       t        Temperature (K)
c       n        The dimension of input arrays.
c
c     OUTPUT:
c
c       qvs      Saturation water vapor specific humidity (kg/kg).
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer n

      real p  (n)
      real t  (n)
      real qvs(n)

c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      DO 10 i=1,n
        qvs(i) = (380. / p(i)) *
     :  exp( 17.27 * (t(i)-273.15)/(t(i) - 35.86) )
10    CONTINUE

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE INTE1D                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE INTE1D(fdat,zdat,ndat,fpro,zpro,npro)

c
c#######################################################################
c
c     PURPOSE:
c
C     Interpolates from sounding data to a vertical profile with nn.form 
c     grid spaces.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     3/17/1991.
c
c     MODIFICATION HISTORY:
c
c     5/02/92 (M. Xue)
c     Added full documentation.
C
c     5/03/92 (M. Xue)
c     Further documentation. 
C
c#######################################################################
c
c     Input :
c
c       fdat     1-d array defined at levels in zdat to be 
c                interpolated to levels defined by zpro.
c       zdat     The height of input data in fdat.
c       ndat     The number of data in input array fdat.
c
c       zpro     The height of grid levels to which data are interpolated 
c       npro     The number of interpolated data levels.
c
c     Output:
c
c       fpro     Output data interpolated to grid levels defined by zpro.
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer ndat,npro

      real fdat(ndat)
      real zdat(ndat)
      real zpro(npro)
      real fpro(npro)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i0,i,j

c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      
      i0=1
      DO 20 J=1,NPRO

        IF(zpro(j).le.zdat(1))THEN
          fpro(j)=fdat(1)
          GOTO 20
        ENDIF

        IF(zpro(j).ge.zdat(ndat))THEN
          fpro(j)=fdat(ndat)
          GOTO 20
        ENDIF

c
c#######################################################################
c
c     Find index for interpolation
c
c#######################################################################
c
        DO 10 i=i0,ndat-1
          IF(zpro(j).ge.zdat(i).and.zpro(j).lt.zdat(i+1))GOTO 15
10      CONTINUE

        GOTO 20
15      i0=i

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

20    CONTINUE

      RETURN
      END
