C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    PROFIN      INGEST OF PROFILER DATA
C   PRGMMR: TRACY SMITH      ORG: FSL/PROFS  DATE: 90-07-19
C
C ABSTRACT: DRIVER  FOR THE INGEST OF PROFILER DATA.
C
C PROGRAM HISTORY LOG:
C   84-??-??  T.L. SMITH    ORIGINAL VERSION FROM R. LIEBERMAN.
C   86-12-29  T.L. SMITH    CHANGED ROUTINE TO GET WINDS TO
C                           QCD_PRF_WINDS (K. BREWSTER ROUTINE)
C   87-02-02  T.L. SMITH    THETA VERSION
C   87-07-13  S.G. BENJAMIN CHANGED TO ASSIGN PRESSURE TO WIND
C                           PROFILER DATA FROM NGM DATA SINCE
C                           ANALYSIS AND QC NEED BOGUS PRES OB.
C   87-09-22  S.G.BENJAMIN  ASSIGN "CHANGED" QC BIT TO ALL LEVELS
C                           BELOW GROUND.
C   87-10-23  S.G.BENJAMIN  PRESSURE OBS SET TO NMCBOGS DON'T ALSO
C                           GET SET TO VERTCHG.
C   88-05-26  T.L. SMITH    CHANGED FOR NEWNAT AND GET_CONFIG
C
C   88-08-31  T.L. SMITH    ADDED NEW DATABASE VARIABLES,
C                           SUBSRC,RPRTTM,PLATID
C
C   89-07-14  T.L. SMITH    TURNED OFF STAPL RADIOMETER DATA,
C                           BY SETTING STAPLE TO FALSE UFN.
C   90-02-19  S.G. BENJAMIN SET PC = SPVAL FOR ALL PROFILER
C   90-11-05  T.L. Smith    Added ability to read network
C                                profilers
C**********************************************************************
C
C USAGE:    CALL PROFIN(ATIME,PMAPS,ZMAPS,g3,ISTA_COUNT,O,QC,OID
C      +                OI,OJ,OTYPE,OTIME,STATUS)
C   INPUT ARGUMENT LIST:
C     ATIME     - CHARACTER   ANALYSIS TIME
C     PMAPS     - REAL ARRAY  PRESSURE FROM PREVIOUS MAPS ANALYSIS.
C     ZMAPS     - REAL ARRAY  HEIGHTS      "
C
C   OUTPUT ARGUMENT LIST:
C     ISTA_COUNT- INTEGER     NUMBER OF OBSERVATIONS
C     O         - REAL ARRAY  OBSERVATIONS (P,Z,THETA,PC,U,V)
C     QC        - INTEGER ARRAY  QUALITY CONTROL INFORMATION
C     OID       - CHARACTER ARRAY STATION IDENTIFICATION
C     OI        - REAL ARRAY  X LOCATION OF OB ON MAPS GRID
C     OJ        - REAL ARRAY  Y LOCATION OF OB ON MAPS GRID
C     OTYPE     - INTEGER ARRAY DATA TYPE OF OBSERVATION
C     OTIME     - INTEGER ARRAY TIME OF OBSERVATION-ANALYSIS TIME
C     STATUS    - INTEGER  RETURN STATUS
C
C   INPUT FILES:
C     FT03F001 - PROFILER NMC DATA FILE.
C
C   OUTPUT FILES:
C     MPSPROF  - MAPS FORMAT PROFILER DATA FILE.
C     FT06F001 - DIAGNOSTIC MESSAGE.
C
C REMARKS:
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN-77
C   MACHINE:  NAS-9000
C$$$
      SUBROUTINE PROFIN(ATIME,PMAPS,ZMAPS,g3,ISTA_COUNT,O,QC,
     +          OID,OCLD, OI,OJ,OTYPE,OTIME,STATUS)

C**********************************************************************
C*     DRIVER FOR INGEST OF PROFILER DATA
C**********************************************************************

      INCLUDE 'IMPLICIT'
      INCLUDE 'MAPSCON'
      INCLUDE 'METCON'

      INTEGER MAXSTA,MAXLEV
      PARAMETER (MAXSTA=500,MAXLEV=300)

      REAL O(NVAR_P,NSTA_P),
     +     OI(NSTA_P),OJ(NSTA_P)

      CHARACTER*(LSTAID_P) OID(NSTA_P)

      INTEGER QC(NVAR_P,NSTA_P),
     +        OTYPE(NSTA_P),
     +        OTIME(NSTA_P),
     +        ISTA_COUNT,
     +        I4TIME,
     +        NLVLS(MAXSTA),
     +        NVLEV(MAXSTA),
     +        NTVLEV(MAXSTA)
      INTEGER STANUM(MAXSTA),
     +        VSTANUM(MAXSTA),
     +        RSTANUM(MAXSTA),
     +        NUM_OF_STAT,
     +        NUM_PROF,
     +        NUM_VAD,
     +        NUM_OF_RASS,
     +        LUNIT,
     +        I,J,K,L,IRASS,KNT,
     +        NREJLST,
     +        STATUS,ISTAT

      INTEGER      OCLD  (NVARCLD_P,mxst_P)

      REAL    P(MAXLEV),
     +        thv(maxlev),
     +        ZW(MAXLEV,MAXSTA),
     +        UU(MAXLEV,MAXSTA),
     +        VV(MAXLEV,MAXSTA),
     +        VZW(MAXLEV,MAXSTA),
     +        VUU(MAXLEV,MAXSTA),
     +        VVV(MAXLEV,MAXSTA),
     +        RMS(MAXLEV,MAXSTA),
     +        TV(MAXLEV,MAXSTA),
     +        TV_HGT(MAXLEV,MAXSTA)
      REAL    UTEMP,
     +        VTEMP,
     +        SINDL,
     +        COSDL,
     +        LAT(MAXSTA),
     +        LON(MAXSTA),
     +        VLAT(MAXSTA),
     +        VLON(MAXSTA),
     +        RLAT(MAXSTA),
     +        RLON(MAXSTA)
      REAL    XLIMN,
     +        XLIMX,
     +        YLIMN,
     +        YLIMX,
     +        XCART,
     +        YCART,
     +        ELEVATION(MAXSTA),
     +        VELEV(MAXSTA),
     +        RELEVATION(MAXSTA)
C --- The following variables are defined for the MAPS analysis levels

C** P/Z and conversion variables
      REAL    ZFG(NZTN_P),PFG(NZTN_P),thvFG(NZTN_P),
     +        FNEWT,TRBL1,dzlimit
      real    table1(nztn_p,nztn_p)
      real    table2(nztn_p,nztn_p)
      real    um(nztn_p)
      real    vm(nztn_p)

      INTEGER II,JJ,IP1,JP1,IP2,JP2,IP3,JP3,ivar,nlev,k1

      INTEGER RJUNIT
      REAL    PMAPS  (NX_P,NY_P,NZTN_P),
     +        ZMAPS  (NX_P,NY_P,NZTN_P)
     +       ,g3     (NX_P,NY_P,NZTN_P,nvar3_p)
      REAL    XC,YC,XX,YY,
     +        XI,XI1,XI2,XI3,
     +        YJ,YJ1,YJ2,YJ3,
     +        FACTOR1,FACTOR2,FACTOR3,FACTOR4,
     +        FACTOR5,FACTOR6,FACTOR7,FACTOR8

      CHARACTER*9          ATIME
      CHARACTER*9          STANAM(MAXSTA),
     +                     VSTANAM(MAXSTA),
     +                     RSTANAM(MAXSTA)

      CHARACTER*5          REJLST(MAXSTA)

C*** start program
      XLIMN=X_MIN_P - DX_P/2.
      XLIMX=X_MAX_P + DX_P/2.
      YLIMN=Y_MIN_P - DY_P/2.
      YLIMX=Y_MAX_P + DY_P/2.

      CALL GETLUN(LUNIT,1)
      OPEN(UNIT = LUNIT,
     +     FILE = 'PROFDAT',
     +     STATUS = 'UNKNOWN',
     +     ACCESS = 'SEQUENTIAL',
     +     FORM = 'FORMATTED')

C*** Get the time in integer format.
      CALL MKI4TM(ATIME,
     +            I4TIME,
     +            STATUS)
      IF(STATUS.NE.SUCCESS_P) THEN
        WRITE(LUNIT,*) 'BAD STATUS FROM MKFNAM, PROFIN',STATUS
        RETURN
      ENDIF

      ISTA_COUNT = 0

      do i=1,nvar_p
      do j=1,maxsta
        qc(i,j)=0
        O (i,j)=spval_p
      enddo
      enddo

      WRITE(LUNIT,112) ATIME
112   FORMAT(1X,'Profiler DATA FOR ',A9)


C*** First, get all wind profiler data
      CALL GTPROF(ATIME,MAXSTA,MAXLEV,
     +            NUM_PROF,STANUM,STANAM,
     +            ZW,UU,VV,NLVLS,
     +            ELEVATION,LAT,LON,STATUS)
      IF (STATUS.NE.SUCCESS_P) THEN
        WRITE(LUNIT,*) 'BAD STATUS FROM GTPROF, STAT=',STATUS
        NUM_PROF=0
      ENDIF

      NUM_OF_STAT=NUM_PROF

C*** next, get all radar derived wind data
      CALL GETVAD(ATIME,MAXSTA,MAXLEV,NUM_VAD,
     +            VSTANUM,VSTANAM,
     +            VZW,VUU,
     +            VVV,NVLEV,
     +            VELEV,VLAT,
     +            VLON,RMS,STATUS)
      IF (STATUS.NE.SUCCESS_P) THEN
        WRITE(LUNIT,*) 'BAD STATUS FROM GETVAD, STAT=',STATUS
        NUM_VAD=0
      ENDIF

      NUM_OF_STAT=NUM_OF_STAT+NUM_VAD

      write(6,*)'                                   '
      write(6,*)'***********************************'
      write(6,*)'No. of stations from profiler/VAD ',
     +    num_prof, num_vad
      write(6,*)'***********************************'
C*** Combine the profiler and VAD winds
      DO I=1,NUM_VAD
        KNT=I+NUM_PROF
        STANUM(KNT)=VSTANUM(I)
        STANAM(KNT)=VSTANAM(I)
        NLVLS(KNT)=NVLEV(I)
        ELEVATION(KNT)=VELEV(I)
        LAT(KNT)=VLAT(I)
        LON(KNT)=VLON(I)
        DO J=1,NVLEV(I)
          ZW(J,KNT)=VZW(J,I)
          UU(J,KNT)=VUU(J,I)
          VV(J,KNT)=VVV(J,I)
        ENDDO
      ENDDO

C*** finally, get RASS temperatures to fill in later
      CALL GTRASS(ATIME,MAXSTA,MAXLEV,
     +            NUM_OF_RASS,TV,TV_HGT,NTVLEV,
     +            RSTANUM,RSTANAM,RLAT,RLON,RELEVATION,STATUS)
      IF (STATUS.NE.SUCCESS_P) THEN
        WRITE(LUNIT,*) 'BAD STATUS FROM GTRASS, STAT=',STATUS
        WRITE(LUNIT,*) 'NO RASS DATA IN PROFILER TEMPS FOR ',ATIME
        NUM_OF_RASS=0
      ENDIF

      IF(NUM_OF_STAT.LT.1) THEN
         STATUS=FAILURE_P
         GOTO 999
      ELSE
         STATUS=SUCCESS_P
      ENDIF
C** RETURN TO MAIN ROUTINE GOTO

C --- Read in reject listed profiler sites
      CALL GETLUN(RJUNIT,1)
      NREJLST=0
      OPEN(UNIT = RJUNIT,
     +     FILE = 'REJCTPRF',
     +     STATUS = 'OLD',
     +     ACCESS = 'SEQUENTIAL',
     +     FORM = 'FORMATTED')

      DO  10 I=1,100
        READ (RJUNIT,887) REJLST(I)
C        WRITE(LUNIT,*) REJLST(I),I
        IF(REJLST(I).EQ.'XXXXX') GOTO 886
        NREJLST=NREJLST+1
10    CONTINUE

886   CLOSE(RJUNIT)
      CALL GETLUN(RJUNIT,2)
887   FORMAT(A5)

      write(6,*)'                                   '
      write(6,*)'***********************************'
      write(6,*)'No. of ob levels before profiler/VAD ',
     +    ista_count
      write(6,*)'***********************************'
      
C*** Start Station Loop
      DO 100  I=1,NUM_OF_STAT

C*** First check to see if this station is blacklisted
        DO 15 J=1,NREJLST
          IF (REJLST(J) .EQ. STANAM(I)(1:5)) THEN
C           WRITE(LUNIT,*) 'FOUND REJLST',STANAM(I)
            GO TO 100
          ENDIF
15      CONTINUE

        IF(NLVLS(I).LE.1) GOTO 100

C*sb Check for lat/lon range
        IF (LAT(I).GT.90. .OR. LAT(I).LE.0.) GO TO 100
        IF (LON(I).GE.180..OR. LON(I).LE.-180.) GO TO 100

C*** Is station in MAPS domain
        CALL W3FB11(LAT(I),360.+LON(I),LAT_LL_P,360.+LON_LL_P,DX_P,
     +              360.+LON_XX_P,LAT_TAN_P,XC,YC)
        XCART=(XC-1)*DX_P
        YCART=(YC-1)*DY_P
        IF     (XCART .LT. XLIMN  .OR.  XCART .GT. XLIMX
     +    .OR.  YCART .LT. YLIMN  .OR.  YCART .GT. YLIMX) then
          write(6,*)'Profiler/VAD station outside of domain ',stanam(i),
     +        lat(i),lon(i)
          GO TO 100
        end if

        WRITE(LUNIT,*)' STATION,ELEV,LAT,LON ',StaNam(I),StaNum(I),
     +          elevation(I),lat(I),lon(I)


C***Get conversion grids for p to z for this location
        II=MAX(1,INT(XC)-1)
        JJ=MAX(1,INT(YC)-1)
        XX=FLOAT(II)
        YY=FLOAT(JJ)

        IP1 = MIN(II+1,NX_P)
        JP1 = MIN(JJ+1,NY_P)
        IP2 = MIN(II+2,NX_P)
        JP2 = MIN(JJ+2,NY_P)
        IP3 = MIN(II+3,NX_P)
        JP3 = MIN(JJ+3,NY_P)

        XI = XC-XX
        XI1 = XC-XX-1.
        XI2 = XC-XX-2.
        XI3 = XC-XX-3.

        YJ = YC-YY
        YJ1 = YC-YY-1.
        YJ2 = YC-YY-2.
        YJ3 = YC-YY-3.

        FACTOR1 = XI1*XI2*XI3/(-6.)
        FACTOR2 = XI*XI2*XI3/(2.)
        FACTOR3 = XI*XI1*XI3/(-2.)
        FACTOR4 = XI*XI1*XI2/(6.)
        FACTOR5 = YJ1*YJ2*YJ3/(-6.)
        FACTOR6 = YJ*YJ2*YJ3/(2.)
        FACTOR7 = YJ*YJ1*YJ3/(-2.)
        FACTOR8 = YJ*YJ1*YJ2/(6.)

        DO  90 K=1,NZTN_P
          PFG(K) =
     +             FACTOR1 * FACTOR5 * PMAPS(II ,JJ ,K)
     +           + FACTOR2 * FACTOR5 * PMAPS(IP1,JJ ,K)
     +           + FACTOR3 * FACTOR5 * PMAPS(IP2,JJ ,K)
     +           + FACTOR4 * FACTOR5 * PMAPS(IP3,JJ ,K)
     +           + FACTOR1 * FACTOR6 * PMAPS(II ,JP1,K)
     +           + FACTOR2 * FACTOR6 * PMAPS(IP1,JP1,K)
     +           + FACTOR3 * FACTOR6 * PMAPS(IP2,JP1,K)
     +           + FACTOR4 * FACTOR6 * PMAPS(IP3,JP1,K)
     +           + FACTOR1 * FACTOR7 * PMAPS(II ,JP2,K)
     +           + FACTOR2 * FACTOR7 * PMAPS(IP1,JP2,K)
     +           + FACTOR3 * FACTOR7 * PMAPS(IP2,JP2,K)
     +           + FACTOR4 * FACTOR7 * PMAPS(IP3,JP2,K)
     +           + FACTOR1 * FACTOR8 * PMAPS(II ,JP3,K)
     +           + FACTOR2 * FACTOR8 * PMAPS(IP1,JP3,K)
     +           + FACTOR3 * FACTOR8 * PMAPS(IP2,JP3,K)
     +           + FACTOR4 * FACTOR8 * PMAPS(IP3,JP3,K)

          ZFG(K) =
     +               FACTOR1 * FACTOR5 * ZMAPS(II ,JJ ,K)
     +             + FACTOR2 * FACTOR5 * ZMAPS(IP1,JJ ,K)
     +             + FACTOR3 * FACTOR5 * ZMAPS(IP2,JJ ,K)
     +             + FACTOR4 * FACTOR5 * ZMAPS(IP3,JJ ,K)
     +             + FACTOR1 * FACTOR6 * ZMAPS(II ,JP1,K)
     +             + FACTOR2 * FACTOR6 * ZMAPS(IP1,JP1,K)
     +             + FACTOR3 * FACTOR6 * ZMAPS(IP2,JP1,K)
     +             + FACTOR4 * FACTOR6 * ZMAPS(IP3,JP1,K)
     +             + FACTOR1 * FACTOR7 * ZMAPS(II ,JP2,K)
     +             + FACTOR2 * FACTOR7 * ZMAPS(IP1,JP2,K)
     +             + FACTOR3 * FACTOR7 * ZMAPS(IP2,JP2,K)
     +             + FACTOR4 * FACTOR7 * ZMAPS(IP3,JP2,K)
     +             + FACTOR1 * FACTOR8 * ZMAPS(II ,JP3,K)
     +             + FACTOR2 * FACTOR8 * ZMAPS(IP1,JP3,K)
     +             + FACTOR3 * FACTOR8 * ZMAPS(IP2,JP3,K)
     +             + FACTOR4 * FACTOR8 * ZMAPS(IP3,JP3,K)

          thvFG(K) =
     +               FACTOR1 * FACTOR5 * g3(II ,JJ ,K,theta_p)
     +             + FACTOR2 * FACTOR5 * g3(IP1,JJ ,K,theta_p)
     +             + FACTOR3 * FACTOR5 * g3(IP2,JJ ,K,theta_p)
     +             + FACTOR4 * FACTOR5 * g3(IP3,JJ ,K,theta_p)
     +             + FACTOR1 * FACTOR6 * g3(II ,JP1,K,theta_p)
     +             + FACTOR2 * FACTOR6 * g3(IP1,JP1,K,theta_p)
     +             + FACTOR3 * FACTOR6 * g3(IP2,JP1,K,theta_p)
     +             + FACTOR4 * FACTOR6 * g3(IP3,JP1,K,theta_p)
     +             + FACTOR1 * FACTOR7 * g3(II ,JP2,K,theta_p)
     +             + FACTOR2 * FACTOR7 * g3(IP1,JP2,K,theta_p)
     +             + FACTOR3 * FACTOR7 * g3(IP2,JP2,K,theta_p)
     +             + FACTOR4 * FACTOR7 * g3(IP3,JP2,K,theta_p)
     +             + FACTOR1 * FACTOR8 * g3(II ,JP3,K,theta_p)
     +             + FACTOR2 * FACTOR8 * g3(IP1,JP3,K,theta_p)
     +             + FACTOR3 * FACTOR8 * g3(IP2,JP3,K,theta_p)
     +             + FACTOR4 * FACTOR8 * g3(IP3,JP3,K,theta_p)

90      CONTINUE

        CALL DTABLE(ZFG,PFG,  TABLE1,NZTN_P,1,TRBL1,NZTN_P)
        CALL DTABLE(ZFG,thvFG,TABLE2,NZTN_P,1,TRBL1,NZTN_P)

C*** Get factors for rotating winds from lat-long to map coord
            CALL ROTFAC(LON(I),SINDL,COSDL,LON_XX_P)

C     ------------------------------------------------------------
C ----------  First, call winder to get winds interpolated
C             to background levels.
C     ------------------------------------------------------------

         if (i.lt.num_prof ) then
           dzlimit = 1200.
         else
           dzlimit = 1000.
         end if
  
         do k = 1,nztn_p
           um(k) = spval_p
           vm(k) = spval_p
         end do
         
         nlev = nlvls(i)


C --- Interpolate u and v (wind components) linearly w.r.t. height
C     ------------------------------------------------------------
         CALL WINDER(Uu(1,i),Vv(1,i),Zw(1,i), Um,Vm,Zfg,nlev)

         do k = 2,nlvls(i)
           if (zw(k,i)-zw(k-1,i) .gt. dzlimit) then
             do k1 = 1,nztn_p
               if (zfg(k1).gt.zw(k-1,i) .and.
     1             zfg(k1).lt.zw(k  ,i) ) then
                 um(k1) = spval_p
                 vm(k1) = spval_p
               end if
             end do
           end if
         end do

         do 987 k = 1,nztn_p
           if (um(k).lt.spval_p) then
        
C*** OK, so we have a good ob from this station, count it
            ISTA_COUNT=ISTA_COUNT+1

C*** fill station with SPVAL_Ps before proceeding
            O(P_P,ISTA_COUNT)=SPVAL_P
            O(H_P,ISTA_COUNT)=SPVAL_P
            O(THETA_P,ISTA_COUNT)=SPVAL_P
            O(PC_P,ISTA_COUNT)=SPVAL_P
            O(U_P,ISTA_COUNT)=SPVAL_P
            O(V_P,ISTA_COUNT)=SPVAL_P
            O(PWOB_P,ISTA_COUNT)=SPVAL_P

C***  ...And rotate winds
            UTEMP=Um(K)
            VTEMP=Vm(K)
            O(U_P,ISTA_COUNT) = COSDL*UTEMP + SINDL*VTEMP
            O(V_P,ISTA_COUNT) =-SINDL*UTEMP + COSDL*VTEMP
          ELSE
            GOTO 987
C** Continue on to next level
          ENDIF


C* Set QC to bogus for derived pressures, and actual heights, so
C* as not to confuse the analysis with the p/z conversion.

            O(P_P,ISTA_COUNT)=Pfg(K)
            CALL SETQC(QC(P_P,Ista_count),MSTRBAD_P,
     +                  QC(P_P,ista_count),ISTAT)
            CALL SETQC(QC(P_P,ista_count),BOGSNMC_P,
     +                 QC(P_P,ista_count),ISTAT)
            CALL SETQC(QC(h_P,ista_count),BOGSNMC_P,
     +                 QC(h_P,ista_count),ISTAT)
            CALL SETQC(QC(H_P,Ista_count),MSTRBAD_P,
     +                 QC(H_P,Ista_count),ISTAT)

          IF (thvfg(K).LT.SPVAL_P) THEN
            O(theta_P,ISTA_COUNT)=thvfg(K)
            CALL SETQC(QC(theta_P,Ista_count),MSTRBAD_P,
     +                  QC(theta_P,ista_count),ISTAT)
            CALL SETQC(QC(theta_P,ista_count),BOGSNMC_P,
     +                 QC(theta_P,ista_count),ISTAT)
          ELSE
            write(6,*) 'Something wrong with PROFIN thv/z,',
     +              thvfg(K),Zfg(K)
            ista_count=ista_count-1
            goto 987
          ENDIF

          O(H_P,ISTA_COUNT)=Zfg(K)


          OTYPE(ISTA_COUNT) = 203
C                        ! network data
          IF(I.GT.NUM_PROF) THEN
            OTYPE(ISTA_COUNT) = 210
C*** subsrc 210 represents VAD winds
          ENDIF


          OID(ista_count) = STANAM(I)
          OI(ista_count) = XC
          OJ(ista_count) = YC
          OTIME(ISTA_COUNT) = 0
          WRITE(LUNIT,535) (O(ivar,ISTA_COUNT), ivar=1,NVAR_P)

987     continue

967     continue


C     ------------------------------------------------------------
C       Now, process original wind data levels
C     ------------------------------------------------------------

C*** Fill in from 1st wind gate up
        DO  80 K=1,NLVLS(I)
          IF (UU(K,I) .LT. SPVAL_P  .AND.  VV(K,I) .LT. SPVAL_P) THEN
C***   Skip level if VAD and above 6200 m.
          IF(I.GT.NUM_PROF .and. zw(k,i).gt.6200.) go to 80
C***   Skip level if both u and v are equal exactly to zero.

          if (uu(k,i).eq.0. .and. vv(k,i).eq.0.) go to 80

C*** OK, so we have a good ob from this station, count it
            ISTA_COUNT=ISTA_COUNT+1

C*** fill station with SPVAL_Ps before proceeding
            O(P_P,ISTA_COUNT)=SPVAL_P
            O(H_P,ISTA_COUNT)=SPVAL_P
            O(THETA_P,ISTA_COUNT)=SPVAL_P
            O(PC_P,ISTA_COUNT)=SPVAL_P
            O(U_P,ISTA_COUNT)=SPVAL_P
            O(V_P,ISTA_COUNT)=SPVAL_P
            O(PWOB_P,ISTA_COUNT)=SPVAL_P

C***  ...And rotate winds
            UTEMP=UU(K,I)
            VTEMP=VV(K,I)
            O(U_P,ISTA_COUNT) = COSDL*UTEMP + SINDL*VTEMP
            O(V_P,ISTA_COUNT) =-SINDL*UTEMP + COSDL*VTEMP
          ELSE
            GOTO 80
C** Continue on to next level
          ENDIF

C*** Assign derived pressure to P and heights to Z
          P(K)=
     +    FNEWT(ZFG,PFG,TABLE1,NZTN_P,1,1,ZW(K,I),TRBL1,NZTN_P)
          IF (P(K).LE.0.0) THEN
            P(K)=SPVAL_P+1
          ENDIF
C*** Assign derived theta-v to profiler level
          thv(K)=
     +    FNEWT(ZFG,thvFG,TABLE2,NZTN_P,1,1,ZW(K,I),TRBL1,NZTN_P)
          IF (thv(K).LE.0.0) THEN
            thv(K)=SPVAL_P+1
          ENDIF

C* Set QC to bogus for derived pressures, and actual heights, so
C* as not to confuse the analysis with the p/z conversion.

          IF (P(K).LT.SPVAL_P) THEN
            O(P_P,ISTA_COUNT)=P(K)
            CALL SETQC(QC(P_P,Ista_count),MSTRBAD_P,
     +                  QC(P_P,ista_count),ISTAT)
            CALL SETQC(QC(P_P,ista_count),BOGSNMC_P,
     +                 QC(P_P,ista_count),ISTAT)
            CALL SETQC(QC(h_P,ista_count),BOGSNMC_P,
     +                 QC(h_P,ista_count),ISTAT)
            CALL SETQC(QC(H_P,Ista_count),MSTRBAD_P,
     +                 QC(H_P,Ista_count),ISTAT)
          ELSE
            write(6,*) 'Something wrong with PROFIN p/z,',P(K),ZW(K,I)
            ista_count=ista_count-1
            goto 80
          ENDIF

          IF (thv(K).LT.SPVAL_P) THEN
            O(theta_P,ISTA_COUNT)=thv(K)
            CALL SETQC(QC(theta_P,Ista_count),MSTRBAD_P,
     +                  QC(theta_P,ista_count),ISTAT)
            CALL SETQC(QC(theta_P,ista_count),BOGSNMC_P,
     +                 QC(theta_P,ista_count),ISTAT)
          ELSE
            write(6,*) 'Something wrong with PROFIN thv/z,',
     +              thv(K),ZW(K,I)
            ista_count=ista_count-1
            goto 80
          ENDIF

          O(H_P,ISTA_COUNT)=ZW(K,I)

          OTYPE(ISTA_COUNT) = 203
C                        ! network data
          IF(I.GT.NUM_PROF) THEN 
            OTYPE(ISTA_COUNT) = 210  
C*** subsrc 210 represents VAD winds
          ELSE
C*** Try to insert RASS data, if any
            DO 200 IRASS=1,NUM_OF_RASS
C  changed 6/25/98 - Stan, John B., Tracy since there are
C    not unique station numbers for RASS but there are unique
C    station NAMES.

              IF(STANaM(I).EQ.RSTANaM(IRASS)) THEN
c              IF(STANUM(I).EQ.RSTANUM(IRASS)
c     +           .AND.STANUM(I).NE.99999) THEN
                Do L=1,NTVLEV(IRASS)
                  IF(TV_HGT(L,IRASS).EQ.ZW(K,I)
     +               .and.TV(L,IRASS).LT.SPVAL_P) THEN

                    O(THETA_P,ISTA_COUNT)=
     +                TV(L,IRASS)*(1000./P(K))**ROVCP_P 
                    QC(THETA_P,ISTA_COUNT)=0

                  ENDIF
                Enddo
              ENDIF
200         CONTINUE
          ENDIF


          OID(ista_count) = STANAM(I)
          OI(ista_count) = XC
          OJ(ista_count) = YC
          OTIME(ISTA_COUNT) = 0
          WRITE(LUNIT,535) (O(ivar,ISTA_COUNT), ivar=1,NVAR_P)
 535     FORMAT(1x,7G15.8)


80      CONTINUE
C* levels loop (k)

      if (i.eq.num_prof) then
      write(6,*)'                                   '
      write(6,*)'***********************************'
      write(6,*)'No. of ob levels after profiler only',
     +    ista_count
      write(6,*)'***********************************'
      end if

100   CONTINUE
C                      !End of station loop (i)

      write(6,*)'                                   '
      write(6,*)'***********************************'
      write(6,*)'No. of ob levels after  profiler/VAD ',
     +    ista_count
      write(6,*)'***********************************'
      
999   CLOSE(UNIT=LUNIT, STATUS='KEEP')
      CALL GETLUN(LUNIT,2)
      RETURN
      END





