      PROGRAM SNDPST
C$$$  MAIN PROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C MAIN PROGRAM: ETA_SNDP
C   PRGMMR: ROGERS           ORG: NP22        DATE: 1999-09-24
C
C ABSTRACT:  THIS ROUTINE POSTS PROFILE DATA AND WRITES
C   OUTPUT IN BUFR FORMAT.  THIS REPLACES CODE THAT USED TO
C   RUN INSIDE OF CHKOUT IN THE ETA MODEL.
C     
C PROGRAM HISTORY LOG:
C   95-07-26  MIKE BALDWIN
C   96-05-07  MIKE BALDWIN - USE SMASK TO SET SOIL VARS TO MISSING
C   96-11-22  MIKE BALDWIN - ADD OPTION OF DOING MONOLITHIC FILE OR
C                            SPLIT OUT FILES OR BOTH
C   97-12-16  MIKE BALDWIN - NEW MULTI-LEVEL PARAMETERS (SUCH
C                            AS SOIL MOISTURE)
C   98-07-23  ERIC ROGERS  - MADE Y2K COMPLIANT
C   98-09-29  MIKE BALDWIN - SET ACC/AVE VARS TO MISSING AT T=0
C   99-04-01  GEOFF MANIKIN - MAJOR CHANGES FEATURING A REMOVAL OF
C                            THE DISTNICTION OF CLASS0 - ALL OUTPUT
C                            IS NOW CLASS1.  ALSO, THE FIELDS OF
C                            VISIBILITY AND CLOUD BASE PRESSURE 
C                            ARE ADDED.
C   99-09-03  JIM TUCCILLO - REDUCED MEMORY REQUIREMENTS BY CHANGING
C                            THE SIZE OF PRODAT AND INTRODUCING LPRO.
C                            ALSO, PRODAT'S STRUCTURE WAS CHANGED TO
C                            PROVIDE STRIDE-1 ACCESS.
C                            NOTE: THIS CODE CAN STILL BE MODIFIED TO
C                            REDUCE MEMORY. THE CHANGES TODAY WILL
C                            NOT AFFECT ITS FUNCTIONALITY BUT WILL
C                            ALLOW IT TO RUN ON A WINTERHAWK NODE
C                            WITHOUT PAGING. THE MEMORY REQUIREMENT
C                            IS NOW AT ABOUT 260 MBs.
C   00-03-10  GEOFF MANIKIN - CODE CHANGED TO BE READY FOR ETA EXTENSION
C                            TO 60 HOURS
C   00--5-15  ERIC ROGERS   - PUT NSOIL AND LM1 IN INCLUDED PARAMETER 
C                             FILE
C     
C USAGE:   
C   INPUT ARGUMENT LIST:
C     NONE    
C     
C   OUTPUT FILES:
C     NONE
C     
C   SUBPROGRAMS CALLED:
C     UTILITIES:
C       CALWXT
C       CALHEL
C       BFRIZE
C     LIBRARY:
C       BUFRLIB     
C     
C   ATTRIBUTES:
C     LANGUAGE: FORTRAN
C     MACHINE : CRAY C-90
C$$$  
C     
C***
C***   7/26/95 M. BALDWIN
C*** 
C***     SOUNDING POST PROCESSOR
C***     PROGRAM TO READ THE PROFILE OUTPUT FILE ON THE CRAY
C***     AND PRODUCE DIAGNOSTIC QUANTITIES AND PACK INTO BUFR
C***
C--------------------------------------------------------------------
C
C    PARMS FOR HOURLY PROFILER OUTPUT
C      LM    - MAX NUMBER OF VERTICAL LEVELS
C      NPNT  - MAX NUMBER OF OUTPUT TIMES    
C      
C          TO HOLD ALL VARIABLES FOR ANY CLASS OF OUTPUT
C          (MAX NO MULTI-LAYER VARIABLES*LM + NO OF SINGLE LAYER VARS)
C      LCL1ML - MAX NUMBER OF MULTI-LAYER VARIABLES FOR CLASS 0 OR 1
C      LCL1SL - MAX NUMBER OF SINGLE LAYER VARIABLES FOR CLASS 0 OR 1
C      LCL1SOIL - MAX NUMBER OF SOIL LAYER VARIABLES FOR CLASS 0 OR 1
C      NSTAT - MAX NUMBER OF STATIONS
C
C        NOTE: THESE NUMBERS WILL BE LARGER THAN THE NUMBERS
C              COMING OUT OF THE MODEL IN THE BINARY FILE SINCE
C              WE ARE COMPUTING SOME ADDITIONAL VARIABLES TO GO INTO
C              BUFR IN THIS PROGRAM.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C	Different philosophy for workstation Eta:
C
C	get actual values for NSTAT, LM, and NPNT.  Reducing these
C	values increases chances that code will not seg fault.  SNDP
C	makes major memory demands on a computer
C
C--------------------------------------------------------------------
        INCLUDE "parmeta"
	INCLUDE "parmsndp"
	INCLUDE "nsta.txt"

     	PARAMETER (SPVAL=-99999.0,SMISS=1.E10
     &, LCL1ML=13,LCL1SL=58,LCL1SOIL=2
     &, LCL1ML1=13,LCL1SL1=50 
     &, NWORD=(LCL1ML+1)*LM+2*LCL1SL+NSOIL*LCL1SOIL
     &, ROG=287.04/9.8)
      PARAMETER (SNOCON=1.4594E5,RAINCON=1.1787E4)
C
      LOGICAL LVLWSE,SEQFLG(8),NEED
      CHARACTER*16 SEQNM1(8), SBSET
      CHARACTER*80 CLIST1(8),FMTO,ASSIGN
      CHARACTER*8 CISTAT
cwas  DIMENSION FPACK(NWORD),PRODAT(NSTAT,NPNT,NWORD)
      DIMENSION FPACK(NWORD)
      REAL PRODAT(NWORD,NSTAT,NPNT)
C
C     THE PURPOSE OF LPRO IS TO HOLD THE VALUES OF RISTAT UNTIL
C     THEY ARE COPIED TO FRODAT. THIS ADDITION WILL ALLOW PRODAT
C     TO BE A REAL(4) ARRAY ( AND SAVE A CONSIDERABLE AMOUNT OF 
C     MEMORY) . PRODAT CAN BE FURTHER REDUCED WITH SOME MORE EFFORT.
C                    JIM TUCCILLO
C
      REAL(8) LPRO(NSTAT,NPNT)
      REAL(8) RISTAT 
cwas  REAL(8) PRODAT(NSTAT,NPNT,NWORD),RISTAT 
      REAL(8) FRODAT(NWORD),WORKK(NWORD)
      DIMENSION P(LM),T(LM),U(LM),V(LM),Q(LM),PINT(LM+1),ZINT(LM+1)
      REAL CWTR(LM),IMXR(LM)
      INTEGER IDATE(3),NP1(8),LLMH(NSTAT),NLVL(2)
      EQUIVALENCE (CISTAT,RISTAT)
C--------------------------------------------------------------------     
C
C     SET OUTPUT UNITS FOR CLASS 1 PROFILE FILE.
C      LCLAS1 - OUTPUT UNIT FOR CLASS 1 BINARY FILE
C      LTBCL1 - INPUT UNIT FOR CLASS 1 BUFR TABLE FILE
C      LUNCL1 - OUTPUT UNIT FOR CLASS 1 BUFR FILE
C
C--------------------------------------------------------------------     
                            I N T E G E R
     & LCLAS1,LTBCL1,LUNCL1,STDOUT
C--------------------------------------------------------------------     
                            L O G I C A L
     & MONOL,BRKOUT
C--------------------------------------------------------------------     
       NAMELIST /MODTOP/ ETOP
       NAMELIST /OPTION/ MONOL,BRKOUT
                            D A T A
     & LCLAS1 / 76 /
     &,LTBCL1 / 32 /
     &,LUNCL1 / 78 /
     &,STDOUT / 6 /
     &,SEQNM1 /'HEADR','PROFILE','SURF','FLUX',
     &         'HYDR','D10M','SLYR','XTRA'/
     &,SEQFLG /.FALSE.,.TRUE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,
     &   .TRUE.,.FALSE./
     &,LVLWSE /.TRUE./
      CALL W3TAGB('ETA_SNDP',1999,0267,0084,'NP22')


      PRODAT=0.
      FMTO='("ln -s ${DIRD}",I5.5,".",I4.4,3I2.2,'//
     &     '"  fort.",I2.2)'
C
C   GET MODEL TOP PRESSURE
C
       PTOP=25.0*100.0
       READ(5,MODTOP,END=12321)
12321  CONTINUE
       PTOP=ETOP*100.0
C
C   READ IN SWITCHES TO CONTROL WHETHER TO DO...
C     MONOL=.TRUE.   DO MONOLITHIC FILE
C     BRKOUT=.TRUE.  DO BREAKOUT FILES
C
       MONOL=.TRUE.
       BRKOUT=.FALSE.
       READ(11,OPTION,END=12322)
12322  CONTINUE
C
        IFCSTL=-99
        JHR=0
C----------------------------------------------------------------------
C---READ STATION DATA--------------------------------------------------
C----------------------------------------------------------------------
      LUNIT=66
      LRECPR=4*(8+9+LCL1ML1*LM1+LCL1SL1)
      OPEN(UNIT=LUNIT,ACCESS='DIRECT',RECL=LRECPR,IOSTAT=IER,
     +		FILE='profilm.c1.t00s')
      NREC=0
 33   CONTINUE
      NREC=NREC+1
      READ(LUNIT,REC=NREC,ERR=999) IHRST,IDATE,IFCST,ISTAT,CISTAT,  
     &   (FPACK(N),N=1,9),(FPACK(N),N=10,FPACK(7))
        IF (IFCST.GT.IFCSTL) THEN
         INUMS=1
         JHR=JHR+1
         IFCSTL=IFCST
        ELSE
         INUMS=INUMS+1
        ENDIF
        IYR=IDATE(3)
        IMON=IDATE(1)
        IDAY=IDATE(2)
!        RLAT=FPACK(1)*DEG
!        RLON=FPACK(2)*DEG
!        ELEV=FPACK(3)
        LLMH(INUMS)=NINT(FPACK(4))
        LMH=NINT(FPACK(4))
        DO 25 L=1,LMH
C   REVERSE ORDER SO THAT P(1) IS THE TOP AND P(LMH) IS THE BOTTOM
         LV=LMH-L+1
         P(LV)=FPACK(L+9)
         T(LV)=FPACK(L+9+LMH)
         U(LV)=FPACK(L+9+LMH*2)
         V(LV)=FPACK(L+9+LMH*3)
         Q(LV)=FPACK(L+9+LMH*4)
C   IF THE CLOUD WATER IS NEGATIVE, IT IS ICE
         CWTR(LV)=FPACK(L+9+LMH*6)
         IF (CWTR(LV).LT.0.) THEN
           IMXR(LV)= -1. * CWTR(LV)
           CWTR(LV)= 0. 
           FPACK(L+9+LMH*6) = 0.
         ELSE
           IMXR(LV) = 0.
         ENDIF
 25     CONTINUE

C  USE SEA MASK TO SET SOIL/SFC VARIABLES TO MISSING VALUES
C    (IF SEA)
C
         SM  =FPACK(13*LMH+54)
         IF (SM.GT.0.5) THEN
C   SMSTAV
          FPACK(13*LMH+15)=SMISS
C   SUBSHX
          FPACK(13*LMH+21)=SMISS
C   SNOPCX
          FPACK(13*LMH+22)=SMISS
C   ACSNOW, SMSTOT, SNO, ACSNOM, SSROFF, BGROFF, SOILTB
          DO LKJ=20,26
            FPACK(13*LMH+LKJ+9)=SMISS
          ENDDO
C   SFCEXC, VEGFRC, CMC, SMC(1:4), STC(1:4)
          DO LKJ=34,44
            FPACK(13*LMH+LKJ+9)=SMISS
          ENDDO
         ENDIF
C
C  GET PPT FOR CALWXT
C
         PPT  =FPACK(13*LMH+16)
C     COMPUTE PINT,ZINT
C
        PINT(1)=PTOP
        DO L=1,LMH
          DP1=P(L)-PINT(L)
          PINT(L+1)=P(L)+DP1
        ENDDO
        ZINT(LMH+1)=FPACK(3)
        DO L=LMH,1,-1
         TV2=T(L)*(1.0+0.608*Q(L))
         ZZ=ROG*TV2*ALOG(PINT(L+1)/PINT(L))
         ZINT(L)=ZINT(L+1)+ZZ
        ENDDO
C
C     CALL PRECIP TYPE SUBROUTINE.
C
      CALL CALWXT(T,Q,P,PINT,LMH,LM,PPT,IWX)
C     
C     DECOMPOSE IWX
C
        CSNO=MOD(IWX,2)
C
        CICE=MOD(IWX,4)/2
C
        CFZR=MOD(IWX,8)/4
C
        CRAI=IWX/8
C
C
C   COMPUTE HELICITY AND STORM MOTION
C
      CALL CALHEL(U,V,P,ZINT,PINT,LMH,LM,HELI,UST,VST)
C
C   COMPUTE VISIBILITY
C   FIRST, EXTRACT THE SNOW RATIO AND SEA LEVEL PRESSURE
      SR=FPACK(13*LMH+58)
      SLP=FPACK(13*LMH+10)

       if (PPT.LT.0.)then
         PPT=0.0
       endif

       SNORATE=(SR/100.)*PPT/3600.
       RAINRATE=(1-(SR/100.))*PPT/3600.
       TERM1=(T(LMH)/SLP)**0.4167
       TERM2=(T(LMH)/(P(LMH)))**0.5833
       TERM3=RAINRATE**0.8333
       QRAIN=RAINCON*TERM1*TERM2*TERM3
       TERM4=(T(LMH)/SLP)**0.47
       TERM5=(T(LMH)/(P(LMH)))**0.53
       TERM6=SNORATE**0.94
       QSNO=SNOCON*TERM4*TERM5*TERM6
       TT=T(LMH)
       QV=Q(LMH)
       QCD=CWTR(LMH)
       QICE=IMXR(LMH)
       PPP=P(LMH)
       CALL CALVIS(QV,QCD,QRAIN,QICE,QSNO,TT,PPP,HOVI)

C   COMPUTE CLOUD BASE PRESSURE
C   FIRST, EXTRACT THE CONVECTIVE CLOUD BASE
       HBOT=FPACK(13*LMH+59)
       CLIMIT =1.0E-06
       NEED = .TRUE.
       CDBP = SMISS
       CBOT = 5000
       DO L=LMH,1,-1
C GSM
C START AT THE FIRST LAYER ABOVE GROUND, AND FIND THE
C   FIRST LAYER WITH A VALUE OF CLOUD WATER GREATER THAN
C   THE SIGNIFICANT LIMIT (VALUE DESIGNATED BY Q. ZHAO).
C   THIS LAYER WILL BE THE CLOUD BOTTOM UNLESS THE BOTTOM
C   OF THE CONVECTIVE CLOUD (HBOT) IS FOUND BELOW IN WHICH
C   CASE HBOT BECOMES THE CLOUD BASE LAYER.

        IF ((CWTR(L)+IMXR(L)).GT.CLIMIT.AND.NEED) THEN
            CBOT=L
            IF (HBOT.GT.CBOT) THEN
              CBOT = HBOT
            ENDIF
            NEED=.FALSE.
          ENDIF
        ENDDO


        IF (CBOT.GT.LMH) THEN
          CDBP=SMISS
        ELSE
          CDBP=P(INT(CBOT))
        ENDIF
         
C
C
C   SET ACC/AVERAGED VARIABLES TO MISSING IF IFCST=0
C
      IF (IFCST.EQ.0) THEN
          DO L=1,LMH
           FPACK(L+9+LMH*7)=SMISS
           FPACK(L+9+LMH*8)=SMISS
          ENDDO
          DO JK=16,29
           FPACK(13*LMH+JK)=SMISS
          ENDDO
          DO JK=32,34
           FPACK(13*LMH+JK)=SMISS
          ENDDO
      ENDIF
C

C      ADD 9 SINGLE LEVEL VARIABLES TO THE OUTPUT
C      TACK THEM ON TO THE END;  WE DON'T NEED CONVECTIVE
C      CLOUD BASE, THOUGH, SO WRITE OVER THAT RECORD
C
          NLEN = FPACK (7) 
          FPACK(NLEN) = CSNO
          FPACK(NLEN+1) = CICE
          FPACK(NLEN+2) = CFZR
          FPACK(NLEN+3) = CRAI
          FPACK(NLEN+4) = UST
          FPACK(NLEN+5) = VST
          FPACK(NLEN+6) = HELI
          FPACK(NLEN+7) = CDBP
          FPACK(NLEN+8) = HOVI
          FPACK(5) = FPACK(5) + 1
          FPACK(6) = FPACK(6) + 8
          FPACK(7) = 9 + FPACK(5)*FPACK(4) + FPACK(6)
C
C           PLACE DATA INTO PRODAT IN PROPER LOCATIONS

          PRODAT (1,INUMS,JHR) = FLOAT(IFCST)
          PRODAT (2,INUMS,JHR) = FLOAT(ISTAT)
C         RISTAT is a REAL(8) variable by virtue of the fact that it 
C         is equivalenced to CISTAT. Everything else stored in PRODAT
C         is REAL(4). We have made PRODAT REAL(4) but need a REAL(8)
C         array for storing RISTAT - that is what LPRO is. Farther
C         down in the code, we will pull values out of LPRO and store
C         in FRODAT ( a REAL(8) array ).
cwas      PRODAT (3,INUMS,JHR) = RISTAT 
          LPRO   (  INUMS,JHR) = RISTAT
          PRODAT (4,INUMS,JHR) = FPACK (1)
          PRODAT (5,INUMS,JHR) = FPACK (2)
          PRODAT (6,INUMS,JHR) = FPACK (3)
          PRODAT (7,INUMS,JHR) = 1 

          DO IJ = 10, 13*LMH+9
            PRODAT (IJ-2,INUMS,JHR) = FPACK (IJ)
          ENDDO
C    TACK ON THE ICE WATER TO THE PROFILE SECTION
C    IT IS CURRENTLY WRITTEN IN REVERSE ORDER.
          DO L=1,LMH
            LV=LMH-L+1
            PRODAT(L+7+LMH*13,INUMS,JHR)=IMXR(LV)
          ENDDO
          DO IJ = 13*LMH+10,NLEN+8
            PRODAT (IJ+LMH-2,INUMS,JHR) = FPACK (IJ)
          ENDDO
        GOTO 33
 999    CONTINUE

	write(6,*) 'past 999'
C
C  WRITE OUT INDIVIDUAL FILES FOR EACH STATION
C
        IF (BRKOUT) THEN
        DO I=1,INUMS
         NLVL(1)=LLMH(I)
         NLVL(2)=NSOIL
C
         DO J=1,JHR
          DO IJ = 1, NWORD
            FRODAT(IJ) = PRODAT (IJ,I,J)
          ENDDO
          FRODAT(3) = LPRO(I,J)
          ISTAT=NINT(FRODAT(2))
  
C
          IF (J.EQ.1) THEN
C
C     INITIALIZE BUFR LISTS SO BFRHDR WILL BE CALLED THE FIRST
C     TIME THROUGH.
C
            WRITE(ASSIGN,FMTO) ISTAT,IYR,IMON,IDAY,IHRST,LCLAS1
            CALL SYSTEM(ASSIGN)
            CLIST1(1)=' '
          ENDIF
C
C           CALL BUFR-IZING ROUTINE
C
           NSEQ = 8
           SBSET = 'ETACLS1'
          CALL BFRIZE(LTBCL1,LCLAS1,SBSET,IYR,IMON,IDAY,IHRST
     1,               SEQNM1,SEQFLG,NSEQ,LVLWSE,FRODAT,NLVL,CLIST1,NP1
     2,               WORKK,IER)
          IF(IER.NE.0)WRITE(6,1080)ISTAT,IER,FRODAT(1)
 1080   FORMAT(' SOME SORT OF ERROR ',2I8,F9.1)
C
C
         ENDDO
C
C   FINISHED, CLOSE UP BUFR FILES
C
        NSEQ = 8
        CALL BFRIZE(0,LCLAS1,SBSET,IYR,IMON,IDAY,IHRST
     1,             SEQNM1,SEQFLG,NSEQ,LVLWSE,FRODAT,NLVL,CLIST1,NP1
     2,             WORKK,IER)
         ENDDO
         ENDIF
         IF (MONOL) THEN
C
C  WRITE OUT ONE FILE FOR ALL STATIONS
C
C     INITIALIZE BUFR LISTS SO BFRHDR WILL BE CALLED THE FIRST
C     TIME THROUGH.
C
            CLIST1(1)=' '
        DO I=1,INUMS
         NLVL(1)=LLMH(I)
         NLVL(2)=NSOIL
C
         DO J=1,JHR
          DO IJ = 1, NWORD
            FRODAT(IJ) = PRODAT (IJ,I,J)
          ENDDO
          FRODAT(3) = LPRO(I,J)
          ISTAT=NINT(FRODAT(2))
C
C           CALL BUFR-IZING ROUTINE
C
           NSEQ = 8
           SBSET = 'ETACLS1'
          CALL BFRIZE(LTBCL1,LUNCL1,SBSET,IYR,IMON,IDAY,IHRST
     1,               SEQNM1,SEQFLG,NSEQ,LVLWSE,FRODAT,NLVL,CLIST1,NP1
     2,               WORKK,IER)
          IF(IER.NE.0)WRITE(6,1080)ISTAT,IER,FRODAT(1)
C
C
         ENDDO
         ENDDO
C
C   FINISHED, CLOSE UP BUFR FILES
C
        NSEQ = 8
        CALL BFRIZE(0,LUNCL1,SBSET,IYR,IMON,IDAY,IHRST
     1,             SEQNM1,SEQFLG,NSEQ,LVLWSE,FRODAT,NLVL,CLIST1,NP1
     2,             WORKK,IER)
         ENDIF
C
        WRITE(STDOUT,*) ' END OF SOUNDING POST '
      CALL W3TAGE('ETA_SNDP')
        STOP
        END
