*********************************************************************   
*********************************************************************   
*****MIXING HEIGHT PROGRAM    
***** 
*********************************************************************   
C  9/24/96 TOM WHITEHURST     
C  8/5/97 CHANGED SEARCH RANGE FOR MIN TEMP TO 02-08 LST    
C         MANY PART TIME STATIONS DON'T START UNTIL 7 AM    
C  TOM WHITEHURST 
C  ADAPTED THE HOLZI FORTRAN PROGRAM FOR USE AT NCDC. 
C  ADDED THE FOLLWING OUTPUT. 
C     
*                     MIXING HEIGHT FORMAT
*     
*RECORD POS       ITEM  
*1-5                              STATION NUMBER
*6-7                              YEAR    
*8-9                              MONTH   
*10                               SEASON  
*11-12                          DAY 
*13                               AM TYPE 
*14-17                          AM MIXING DEPTH 
*18-20                          AM AVERAGE WIND SPEED THRU MIXING 
*21-23                          AM AVERAGE SURFACE WIND SPEED     
*24                               PM TYPE 
*25-28                          PM MIXING DEPTH 
*29-31                          PM AVERAGE WIND SPEED THRU MIXING 
*32-34                          PM AVERAGE SURFACE WIND SPEED     
*     
*SEASON: 1=DECEMBER, JANUARY, FEBRUARY    
*             2=MARCH, APRIL, MAY   
*             3=JUNE, JULY, AUGUST  
*             4=SEPTEMBER, OCTOBER, NOVEMBER    
*     
* TYPE: 1= NO PRECIPITATION   
*          2= PRECIPITATION OCCURRED AT TWO OR MORE OBSERVATION TIMES WITH    
*          THE INTENSITY OF LIGHT, OR AT ONE OBSERVATION WITH THE INTENSITY   
*          OF MODERATE OR HEAVY DURING THE HOURS 1000-2100 LST FOR THE  
*          AFTERNOON OR DURING THE HOURS 2200(PREVIOUS DAY )-0900 LST FOR     
*          THE MORNING. 
*          3=COLD AIR ADVECTION     
*               (A). PM TYPE; THE MAXIMUM TEMPERATURE SELECTED BETWEEN  
*               1200-1600 LST IS LESS THAN THE SURFACE TEMPERATURE OF THE     
*               1200Z RAOB TEMPERATURE.   
*               (B). AM TYPE; THE MINIMUM TEMPERATURE SELECTED BETWEEN  
*               0200-0800 LST PLUS 5 DEGREES CELSIUS IS LESS THAN THE   
*               SURFACE TEMPERATURE OF THE 1200Z RAOB TEMPERATURE.
*          4= MISSING   
*     
*MIXING HGT: MIXING HEIGHT ABOVE SURFACE (WHOLE METERS).    
*               (A). AM: MIXING DEPTH DETERMINED FROM THE INTERSECTION OF THE 
*               1200Z RAOB TEMPERATURE PROFILE AND THE
*POTENTIAL TEMPERATURE BASED ON THE SURFACE RAOB PRESSURE AT 1200 
*MINIMUM    
*HOURLY SURFACE TEMPERATURE   
*DURING THE 0200-0800 LST PLUS 5 DEGREES CELSIUS.     
*                (B). PM: MIXING DEPTH DETERMINED FROM THE INTERS 
*               1200Z RAOB TEMPERATURE PROFILE AND THE
*POTENTIAL TEMPERATURE BASED ON THE SURFACE RAOB PRESSURE AT 1200 
*MAXIMUM HOURLY SURFACE TEMPERATURE 
*DURING THE 1200-1600 LST.    
*               (C). FOR TYPE 3: NO MIXING DEPTH COMPUTED   
*     
*AVERAGE WIND SPEED IN LAYER: WIND SPEED IN METERS PER SECOND AND 
*AVERAGED THROUGH THE MIXING LAYER (SURFACE WIND ONLY IF MIXING D 
*METERS) THE WIND SPEED AT THE TOP OF THE MIXING LAYER IS DETERMI 
*INTERPOLATION IF IT FALLS BETWEEN TWO RAOB LEVELS.   
*               (A).  AM - THE AVERAGE IS DETERMINED FROM THE 1200Z RAOB.     
*               THE AVERAGE OF THE HOURLY SURFACE WIND SPEED DURING 0200-0800 
*               LST WAS USED IN PLACE OF THE SURFACE WIND SPEED IN THE RAOB.  
*               (B).  PM - THE AVERAGE IS DETERMINED FROM THE 0000Z RAOB OF   
*               THE FOLLOWING GREENWICH DAY. THE AVERAGE OF THE HOURLY  
*               SURFACE WIND SPEED DURING 1200-1600 LST WAS USED IN PLACE OF  
*               THE SURFACE WIND SPEED IN THE RAOB.   
*               (C).   FOR TYPE 3 - WIND SPEEDS WERE AVERAGED FROM THE  
*               SURFACE THROUGH 4,000 METERS.   
*     
*AVERAGE WIND SPEED - AVERAGE SURFACE WIND SPEED IN METERS PER SE 
*TENTHS.    
*               (A).  AM - THE AVERAGE OF THE HOURLY SURFACE WIND SPEED 
*               DURING 0200-0800 LST
*               (B).  PM - THE AVERAGE OF THE HOURLY SURFACE WIND SPEED 
*               DURING 1200-1600 LST
*     
*THE MIXING DEPTHS CAN BE COMPUTED WITH 00Z OR 12Z SOUNDING.  THI 
*CONTROLLED BY USER INPUT.    
*     
*IF 12Z OR 00Z SOUNDING IS NOT AVAILABLE THE PROGRAM WILL LOOK FO 
*FOR 12 AND 23Z OR 01Z FOR 00.  IF THE NEXT     
*BEST SOUNDING IS 23Z PROGRAM WILL CONVERT TO 00Z OF NEXT DAY.    
C     
C MODIFIED 4/14/93...BY WILLIAM PUCKETT...
C ADDED VARIABLE TO ALLOW USER TO INPUT FACTOR TO ADD TO MINIMUM TEMPER-
C TURE TO PRODUCE MIN. TEMPERATURE FOR COMPUTING MORNING MIXING HEIGHT  
C ADDED ALGORITHM FOR COMPUTING POTENTIAL TEMPERATURE AT SURFACE..
C ALSO PROGRAM NOW USES FIRST PRESSURE LEVEL (SURFACE) FROM SOUNDING    
C AS STATION PRESSURE RATHER THAN FROM SURFACE DATA (WHICH IS ERRON-    
C EOUSLY GIVING SEA LEVEL PRESSURE IN OUR DATA SETS -- THIS IS A  
C PROBLEM IN HOW OUR DATA SOURCE (WEATHER NETWORK) RETREIVES DATA).     
C MODIFIED 1/27/94...BY WILLIAM PUCKETT...
C CHANGED IWORK2(1,1) TO UAOBS(1,1,1) SO THAT POTENTIAL MINIMUM   
C AND POTENTIAL MAXIMUM TEMPERATURES WERE COMPUTED USING STATION  
C PRESSURE OF CURRENT DAY'S SOUNDING RATHER THAN THE FOLLOWING DAY'S    
C SOUNDING.  THIS ERROR WAS CAUGHT BY STEVEN PREVATT OF TRIGON    
C ENGINEERING CONSULTANTS, INC. (919) 859-1168. 
*AB  5/98.  MARADA Y2K TEAM.  ANN BARRETT 
*AB  INCREASED YEAR FROM I2 TO I4.  INCREASED REC LENGTH FROM 38 TO 40  
*AB  ON TAPE OUTPUT, SYS139.  ALSO MODIFIED PRINTED OUTPUT, MIXPRT.     
*     
C=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
C  PROGRAM HOLZI  
C     
C  PURPOSE  
C     THIS ROUTINE IS USED TO COMPUTE MORNING AND EVENING MIXING  
C     HEIGHTS ACCORDING TO THE TECHNIQUE PRESENTED BY HOLZWORTH (1967), 
C     J. APPL. METEOR., 6, 1039-1044.     
C     
C     THE PROGRAM READS UPPER AIR SOUNDING DATA AND SURFACE DATA  
C     (THE FORMAT IS THE OUTPUT FORMAT PRODUCED  BY THE METEOROLOGICAL  
C     PROCESSOR FOR REGULATORY MODELS (MPRM) STAGE 1 PROCESSOR),  
C     ALIGNS THE UPPER AIR DATA AND SURFACE DATA BY DATE, AND COMPUTES  
C     THE MORNING AND EVENING MIXING HEIGHTS.   
C     
C  I/O UNITS:     5      = USER INPUT (SEE BELOW)     
C                 DEV13  = INPUT UNIT FOR UPPER AIR DATA (INITIALIZED TO
C                          13 IN THE COMMON BLOCKS)   
C                 DEV22  = INPUT UNIT FOR SURFACE DATA (INITIALIZED TO  
C                          22 IN THE COMMON BLOCKS)   
C                 DEV40  = OUTPUT UNIT (INITIALIZED TO 40 IN THE COMMON 
C                          BLOCKS)  
C     
C  FORMATS:       INPUT : REFER TO THE MPRM USER'S GUIDE    
C                 OUTPUT: (A5,3I2,1X,I5,13X,I5) 
C     
C  LIBRARY CALLS: CHROND, GMTLST, GREG, ICHRND, JULIAN,     
C                 FLIWK1, FLIWK2, FLSDG, FLSFC, FLZI, ERROR 
C     
C  SUBROUTINES ATTACHED TO MAIN PGM: UAMOVE     
C     
C  VERSION DATE:  04 JANUARY 1989   
C  REVISION 1  :  28 JUNE 1989    (CORRECT AN ERROR IN DETERMINING THE  
C                                  MINIMUM AND MAXIMUM TEMPERATURE)     
C  REVISION 2  :  14 JULY 1989    (ADDITIONAL COMMENTS)     
C     
C-----------------------------------------------------------------------
****************************************************************  
********* DATA DICTIONARY  *************************************  
****NAME           TYPE     DESCRIPTION   
****AMMIX          REAL     AM MIXING HEIGHT    
****BLNK           CHAR     STRING OF BLANKS    
*     
****CHSTR1         CHAR     INPUT CHARACTER STRING    
****DEV13          INT      FILE POINTER FOR UPPER AIR INPUT
****DEV21          INT      FILE POINTER FOR PRINT FILE     
****DEV22          INT      FILE POINTER FOR SURFACE INPUT  
****DEV40          INT      FILE POINTER FOR EPA DISK OUTPUT
****EOFFLG         LOGI     TEST END OF FILE INPUTS   
****FCTR           REAL     NUMBER TO ADD TO MAX AND MIN TEMPERATURE    
****FIL13          CHAR     UPPER AIR INPUT FILE NAME 
****FIL22          CHAR     SURFACE INPUT FILE NAME   
****FIL40          CHAR     EPA OUTPUT FILE NAME
****I              INT      COUNTER 
****IAMMIX         INT      AM MIXING HEIGHT    
****IATYP          INT      AM TYPE 
****IAWSD          INT      AM AVERAGE UPPER AIR WINDS
****IAWSS          INT      AM AVERAGE SURFACE WINDS  
****ICAHP          INT      AM INDICATOR FOR HEAVY PRECIP   
****ICALP          INT      AM INDICATOR FOR LIGHT PRECIP   
****ICNHP          INT      OVERNITE INDICATOR FOR HEAVY PRECIP   
****ICNLP          INT      OVERNITE INDICATOR FOR LIGHT PRECIP   
****ICPHP          INT      PM INDICATOR FOR HEAVY PRECIP   
****ICPLP          INT      PM INDICATOR FOR LIGHT PRECIP   
****IFLAG          INT      INDICATES BAD MIXING HEIGHT     
****ILEV           INT      COUNTER 
****IOST13         INT      I/O STATUS FILE 13  
****IOST22         INT      I/O STATUS FILE 22  
****IPMMIX         INT      PM MIXING HEIGHT    
****IPTYP          INT      PM TYPE 
****IPWSD          INT      PM AVERAGE OF UPPER AIR WINDS   
****IPWSS          INT      PM AVERAGE SURFACE WINDS  
****ISEASN         INT      SEASON INDICATOR    
****ISWAN          INT      AM NUMBER OF WINDS  
****ISWPN          INT      PM NUMBER OF WINDS  
****IVBL           INT      COUNTER 
****IWORK1(1000)   INT      WORK ARRAY FOR CONTROL AND SURFACE DATA     
****IWORK2(200,6)  INT      WORK ARRAY FOR UPPER AIR DATA   
****JULIAN         INT      JULIAN DAY    
****KDAY           INT      COUNTER 
****LEV            INT      COUNTER 
****LSTGMT         INT      LST TO GMT CONVERSION FACTOR    
****LYRWND         REAL     AVERAGE WIND SPEED IN LAYER     
****MESS           CHAR     ERROR STRING  
****MHR            INT      UPPER AIR HOUR TO COMPUTE MIXING HEIGHT     
****MRDAY1         INT      BEGINNING OF PERIOD TO COMPUTE MIXING HEIGHTS JULIAN    
****MRDAY2         INT      ENDING OF PERIOD TO COMPUTE MIXING HEIGHTS    JULIAN    
****MRGDY1         INT      BEGINNING OF PERIOD TO COMPUTE MIXING HEIGHTS     
****MRGDY2         INT      ENDING OF PERIOD TO COMPUTE MIXING HEIGHTS  
****MRGMO1         INT      BEGINNING OF PERIOD TO COMPUTE MIXING HEIGHTS     
****MRGMO2         INT      ENDING OF PERIOD TO COMPUTE MIXING HEIGHTS  
****MRYR1          INT      BEGINNING OF PERIOD TO COMPUTE MIXING HEIGHTS     
****MRYR2          INT      ENDING OF PERIOD TO COMPUTE MIXING HEIGHTS  
****NV             INT      COUNTER 
****P1             REAL     AM SURFACE PRESSURE 
****P2             REAL     PM SURFACE PRESSURE 
****PMMIX          REAL     PM MIXING HEIGHT    
****PTEM(200)      REAL     ARRAY OF POTENTIAL TEMPERTURES  
****RATIO          REAL     USED FOR LINEAR INTERPOLATION   
****RUNTYP         INT      OUTPUT CONTROL FLAG 
****SETSEA         INT      FNCTION TO SET THE SEASON TYPE  
****SFCTH          REAL     SURFACE POTENTIAL TEMPERATURE   
****SFDAT          LOGI     TEST FOR SURFACE DATE     
****SFDAYC         INT      CURRENT SURFACE OBSERVATION IN JULIAN 
****SFEOF          LOGI     TEST FOR SURFACE END OF FILE    
****SFGDY          INT      DATE OF CURRENT OBSERVATION     
****SFGHR          INT      DATE OF CURRENT OBSERVATION     
****SFGMO          INT      DATE OF CURRENT OBSERVATION     
****SFGYR          INT      DATE OF CURRENT OBSERVATION     
****SFRED          INT      RECORD COUNTER FOR SURFACE FILE 
****STNID          CHAR     STATION NUMBER
****STNNS          CHAR     SURFACE STATION NAME
****STNNU          CHAR     UPPER AIR STATION NAME    
****STNS           CHAR     SURFACE STATION NUMBER    
****STNU           CHAR     UPPER AIR STATION NUMBER  
****SUM4WD         REAL     FUNCTION RETURNS AVERAGE WIND IN LOWER 4000 FEET  
****SWAS           REAL     AM SUM OF SURFACE WINDS   
****SWPS           REAL     PM SUM OF SURFACE WINDS   
****TKELV          REAL     KELVIN TEMPERATURE  
****TMAX           REAL     MAXIMUM TEMPERATURE FROM OBSERVATIONS 
****TMIN           REAL     MINIMUM TEMPERATURE FROM OBSERVATIONS 
****TMIN5          REAL     MINIMUM TEMPERATURE ADJUSTED BY 5     
****TTEMPN         REAL     MINIMUM KELVIN SURFACE TEMPERATURE    
****TTEMPX         REAL     MAXIMUM KELVIN SURFACE TEMPERATURE    
****UADAT          LOGI     TEST FOR UPPER AIR DATA FOR DATE
****UADAYC         INT      CURRENT UPPER AIR DATE IN JULIAN
****UAEOF          LOGI     TEST FOR UPPER AIR FILE END OF FILE   
****UAGDY          INT      CURRENT UPPER AIR DATE    
****UAGHR          INT      CURRENT UPPER AIR DAT     
****UAGMO          INT      CURRENT UPPER AIR DAT     
****UAGYR          INT      CURRENT UPPER AIR DAT     
****UALEV(2)       INT      NUMBER OF UPPER AIR OBS   
****UAOBS(2,200,6) INT      UPPER AIR WORK ARRAY
****UARED          INT      UPPER AIR FILE READ COUNT 
****UAWND          LOGI     TEST FOR UPPER AIR WIND   
****WNDP(200)      REAL     PM WIND PROFILE     
****WNDS(200)      REAL     AM WIND PROFILE     
****Z(200)         REAL     AM HEIGHT PROFILE   
****ZP(200)        REAL     PM HEIGHT PROFILE   
****************************************************************  
      PROGRAM HOLZI     
C     
C     
C     *** NOTE, MXL MUST BE SET TO BE EQUAL TO UAML   
C         AS THIS PROVIDES IN THE WORK ARRAYS PTEMP AND Z   
C         AS MANY LEVELS TO WORK WITH AS THERE ARE IN THE   
C         UA SOUNDINGS. 
C     
      CHARACTER*32 FIL13, FIL22, FIL40    
      CHARACTER*80 CHSTR1     
      CHARACTER*5  STNID,STNU,STNS  
      CHARACTER BLNK40*40,MESS*40,STNNU*30,STNNS*30   
      INTEGER      JULIAN, IOST13,IOST22, NV, LEV,    
     &             IVBL, KDAY, ILEV,MHR,RUNTYP  
      INTEGER      UARED, SFRED, LSTGMT, IAMMIX, IPMMIX, IFLAG    
      INTEGER      IATYP,IPTYP,SETSEA,IWORK1(1000),IWORK2(200,6)  
      INTEGER      I,DEV21,DEV13,DEV40,DEV22,UAOBS(2,200,6) 
      INTEGER      MRDAY1,MRDAY2,MRYR1,MRYR2,MRGDY1,MRGDY2, 
     1             MRGMO1,MRGMO2    
      INTEGER      SFGDY,SFGMO,SFGYR,SFGHR,SFDAYC     
      INTEGER      UAGYR,UAGHR,UADAYC,UAGDY,UAGMO,UALEV(2)  
      REAL         WNDS(200),LYRWND,SUM4WD,ZP(200),WNDP(200)
      REAL         TMIN,TMAX,TMIN5,P1,P2,RATIO,AMMIX,PMMIX, 
     &             PTEM(200),Z(200),TKELV,TTEMPN,TTEMPX,FCTR,SFCTH
      LOGICAL      SFEOF,UAEOF, UADAT,UAWND,SFDAT,EOFFLG    
      PARAMETER (DEV13=13,DEV22=22,DEV21=21,DEV40=40) 
C     
        DATA BLNK40/'                                        '/   
C     
C     
C---  INITIALIZE VARIABLES    
      DATA SFEOF/.FALSE./,UAEOF/.FALSE./,UADAT/.FALSE./,    
     &SFDAT/.FALSE./,MHR/12/,RUNTYP/2/,EOFFLG/.FALSE./
      DATA TMIN/9999.0/, TMIN5/9999.0/, TMAX/-9999.0/ 
      DATA IAMMIX/-9999/, IPMMIX/-9999/, PTEM/200*-9999.0/  
      DATA UARED /0/, SFRED/0/,SWAS/0.0/,SWPS/0.0/    
C     
C-----------------------------------------------------------------------
C- THE FOLLOWING IS A MAP OF THE WORK ARRAYS AND THE MET. DATA    
C   ARE IN THE ARRAY ELEMENTS 
C    IWORK1(10) = UPPER AIR CHRONOLOGICAL DAY (COMPUTATIONAL VBL.)
C    IWORK1(11) = NUMBER OF LEVELS IN A SOUNDING
C    IWORK1(100) = SURFACE OBS. CHRONOLOGICAL DAY (COMPUTATIONAL) 
C    IWORK1(101:122) = SURFACE OBS. DATA, ONE HOUR, TEMPORARY     
C    IWORK1(900:901) = CURRENT COMPUTATIONAL YEAR : JULIAN DAY    
C    IWORK1(910:911) = BEGINNING : ENDING CHRONOLOGICAL DAYS FOR ZI     
C    IWORK1(915) = CURRENT MONTH, RETURNED FROM ROUTINE 'GREG'    
C    IWORK1(916) = CURRENT DAY, RETURNED FROM ROUTINE 'GREG'
C     
C    IWORK2(1:UAML,1:6) = UPPER AIR DATA, ONE SOUNDING, TEMPORARY 
C     
C    IWORK1(20)   = NUMBER OF UPPER AIR OBS. FOR THE DAY    
C    IWORK1(21)   = NUMBER OF SFC. OBS. FOR THE DAY   
C     
C *** OTHER VARIABLE DEFINITIONS    
C     
C    IOST13,-22,    = IOSTAT VALUES RETURNED FROM READING UNITS 13, 22  
C    KDAY           = LOOP INDEX BEGINNING WITH CHRONOLOGICAL DAY 
C                      IWORK1(910) AND ENDING WITH DAY IWORK1(911);     
C    UARED,SFRED    = COUNTERS TRACKING THE NUMBER OF READS OF THE FILES
C     
C    UAEOF, SFOEF   = END-OF-FILE FLAGS   
C    UADAT          = INDICATES IF A 12Z SOUNDING IS AVAILABLE    
C     
C=======================================================================
C- OPEN FILES; FLUSH ALL ARRAYS; INITIALIZE VARIABLES; ENTER REQUIRED   
C   PROCESSING INFORMATION    
C     
C     
      DO 30 K = 1,2     
       DO 35 I = 1,200  
        DO 42 J =1,6    
         UAOBS(K,I,J) = -9999 
   42   CONTINUE  
   35  CONTINUE   
   30 CONTINUE    
C     
      DO 132 I = 1,200  
       DO 146 J = 1,6   
        IWORK2(I,J) = -9999   
  146  CONTINUE   
  132 CONTINUE    
C     
C-----------------------------------------------------------------------
      DO 152 I = 1,1000 
        IWORK1(I) = -9999     
  152 CONTINUE    
C     
C     
C- ENTER THE 5 CHARACTER STATION ID (USED ONLY WHEN WRITING THE OUPUT   
C   RECORD IN THIS PROGRAM, BUT IT IS CHECKED BY RAMMET), THE YEAR
C   (2-DIGIT), MONTH, AND DAY TO START AND END THE COMPUTATIONS,  
C   AND THE NUMBER OF HOURS FROM GREENWICH (+ FOR STATIONS WEST OF
C   GREENWICH).  THIS LAST ENTRY IS REQUIRED TO CONVERT THE DATA TO     
C   GMT. THEREFORE, IF THE UPPER AIR DATA ARE GIVEN IN LOCAL STANDARD   
C   TIME, ENTER THE APPROPRIATE VALUE; HOWEVER IF THE DATA ARE ALREADY  
C   GIVEN IN GMT, ENTER A ZERO (0) FOR THIS VALUE.    
C     
       PRINT*,' ENTER TYPES OF OUTPUTS 1=EPA 2=NCDC ' 
       READ(5,*)RUNTYP  
C     
       PRINT *,' ENTER YEAR, MONTH, DAY (FREE FORMAT, INTEGER) TO START'
       READ(5,*) MRYR1,MRGMO1,MRGDY1
C     
       PRINT *,' ENTER YEAR, MONTH, DAY (FREE FORMAT, INTEGER) TO STOP' 
       READ(5,*) MRYR2,MRGMO2,MRGDY2
C     
       PRINT *,' ENTER NUMBER OF HOURS (INTEGER) FROM GREENWICH -'
       PRINT *,'  FOR WEST LONGITUDE ENTER A POSITIVE NUMBER'     
       READ(5,*) LSTGMT 
       PRINT *,'LST TO GMT = ',LSTGMT     
       IF(RUNTYP.EQ.1)THEN    
       PRINT *,' ENTER THE STATION ID (5 CHARACTERS, USE QUOTES)' 
       READ(5,*) STNID  
C     
C      READ FACTOR FOR ADDING TO 02 - 08 LST TEMPERATURE    
C     
       PRINT *,' ENTER FACTOR (DEGREES C) TO ADD TO 02-08 TEMPERATURE ' 
       READ(5,*) FCTR   
C     
       WRITE (*, *) ' ENTER FILE CONTAINING UPPER AIR DATA '
       READ(5, '(A32)') FIL13 
       WRITE (*, *) ' ENTER FILE CONTAINING SURFACE DATA '  
       READ(5, '(A32)') FIL22 
       WRITE (*, *) ' ENTER EPA    FILE CONTAINING MIXING HEIGHTS '     
       READ(5, '(A32)') FIL40 
       OPEN (DEV40,FILE=FIL40,STATUS='UNKNOWN') 
       OPEN (DEV13,FILE=FIL13,STATUS='OLD')     
       OPEN (DEV22,FILE=FIL22,STATUS='OLD')     
       MHR=12     
       ELSEIF(RUNTYP.EQ.2)THEN

      PRINT*,'ENTER THE U/A STATION NUMBER AND NAME'  
      READ(5,'(A5,1X,A30)')STNU,STNNU     
      PRINT*,'ENTER THE SFC STATION NUMBER AND NAME'  
      READ(5,'(A5,1X,A30)')STNS,STNNS     
       PRINT *,' CHOOSE GMT HOUR TO COMPUTE MIXING DEPTH 0 OR 12' 
       READ(5,*) MHR    
       FCTR=5     
       OPEN (DEV21,FILE='MIXPRT',STATUS='UNKNOWN')
       OPEN (23,FILE='MIXOUT',STATUS='UNKNOWN',MRECL=10)    
*AB       OPEN (24,TYPE='ANSI',MRECL=38,BLOCK=380,RFORM='FB')     
      OPEN (24,TYPE='ANSI',MRECL=40,BLOCK=400,RFORM='FB')   
       OPEN (DEV13,FILE='UADATA',STATUS='OLD')  
       OPEN (DEV22,FILE='SFDATA',STATUS='OLD')  
       ELSE 
       STOP 'OUTPUT CHOOSE ERROR'   
       ENDIF
       
C     
C- COMPUTE THE BEGINNING AND ENDING JULIAN DAYS AND CHRONOLOGICAL 
C-   DAYS (IWORK1(910:911))   
C     
C     
       MRDAY1 = JULIAN(MRYR1,MRGMO1,MRGDY1)     
       CALL CHROND( MRYR1,MRDAY1,IWORK1(910) )  
       IF (IWORK1(910) .LT. 0) THEN 
         MESS = BLNK40  
         WRITE(MESS,142)
         PRINT *,' CHROND ERROR COMPUTING START DATE: ',MESS
         STOP     
       ENDIF
C     
       MRDAY2 = JULIAN(MRYR2,MRGMO2,MRGDY2)     
       CALL CHROND( MRYR2,MRDAY2,IWORK1(911) )  
       IF (IWORK1(911) .LT. 0) THEN 
         MESS = BLNK40  
         WRITE(MESS,142)
         PRINT *,' CHROND ERROR COMPUTING STOP DATE: ',MESS 
         STOP     
       ENDIF
C     
C=======================================================================
C- READ THE HEADERS FROM EACH FILE; A HEADER RECORD IS DEFINED BY AN    
C   ASTERISK (*) IN COLUMN 1  
C     
  120 READ (DEV13,40,IOSTAT=IOST13,ERR=13008,END=13030) CHSTR1    
   40 FORMAT(A80) 
      IF (CHSTR1(1:1) .EQ. '*') THEN
        GO TO 120 
      ELSE  
        BACKSPACE 13    
      ENDIF 
C     
  121 READ (DEV22,40,IOSTAT=IOST22,ERR=22008,END=22030) CHSTR1    
      IF (CHSTR1(1:1) .EQ. '*') THEN
        GO TO 121 
      ELSE  
        BACKSPACE 22    
      ENDIF 
C     
C====================================================================== 
C- READ EACH FILE AND RETRIEVE THE FIRST OBSERVATION IN EACH FILE 
C     
C  PROCESS THE FIRST OBSERVATION OF UPPER AIR DATA    
C     
       READ(DEV13,1350,IOSTAT=IOST13,ERR=13010,END=13030)   
     &      UAGYR, UAGMO, UAGDY, UAGHR, (IWORK1(NV),NV=11,13)     
       DO 131 LEV = 1,IWORK1(11)    
        READ(DEV13,1355,IOSTAT=IOST13,ERR=13020,END=13030)  
     &    (IWORK2(LEV,IVBL),IVBL=1,6)     
  131  CONTINUE   
       UARED = UARED + 1
C     
C- COMPUTE THE JULIAN DAY OF THE RECORD JUST READ AND CONVERT IT TO THE 
C   CHRONOLOGICAL DAY (1 = 1 JANUARY 1900)
       UADAYC = JULIAN(UAGYR, UAGMO, UAGDY)     
       CALL CHROND( UAGYR,UADAYC,IWORK1(10) )   
       IF (IWORK1(10).LT.0) THEN    
         MESS = BLNK40  
         WRITE(MESS,140) 13   
         PRINT *,' CHROND ERROR AT FIRST UA OBS: ',MESS     
         STOP     
       ENDIF
C     
       GO TO 22   
C     
C---------------------------------------------------------------------- 
C- BRANCHING FOR ERROR AND END-OF-FILE CONDITIONS ON 1ST UPPER AIR READ 
C     
13008  PRINT *,'  ERROR READING THE UA FILE HEADERS, I/O STATUS=',IOST13
       STOP 
13010  MESS = BLNK40    
       WRITE(MESS,255)  IOST13
       GO TO 13025
13020  MESS = BLNK40    
       WRITE(MESS,256)  IOST13
13025  PRINT *,' ERROR AT FIRST UA OBS: ',MESS  
       STOP 
C     
13030  UAEOF = .TRUE.   
C     
C=======================================================================
C- PROCESS THE FIRST OBSERVATION OF SURFACE DATA
C     
22     READ(DEV22,2265,IOSTAT=IOST22,ERR=22010,END=22030)   
     &      SFGYR,SFGMO, SFGDY, SFGHR, (IWORK1(100+IVBL),IVBL=1,22)     
       SFRED = SFRED + 1
C     
C- COMPUTE THE JULIAN DAY OF THE RECORD JUST READ AND CONVERT IT TO THE 
C   CHRONOLOGICAL DAY (1 = 1 JANUARY 1900)
       SFDAYC = JULIAN(SFGYR,SFGMO,SFGDY) 
       IF (SFGHR .EQ. 0) THEN 
         CALL GREG(SFGYR,SFDAYC,SFGMO,SFGDY)    
       ENDIF
       CALL CHROND( SFGYR,SFDAYC,IWORK1(100) )  
       IF (IWORK1(100) .LT. 0) THEN 
         MESS = BLNK40  
         WRITE(MESS,140) 22   
         PRINT *,' CHROND ERROR AT FIRST SF OBS: ',MESS     
         STOP     
       ENDIF
C     
       GO TO 9001 
C     
C-----------------------------------------------------------------------
C- BRANCHING FOR ERROR AND END-OF-FILE CONDITIONS FOR 1ST SFC. OBS. READ
22008  PRINT *,' ERROR READING THE SF FILE HEADERS, I/O STATUS=',IOST22 
       STOP 
22010  MESS = BLNK40    
       WRITE(MESS,267)  IOST22
       PRINT *,' ERROR AT FIRST SF OBS: ',MESS  
       STOP 
C     
22030  SFEOF = .TRUE.   
C     
C=======================================================================
C- BEGIN ITERATION ON DAYS    
9001  DO 9000 KDAY = IWORK1(910),IWORK1(911)    
C     
C-----------------------------------------------------------------------
C-    PROCESS THE UPPER AIR DATA    
C     
C-    THE UPPER AIR DAY IS BEFORE THE COMPUTATIONAL DAY, SO JUST READ   
C     THE SOUNDINGS UNTIL WE ENCOUNTER THE COMPUTATIONAL DAY (KDAY) OR  
C     LATER (THE 13001 LOOP)  
C     
13000 IF (.NOT. UAEOF) THEN   
13001  IF (IWORK1(10) .LT. KDAY) THEN     
C     
        IWORK1(11) = 0  
        READ(DEV13,1350,IOSTAT=IOST13,ERR=13011,END=13031)  
     &      UAGYR, UAGMO, UAGDY, UAGHR, (IWORK1(NV),NV=11,13)     
        DO 134 LEV = 1,IWORK1(11)   
         READ(DEV13,1355,IOSTAT=IOST13,ERR=13021,END=13032) 
     &    (IWORK2(LEV,IVBL),IVBL=1,6)     
  134   CONTINUE  
        UARED = UARED + 1     
C     
        UADAYC = JULIAN(UAGYR, UAGMO, UAGDY)    
        CALL CHROND( UAGYR,UADAYC,IWORK1(10) )  
        IF (IWORK1(10) .LT. 0) THEN 
         MESS = BLNK40  
         WRITE(MESS,140) 13   
         PRINT *,' CHROND ERROR AT UA OBS. #:',UARED,' : ',MESS   
         STOP     
        ENDIF     
C     
        GO TO 13001     
C     
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C- THE UPPER AIR DAY IS EQUAL TO THE COMPUTATIONAL DAY; IF IT'S 1200 GMT
C   THEN IT'S THE ONE WE WANT.  USE THE SUBROUTINE UAMOVE   
C   DATA FROM TEMPORARY WORK ARRAYS TO THE UPPER AIR ARRAYS IF THERE ARE
C   ANY LEVELS OF DATA (STORED IN IWORK1(11))   
C     
       ELSE IF (IWORK1(10) .EQ. KDAY) THEN
C     
      IGMT= UAGHR+LSTGMT
      IF(IGMT.EQ.24)THEN
         IGMT=0   
      ELSEIF(IGMT.LT.0)THEN   
         IGMT=24+IGMT   
      ELSEIF(IGMT.GT.24)THEN  
         IGMT=IGMT-24   
      ENDIF 
       IF(MHR.EQ.12)THEN
        IF ((IGMT).GE.10.AND.(IGMT).LE.12)THEN  
          IF ( IWORK1(11) .GT. 0 ) THEN   
            UADAT = .TRUE.    
            IWORK1(20) =  1   
            CALL UAMOVE(UALEV,IWORK1,IWORK2,UAOBS)    
          ELSE    
            UADAT = .FALSE.   
          ENDIF   
        ELSEIF ((IGMT).LE.15.AND..NOT.UADAT) THEN     
          IF ( IWORK1(11) .GT. 0 ) THEN   
            UADAT = .TRUE.    
            IWORK1(20) =  1   
            CALL UAMOVE(UALEV,IWORK1,IWORK2,UAOBS)    
          ELSE    
            UADAT = .FALSE.   
          ENDIF   
      ELSEIF ((IGMT).EQ.23.OR.(IGMT).EQ.0)THEN  
          IF ( IWORK1(11) .GT. 0 ) THEN   
            UAWND = .TRUE.    
            IWORK1(20) =  2   
            CALL UAMOVE(UALEV,IWORK1,IWORK2,UAOBS)    
          ELSE    
            UAWND = .FALSE.   
          ENDIF   
        ELSEIF ((IGMT).LE.3.AND..NOT.UAWND) THEN
          IF ( IWORK1(11) .GT. 0 ) THEN   
            UAWND = .TRUE.    
            IWORK1(20) =  2   
            CALL UAMOVE(UALEV,IWORK1,IWORK2,UAOBS)    
          ELSE    
            UAWND = .FALSE.   
          ENDIF   
        ENDIF     
       ELSEIF(MHR.EQ.0)THEN   
        IF ((IGMT).EQ.23.OR.(IGMT).EQ.0)THEN    
          IF ( IWORK1(11) .GT. 0 ) THEN   
            UADAT = .TRUE.    
            IWORK1(20) =  1   
            CALL UAMOVE(UALEV,IWORK1,IWORK2,UAOBS)    
          ELSE    
            UADAT = .FALSE.   
          ENDIF   
        ELSEIF ((IGMT).LE.3.AND..NOT.UADAT) THEN
          IF ( IWORK1(11) .GT. 0 ) THEN   
            UADAT = .TRUE.    
            IWORK1(20) =  1   
            CALL UAMOVE(UALEV,IWORK1,IWORK2,UAOBS)    
          ELSE    
            UADAT = .FALSE.   
          ENDIF   
      ELSEIF ((IGMT).GE.10.AND.(IGMT).LE.12)THEN
          IF ( IWORK1(11) .GT. 0 ) THEN   
            UAWND = .TRUE.    
            IWORK1(20) =  2   
            CALL UAMOVE(UALEV,IWORK1,IWORK2,UAOBS)    
          ELSE    
            UAWND = .FALSE.   
          ENDIF   
        ELSEIF ((IGMT).LE.15.AND..NOT.UAWND) THEN     
          IF ( IWORK1(11) .GT. 0 ) THEN   
            UAWND = .TRUE.    
            IWORK1(20) =  2   
            CALL UAMOVE(UALEV,IWORK1,IWORK2,UAOBS)    
          ELSE    
            UAWND = .FALSE.   
          ENDIF   
        ENDIF     
       ENDIF
C     
C- READ THE NEXT SOUNDING     
        IWORK1(11) = 0  
        READ(DEV13,1350,IOSTAT=IOST13,ERR=13011,END=13031)  
     &       UAGYR, UAGMO, UAGDY, UAGHR, (IWORK1(NV),NV=11,13)    
        DO 137 LEV = 1,IWORK1(11)   
         READ(DEV13,1355,IOSTAT=IOST13,ERR=13021,END=13032) 
     &    (IWORK2(LEV,IVBL),IVBL=1,6)     
  137   CONTINUE  
        UARED = UARED + 1     
C     
        UADAYC = JULIAN(UAGYR, UAGMO, UAGDY)    
        CALL CHROND( UAGYR,UADAYC,IWORK1(10) )  
        IF (IWORK1(10) .LT. 0) THEN 
          MESS = BLNK40 
          WRITE(MESS,140) 13  
          PRINT *,' CHROND ERROR AT UA OBS.',UARED,' : ',MESS     
          STOP    
        ENDIF     
C     
C- THIS WILL KEEP READING UNTIL THE SOUNDING IS FOR A SUBSEQUENT DAY    
        GO TO 13001     
C     
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C-  THE UPPER AIR DAY IS AFTER THE COMPUTATIONAL DAY  
       ELSE IF (IWORK1(10) .GT. KDAY) THEN
          GO TO 22000   
       ENDIF
      ENDIF 
      GO TO 22000 
C-----------------------------------------------------------------------
C- PROCESSING CONTINUES HERE IF THERE WAS AN ERROR READING THE UPPER    
C-     AIR FILE   
C      (13011 = HDR ERROR, 13021 = LVL ERROR, 13032 = EOF LVL)    
C     
13011  MESS = BLNK40    
       WRITE(MESS,255)  IOST13
       GO TO 13026
13021  MESS = BLNK40    
       WRITE(MESS,256)  IOST13
       GO TO 13026
13032  MESS = BLNK40    
       WRITE(MESS,257) IOST13 
13026  PRINT *,' ERROR READING UA DATA AT OBS.',UARED,' : ',MESS  
       STOP 
C     
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C- PROCESSING CONTINUES HERE IF AN END OF FILE WAS ENCOUNTERED    
C     
13031  UAEOF  = .TRUE.  
C     
C=======================================================================
C- PROCESS THE SURFACE DATA   
C     
C     THE SURFACE OBS. DAY IS BEFORE THE COMPUTATIONAL DAY, SO JUST READ
C     THE  HOURLY OBSERVATIONS UNTIL THE COMPUTATIONAL DAY (KDAY) OR    
C     LATER IS ENCOUNTERD (THE 22001 LOOP)
C     
22000 IF (.NOT. SFEOF) THEN   
22001  IF (IWORK1(100) .LT. KDAY) THEN    
****************************************************************  
********** SOFTWARE ADDED BY TOM WHITEHURST*********************  
       IF (IWORK1(100).EQ.KDAY-1)THEN     
           CALL SETTYP(IWORK1(115),SFGHR,ICALP,ICAHP,ICPLP,ICPHP, 
     &ICNLP,ICNHP)
       ENDIF
***   
        READ(DEV22,2265,IOSTAT=IOST22,ERR=22011,END=22031)  
     &       SFGYR,SFGMO, SFGDY, SFGHR, (IWORK1(100+IVBL),IVBL=1,22)    

C     
        SFRED = SFRED + 1     
        SFDAYC = JULIAN(SFGYR,SFGMO,SFGDY)
        CALL CHROND( SFGYR,SFDAYC,IWORK1(100) ) 
C     
        IF (IWORK1(100) .LT. 0) THEN
          MESS = BLNK40 
          WRITE(MESS,140) 22  
          PRINT *,' CHROND ERROR AT SF OBS.',SFRED,' : ',MESS     
          STOP    
        ENDIF     
C     
        GO TO 22001     
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       ELSE IF (IWORK1(100) .EQ. KDAY) THEN     
C- ONLY THE DATE AND HOUR AND TEMPERATURE AND PRESSURE INFORMATION ARE  
C   NEEDED.  TEMPERATURE MUST BE IN DEGREES CENTIGRADE, PRESSURE IN MB. 
C   MISSING VALUES ARE DENOTED BY -9999   
        SFDAT              = .TRUE. 
        IWORK1(21)         = IWORK1(21) + 1     
****************************************************************  
********** SOFTWARE ADDED BY TOM WHITEHURST*********************  
           CALL SETTYP(IWORK1(115),SFGHR,ICALP,ICAHP,ICPLP,ICPHP, 
     &ICNLP,ICNHP)
***   
C     
C-  DETERMINE THE MINIMUM TEMPERATURE (TTEMPN) FOR THE DAY BETWEEN
C    THE HOURS OF 0200 LST AND 0800 LST, INCLUSIVE    
C     
        IF ((SFGHR .GE. 2) .AND. (SFGHR .LE. 8)) THEN 
          IF (IWORK1(117) .GT. -9000.0) THEN    
            TTEMPN = IWORK1(117)/10.0 +273.16   
            IF (TTEMPN .LT. TMIN) THEN    
              TMIN = TTEMPN   
              P1   = IWORK1(103)/10.0     
            ENDIF 
          ENDIF   
*********************************************************************   
**************BEGIN**************************************************   
*****SOFTWARE DEVELOPED BY TOM WHITEHURST 
***** SUM THE SURFACE WINDS   
*********************************************************************   
       IF(IWORK1(122).GT.-9000)THEN 
           SWAS =  SWAS + IWORK1(122)     
          ISWAN = ISWAN + 1   
       ENDIF
*********************************************************************   
**************END****************************************************   
*****SOFTWARE DEVELOPED BY TOM WHITEHURST 
***** 
*********************************************************************   
C     
C-  DETERMINE THE MAXIMUM TEMPERATURE (TTEMPX) FOR THE DAY BETWEEN
C    THE HOURS OF 1200 LST AND 1600 LST, INCLUSIVE    
C     
        ELSE IF ((SFGHR .GE. 12) .AND. (SFGHR.LE.16)) THEN  
          IF (IWORK1(117) .GT. -9000.0) THEN    
            TTEMPX = IWORK1(117)/10.0 + 273.16  
            IF (TTEMPX .GT. TMAX) THEN    
              TMAX = TTEMPX   
              P2   = IWORK1(103)/10.0     
            ENDIF 
          ENDIF   
*********************************************************************   
***************BEGIN*************************************************   
*****SOFTWARE DEVELOPED BY TOM WHITEHURST 
***** SUM THE SURFACE WINDS   
*********************************************************************   
       IF(IWORK1(122).GT.-9000)THEN 
           SWPS =  SWPS + IWORK1(122)     
          ISWPN = ISWPN + 1   
       ENDIF
*********************************************************************   
****************END**************************************************   
*****SOFTWARE DEVELOPED BY TOM WHITEHURST 
***** 
*********************************************************************   
        ENDIF     
C     
C- READ THE NEXT OBSERVATION AND GO BACK UP TO THE 22001 LOOP TO CHECK  
C   THE CHRONOLOGICAL DAY RELATIVE TO THE COMPUTATIONAL DAY.
C     
        READ(DEV22,2265,IOSTAT=IOST22,ERR=22011,END=22031)  
     &       SFGYR,SFGMO, SFGDY, SFGHR, (IWORK1(100+IVBL),IVBL=1,22)    
        SFRED = SFRED + 1     
C     
        SFDAYC = JULIAN(SFGYR,SFGMO,SFGDY)
        CALL CHROND( SFGYR,SFDAYC,IWORK1(100) ) 
        IF (IWORK1(100) .LT. 0) THEN
          MESS = BLNK40 
          WRITE(MESS,140) 22  
          PRINT *,' CHROND ERROR AT SF OBS.',SFRED,' : ',MESS     
          STOP    
        ENDIF     
C     
        GO TO 22001     
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       ELSE IF (IWORK1(100) .GT. KDAY) THEN     
        GO TO 32000     
       ENDIF
      ENDIF 
      GO TO 32000 
C     
C-----------------------------------------------------------------------
C- PROCESSING CONTINUES HERE IF THERE WAS AN ERROR READING THE SFC. FILE
22011  MESS = BLNK40    
       WRITE(MESS,267)  IOST22
         PRINT *,' ERROR READING SF OBS.',SFRED,' : ',MESS  
         STOP     
C     
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C- PROCESSING CONTINUES HERE IF AN END OF FILE WAS ENCOUNTERED ON THE   
C-     SURFACE FILE     
22031  SFEOF  = .TRUE.  
C     
C=======================================================================
C- IF THERE ARE UPPER AIR DATA FOR THE DAY, COMPUTE THE POTENTIAL 
C   TEMPERATURE PROFILE THEN INTERPOLATE THE MINIMUM AND MAXIMUM  
C   (POTENTIAL) TEMPERATURES TO DETERMINE THE MORNING AND EVENING 
C   MIXING HEIGHTS.     
C     
C   UAOBS(1,ILEV,1) = PRESSURE FOR LEVEL 'ILEV' IN INTEGER FORMAT 
C   UAOBS(1,ILEV,2) = HEIGHT, IN METERS, OF LEVEL 'ILEV'    
C   UAOBS(1,ILEV,3) = TEMPERATURE FOR LEVEL 'ILEV' IN INTEGER FORMAT    
C   NOTE: MISSING UPPER AIR VALUES ARE DENOTED BY -9999     
C     
32000 IF ( UADAT.AND.SFDAT) THEN    
        DO 32010 ILEV = 1,UALEV(1)  
          IF ((UAOBS(1,ILEV,1) .GT. -9000) .AND.
     &        (UAOBS(1,ILEV,3) .GT. -9000)) THEN
            RATIO      = 1000.0/(UAOBS(1,ILEV,1)/10.0)
            TKELV      = (UAOBS(1,ILEV,3)/10.0) + 273.16    
            PTEM(ILEV) = TKELV*(RATIO**0.286)   
            Z(ILEV)    = UAOBS(1,ILEV,2)  
************ ADDED BY TOM WHITEHURST ***************************  
            WNDS(ILEV) = UAOBS(1,ILEV,6)  
****************************************************************  
          ELSE    
            PTEM(ILEV) = -9999.0    
          ENDIF   
32010   CONTINUE  
      ELSE  
*       CALL GREG (UAGYR,UADAYC-1,M1,D1)  
*       PRINT *,' NO DATA FOR DAY ',M1,D1 
        IAMMIX = -9000  
        IPMMIX = -9000  
      CALL ICHRND(KDAY,IWORK1(900),IWORK1(901)) 
      IF ((IWORK1(900).LT.0) .OR. (IWORK1(900).LT.0))THEN   
        MESS = BLNK40   
        WRITE(MESS,141) 40    
        PRINT *,' ERROR COMPUTING GREGORIAN DAY ON',KDAY,', ',MESS
        STOP
      ENDIF 
      CALL GREG (IWORK1(900),IWORK1(901),IWORK1(915),IWORK1(916)) 
*****   PRINT *,' NO DATA FOR DAY ',IWORK1(915),IWORK1(916) 
        GO TO 8991
      ENDIF 
***** 
***** 
            IF(IWORK1(20).GT.1)THEN 
              DO 32011 NLEV=1,UALEV(2)    
                 ZP(NLEV)=UAOBS(2,NLEV,2) 
                 WNDP(NLEV)=UAOBS(2,NLEV,6)     
32011         CONTINUE  
            ENDIF 
******
C     
C- INTERPOLATE THE MIN AND MAX TEMPERATURES TO FIND THE MIXING HTS.     
C   NOTE: AN ERROR CONDITION WILL OCCUR IF THE POTENTIAL TEMPERATURE    
C   AT THE FIRST LEVEL IS LESS THAN THE SURFACE POTENTIAL TEMPERATURE   
C   IN THE SOUNDING; AN ERROR CONDITION WILL OCCUR IF THE MAXIMUM 
C   POTENTIAL TEMPERATURE IS GREATER THAN ANY POTENTIAL TEMPERATURE     
C   IN THE SOUNDING, I.E., THE DRY ADIABAT DEFINED BY THE SURFACE 
C   TEMPERATURE DOES NOT INTERSECT THE SOUNDING.
C     
      IF ((TMIN.LT.9000.0).AND.(UAOBS(1,1,1).GT.-900.0)) THEN     
        TMIN5 = TMIN + FCTR   
        SFCTH = TMIN5*(1000.0/(UAOBS(1,1,1)/10.0))**0.286   
        CALL INTRP2(UALEV(1),PTEM,Z,SFCTH,AMMIX,IFLAG)
        IF (IFLAG .EQ. 1) THEN
******    PRINT *,' A.M. INTERPOLATION PROBLEM ON DAY ',KDAY
          IAMMIX = AMMIX*1    
        ENDIF     
      ELSE  
        IAMMIX = -9001  
        AMMIX = -9001.  
      ENDIF 
C     
      IF ((TMAX.GT.-900.0).AND.(UAOBS(1,1,1).GT.-900.0)) THEN     
      SFCTH = TMAX*(1000.0/(UAOBS(1,1,1)/10.0))**0.286
        CALL INTRP2(UALEV(1),PTEM,Z,SFCTH,PMMIX,IFLAG)
        IF (IFLAG .EQ. 1) THEN
******    PRINT *,' P.M. INTERPOLATION PROBLEM ON DAY ',KDAY
          IPMMIX = PMMIX*1    
        ENDIF     
      ELSE  
        PMMIX  = -9002. 
        IPMMIX = -9002  
      ENDIF 
C     
8990  IF (AMMIX .GT. 0.0) IAMMIX = AMMIX +0.5   
      IF (PMMIX .GT. 0.0) IPMMIX = PMMIX +0.5   
C=======================================================================
C- WRITE THE DATA TO THE OUTPUT FILE, CONVERTING THE CHONOLOGICAL DAY   
C   TO THE GREGORIAN DAY
C     
8991  CALL ICHRND(KDAY,IWORK1(900),IWORK1(901)) 
      IF ((IWORK1(900).LT.0) .OR. (IWORK1(900).LT.0))THEN   
        MESS = BLNK40   
        WRITE(MESS,141) 40    
        PRINT *,' ERROR COMPUTING GREGORIAN DAY ON',KDAY,', ',MESS
        STOP
      ENDIF 
      CALL GREG (IWORK1(900),IWORK1(901),IWORK1(915),IWORK1(916)) 
C     
      IF(RUNTYP.EQ.1)THEN     
      WRITE(40,4000) STNID,IWORK1(900),IWORK1(915),IWORK1(916),   
     &               IAMMIX,IPMMIX  
*AB  4000 FORMAT(A5,3(1X,I2),1X,I5,13X,I5)
 4000 FORMAT(A5,1X,I4,2(1X,I2),1X,I5,13X,I5)    
      ELSEIF(RUNTYP.EQ.2)THEN 
*     
*********************************************************************   
*********************************************************************   
*****SOFTWARE DEVELOPED BY TOM WHITEHURST 
***** 
*********************************************************************   
********CALL FUNCTION TO SET SEASON INDICATOR*******************  
      ISEASN=SETSEA(IWORK1(915))    
********FIND THE AM AND PM AVERAGE SUFRACE WINDS****************  
       IF(SWAS.GE.0.AND.ISWAN.GT.0)THEN   
           IAWSS=SWAS/ISWAN + .5    
       ELSE 
           IAWSS=-9999  
       ENDIF
       IF(SWPS.GE.0.AND.ISWPN.GT.0)THEN   
           IPWSS=SWPS/ISWPN + .5    
       ELSE 
           IPWSS=-9999  
       ENDIF
********SET THE TYPE********************************************  
       IF(IAMMIX.EQ.-9999.OR.IAMMIX.EQ.-9000.OR.IAMMIX.EQ.-9001)THEN    
        IATYP=4   
      ELSEIF(IAMMIX.EQ.-7777)THEN   
        IATYP=3   
      ELSEIF(ICAHP.EQ.1.OR.ICALP.GE.2)THEN
        IATYP=2   
      ELSE  
        IATYP=1   
      ENDIF 

      IF(IPMMIX.EQ.-9999.OR.IPMMIX.EQ.-9000.OR.IPMMIX.EQ.-9002)THEN     
        IPTYP=4   
      ELSEIF(IPMMIX.EQ.-8888)THEN   
        IPTYP=3   
      ELSEIF(ICPHP.EQ.1.OR.ICPLP.GE.2)THEN
        IPTYP=2   
      ELSE  
        IPTYP=1   
      ENDIF 
******************COMPUTE AVERAGE UA WINDS *********************  
      IF(IATYP.EQ.1.OR.IATYP.EQ.2)THEN    
        IF(IAMMIX.GT.150)THEN 
           CALL UALWND(UALEV(1),WNDS,Z,LYRWND,AMMIX,IFLAG)  
        ELSE
           LYRWND=IAWSS 
        ENDIF     
      ELSEIF(IATYP.EQ.3)THEN  
       LYRWND=SUM4WD(UALEV(1),WNDS,Z)     
       ELSEIF(IATYP.EQ.4)THEN 
       LYRWND=-9999     
       ENDIF
       IF(LYRWND.NE.-9999)THEN
         IAWSD=LYRWND + .5    
       ELSE 
         IAWSD=LYRWND   
       ENDIF
      IF(IPTYP.EQ.1.OR.IPTYP.EQ.2)THEN    

        IF(IPMMIX.GT.150)THEN 
           CALL UALWND(UALEV(2),WNDP,ZP,LYRWND,PMMIX,IFLAG) 
        ELSE
           LYRWND=IPWSS 
        ENDIF     
      ELSEIF(IPTYP.EQ.3)THEN  
       LYRWND=SUM4WD(UALEV(2),WNDP,ZP)    
       ELSEIF(IPTYP.EQ.4)THEN 
       LYRWND=-9999     
       ENDIF
       IF(LYRWND.NE.-9999)THEN
         IPWSD=LYRWND + .5    
       ELSE 
         IPWSD=LYRWND   
       ENDIF
***************************************************************** 
C     
      IF(UAEOF.OR.SFEOF)EOFFLG=.TRUE.     
      IF(KDAY.EQ.IWORK1(911))EOFFLG=.TRUE.
      CALL MIXRPT(STNU,IWORK1(900),IWORK1(915),ISEASN,IWORK1(916),
     &IATYP,IAMMIX,IAWSD,IAWSS,IPTYP,IPMMIX,IPWSD,IPWSS)    
       CALL MIXFMT(IWORK1(900),IWORK1(915),IWORK1(916),     
     &            IATYP,IAMMIX,IAWSD,IAWSS,IPTYP,IPMMIX,    
     &            IPWSD,IPWSS,EOFFLG,STNS,STNU,STNNS,STNNU) 
*     
       ISEASN=0   
       ICALP=ICNLP
       ICAHP=ICNHP
       ICPLP=0    
       ICPHP=0    
       ICNLP=0    
       ICNHP=0    
       SWAS=0     
       SWPS=0     
       ISWAN=0    
       ISWPN=0    
       IATYP=-9999
       IAWSD=-9999
       IAWSS=-9999
       IPTYP=-9999
       IPWSD=-9999
       IPWSS=-9999
***************************************************************   
         ENDIF    
C     
C- IF AN END-OF-FILE WAS ENCOUNTERD ON EITHER INPUT FILE, THEN    
C   PROCESSING STOPS    
C     
      IF (UAEOF .OR. SFEOF) THEN    
        PRINT *,' E-O-F ENCOUNTERED; PROCESSING STOPS'
        STOP
      ENDIF 
C     
C- FLUSH (REINITIALIZE) THE ARRAYS AND REINITIALIZE THE SCALARS   
      DO 39 K = 1,2     
       DO 38 I = 1,200  
        DO 49 J =1,6    
         UAOBS(K,I,J) = -9999 
   49   CONTINUE  
   38  CONTINUE   
   39 CONTINUE    
      IAMMIX = -9999    
      IPMMIX = -9999    
      TMIN =  9999.0    
      TMAX = -9999.0    
      DO 5000 I = 1,40  
        PTEM(I) = -9999.0     
        Z(I) = -9999.0  
        WNDS(I)= -9999.0
        ZP(I) = -9999.0 
        WNDP(I)= -9999.0
 5000 CONTINUE    
      UADAT = .FALSE.   
      UAWND = .FALSE.   
      SFDAT = .FALSE.   
      UALEV(1) = 0
      IWORK1(21) = 0    
      IWORK1(22) = 0    
C     
 9000 CONTINUE    
C     
C=======================================================================
C-FORMAT STATEMENTS     
C     
  140 FORMAT(' CHRONOLOGICAL DAY <0 ON UNIT',I2,', STOP')   
*AB  141 FORMAT(' YEAR OR JULIAN DAY <0 ON UNIT',I2,', STOP')     
  141 FORMAT(' YEAR OR JULIAN DAY <0 ON UNIT',I4,', STOP')  
  142 FORMAT(' INPUT CHRONOLOGICAL DAY <0, STOP')     
  255 FORMAT(' ERROR READING SDG. HDR., IOST=',I8)    
  256 FORMAT(' ERROR READING SDG. LEVEL, IOST=',I8)   
  257 FORMAT(' EOF READING LVL OF DATA, IOST=',I8)    
  267 FORMAT(' ERROR READING SFC OBS., IOST=',I8)     
*AB 1350 FORMAT(1X,4I2,I5,2(1X,I5)) 
 1350 FORMAT(1X,I4,3I2,I5,2(1X,I5)) 
 1355 FORMAT(6(1X,I5))  
C     
*AB  2265 FORMAT(1X,4I2,4(1X,I5),6(1X,I5.5),/,8X,4(1X,I5.5),8(1X,I5))   
 2265 FORMAT(1X,I4,3I2,4(1X,I5),6(1X,I5.5),/,8X,4(1X,I5.5),8(1X,I5))    
C     
      END   
C     
C=======================================================================
C***** SUBROUTINE UAMOVE(UALEV,IWORK1,IWORK2,UAOBS)   
C     
C   THIS SUBROUTINE MOVES THE UPPER AIR DATA FROM THE TEMPORARY   
C   ARRAY TO THE PERMANENT ARRAY    
C-----------------------------------------------------------------------
      SUBROUTINE UAMOVE(UALEV,IWORK1,IWORK2,UAOBS)    
       INTEGER LEV1, IVBL,UALEV(2),IWORK1(1000) 
     &,IWORK2(200,6),UAOBS(2,200,6) 
       UALEV(IWORK1(20))  = IWORK1(11)    
       DO 138 LEV1 = 1,UALEV(IWORK1(20))  
        DO 139 IVBL = 1,6     
         UAOBS(IWORK1(20),LEV1,IVBL) = IWORK2(LEV1,IVBL)    
  139   CONTINUE  
  138  CONTINUE   
      RETURN
C     
      END   
C     
C=======================================================================
      SUBROUTINE INTRP2(NTEMPS,TEMPS,Z,SFCT,VALUE,FLAG)     
C     
C---  THIS INTERPOLATION ROUTINE LINEARLY INTERPOLATES A PARAMETER TO   
C       A SPECIFIED LEVEL.  MISSING DATA ARE ALLOWED IN THE ARRAYS.     

C---  TO PERFORM THE INTERPOLATION FOR THE HOLZWORTH ESTIMATE OF THE    
C       MIXING HEIGHTS, THE TEMPS AND Z ARRAY FUNCTIONS ARE INTER-
C       CHANGED   
C     
C     
C     *** NOTE, UAMXL MUST BE SET TO BE EQUAL TO UAML 
C         AS THIS PROVIDES IN THE WORK ARRAYS PTEMP AND Z   
C         AS MANY LEVELS TO WORK WITH AS THERE ARE IN THE   
C         UA SOUNDINGS. 
C     
      REAL TEMPS(200),Z(200),SFCT,VALUE   
      INTEGER NTEMPS,FLAG,TOP,BOTTOM
      REAL MISS,MSSCHK,MLOW,MHIGH   
      DATA MISS/-9999./, MSSCHK/-9000.0/  
      DATA MLOW/-8888./, MHIGH/-7777./    
C     
C---  COMPARE THE SPECIFIED INTERPOLATION TEMPERATURE (SFCT) TO  THE    
C      LOWEST HEIGHT IN THE TEMPERATURE ARRAY (TEMPS(1)).   
C      IF SFCT < TEMPS(1), SET THE FLAG TO 1 AND RETURN,    
C      OTHERWISE CONTINUE.    
C     
      IF (SFCT .LT. TEMPS(1)) THEN  
        FLAG  = 1 
        VALUE = MLOW    
        RETURN    
      END IF
C     
C---  LOCATE THE LEVEL OF NON-MISSING HEIGHT AND TEMPERATURE DATA ABOVE 
C       THE INTERPOLATION HEIGHT, BEGINNING WITH THE BOTTOM LEVEL.
C     
      DO 200 I=2,NTEMPS 
        IF ((TEMPS(I).GE.SFCT) .AND. (TEMPS(I).GT.MSSCHK) .AND.   
     &      (Z(I).GT.MSSCHK)) THEN  
C     
C---  ALL NECESSARY DATA PRESENT ABOVE THE INTERPOLATION LEVEL;   
C      LOCATE THE LEVEL OF NON-MISSING HEIGHT AND TEMPERATURE DATA
C      BELOW THE INTERPOLATION HEIGHT.    
C     
            DO 210 II=I-1,1,-1
              IF ((TEMPS(II).GT.MSSCHK) .AND. (Z(II).GT.MSSCHK)) THEN   
C     
C---  ALL REQUIRED DATA ARE AVAILABLE - INTERPOLATE   
C     
                TOP=I   
                BOTTOM=II     
                TOTDIF = TEMPS(TOP) - TEMPS(BOTTOM)   
                PRTDIF = SFCT - TEMPS(BOTTOM)   
                FRAC   = PRTDIF/TOTDIF    
                PRMDIF = Z(TOP) - Z(BOTTOM)     
                VALUE  = Z(BOTTOM) + FRAC*PRMDIF
                FLAG   = 0    
                RETURN  
              ELSE
C     
C--- THE LOWER ARRAY HEIGHT OR TEMPERATURE IS MISSING 
C     
                CONTINUE
              ENDIF     
C     
210         CONTINUE    
C     
C--- NO LEVEL OF NON-MISSING DATA WAS LOCATED BELOW THE INTERPOLATION   
C      HEIGHT     
C     
            VALUE = MISS
            FLAG = 1    
            RETURN
        ELSE
C     
C--- THE ARRAY HEIGHT IS MISSING OR IS BELOW THE INTERPOLATION HEIGHT   
C     
          CONTINUE
        ENDIF     
C     
200   CONTINUE    
C     
C---  NO LEVEL OF NON-MISSING DATA WAS LOCATED ABOVE THE    
C       INTERPOLATION HEIGHT  
C     
      VALUE = MHIGH     
      FLAG = 1    
      RETURN
C     
      END   
C     
C=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
C=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
C  SUBROUTINE GMTLST ( YEAR,MONTH,DAY,HOUR,ZONE FACTOR )    
C     
C  PURPOSE  
C    COMPUTES LOCAL STANDARD TIME FROM GREENWICH MEAN TIME  
C     (LST FROM GMT).   
C     
C  CALLED BY:    GETSDG, MERGE
C     
C  CALLS TO:     GREG, JULIAN 
C     
C  VERSION DATE: 12 AUGUST 1987     
C     
C=======================================================================
      SUBROUTINE GMTLST( JYR,JMO,JDA,JHR,JZONE )
C     
      INTEGER JULIAN, JYR, JMO, JDA, JZONE, JULHR, JULDA, NDAY    
      INTEGER JHR, JULHRC, JDATE    
C     
C  JYR        = INPUT/OUTPUT YEAR (4 DIGITS)    
C  JMO        = INPUT/OUTPUT MONTH  
C  JDA        = INPUT/OUTPUT DAY    
C  JHR        = INPUT/OUTPUT HOUR   
C  JZONE      = INPUT TIME ZONES FROM GREENWICH 
C               + ==> WEST LONGITUDE
C               - ==> EAST LONGITUDE
C  JDATE      = TEMPORARY JULIAN DATE IN COMPUTATION  
C  JULDA      = JULIAN DATE   
C  JULHR      = JULIAN HOUR FROM JULIAN DATE    
C  JULHRC     = JULIAN HOUR CORRECTED FROM GMT TO LST 
C  NDAY       = NUMBER OF DAYS IN THE YEAR
C     
C-----------------------------------------------------------------------
C     
C---  COMPUTE JULIAN DAY
      JULDA = JULIAN(JYR,JMO,JDA)   
C     
C---  COMPUTE JULIAN HOUR     
      JULHR = (JULDA-1)*24. + JHR   
C     
C---  SUBTRACT THE TIME ZONE FACTOR TO GET CORRECTED JULIAN HOUR  
      JULHRC = JULHR - JZONE  
C     
C---  DETERMINE THE JULIAN DAY FROM THE CORRECTED JULIAN HOUR     
      JDATE = JULHRC/24 + 1   
C     
C---  DETERMINE THE NUMBER OF DAYS IN THE YEAR  
      IF(MOD(JYR,4) .EQ. 0)THEN     
        NDAY = 366
      ELSE  
        NDAY = 365
      ENDIF 
C     
C---  COMPUTE THE YEAR, MONTH, DAY AND HOUR FOR THE CORRECTED     
C---  JULIAN HOUR 
      IF(JULHRC .LT. 0)THEN   
        JYR   = JYR - 1 
        JMO   = 12
        JDA   = 31
        JHR   = 24 + JHR - JZONE    
      ELSE IF(JDATE .GT. NDAY)THEN  
        JDATE = JDATE - NDAY  
        JYR   = JYR + 1 
        CALL GREG(JYR,JDATE,JMO,JDA)
        JHR   = MOD(JULHRC,24)
      ELSE  
        JYR   = JYR     
        CALL GREG(JYR,JDATE,JMO,JDA)
        JHR   = MOD(JULHRC,24)
      ENDIF 
C     
      RETURN
      END   
C=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
C  SUBROUTINE GREG (YEAR,JULIAN DAY,MONTH,DAY)  
C     
C  PURPOSE  
C    THIS ROUTINE COMPUTES THE GREGORIAN MONTH AND DAY GIVEN THE YEAR   
C     AND JULIAN DAY    
C     
C  CALLED BY:     ANY ROUTINE THAT NEEDS TO CONVERT JULIAN TO     
C                 GREGORIAN DAY     
C     
C  CALLS TO:      -NONE-
C     
C  VERSION DATE:  12 AUGUST 1987    
C     
C=======================================================================
      SUBROUTINE GREG( YEAR,JDAY,MNTH,MDAY )    
C     
      INTEGER YEAR,JDAY,MNTH,MDAY,L,J     
C     
C      YEAR = CALENDAR YEAR   
C      JDAY = DAY OF YEAR -- 1, 365, OR 366     
C      MNTH = MONTH OF YEAR -- 1, 12
C      MDAY = DAY OF MONTH -- 1, 31 
C     
C-----------------------------------------------------------------------
      L=365 
      IF(MOD(YEAR,4).EQ.  0) L=366  
      IF(MOD(YEAR,100).EQ.0) L=365  
      IF(MOD(YEAR,400).EQ.0) L=366  
      J=MOD(JDAY+305,L) 
      J=MOD(J,153)/61+(J/153)*2+J   
      MNTH=MOD(J/31+2,12)+1   
      MDAY=MOD(J,31)+1  
      RETURN
      END   
C     
C=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
C  SUBROUTINE JULIAN ( YEAR,MONTH,DAY )           (A SUBFUNCTION) 
C     
C  PURPOSE  
C     ROUTINE CONVERTS GREGORIAN MONTH AND DAY  
C     TO JULIAN DAY, FOR GIVEN YEAR.
C     
C        LIMITATIONS:  ROUTINE WORKS ONLY FOR YEARS 1901 - 1999   
C                      INCLUSION OF TEST FOR 100 AND 400 YEAR     
C                      CORRECTIONS COULD EXTEND THIS ROUTINE
C                      THOSE COUNTRIES USING OUR CURRENT    
C                      CALENDAR.  THIS WOULD REQUIRE YEARS TO     
C                      BE EXPRESSED FULLY (E.G. 1936 NOT 36).     
C                      CURRENT CALENDAR CAME INTO     
C                      OFFICIAL USE IN THE WORLD IN 1582, IN
C                      GREAT BRITAIN (AND ASSOCIATD COLONIES
C                      IN 1752.).   
C     
C  CALLED BY:     ANY ROUTINE THAT NEEDS TO CONVERT GREGORIAN     
C                 TO JULIAN DAY     
C     
C  CALLS TO:      -NONE-
C     
C  VERSION DATE:  01 SEPTEMBER 1987 
C     
C=======================================================================
C     
        FUNCTION JULIAN( YR1,MO1,DY1 )    
C     
C    LOCAL VARIABLES    
C     
        INTEGER YR1,MO1,DY1,JULIAN,IWORK1,IWORK2,IWORK3     
C        YR1        LAST 2-DIGITS CALENDAR YEAR 
C        MO1        MONTH OF YEAR (1-12)  
C        DY1        DAY OF MONTH (1-31)   
C-----------------------------------------------------------------------
C     
        IF( YR1.LE.0 .OR. MO1.LE.0 .OR. DY1.LE.0 ) THEN     
         JULIAN = 0     
        ELSE
C     
         IWORK1 = MOD( (MO1+9),12 ) 
         IWORK2 = (IWORK1*153+2)/5 + DY1 + 58   
C     
         IF( MOD(YR1,4).EQ.0 ) THEN 
          IWORK3 = 366  
          IWORK2 = IWORK2 + 1 
         ELSE     
          IWORK3 = 365  
         END IF   
C     
         JULIAN = 1 + MOD( IWORK2,IWORK3 )
C     
        END IF    
C     
        RETURN    
        END 
C     
C     
C=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
C  SUBROUTINE CHROND ( IYR,IJD,IYDY )     
C     
C  PURPOSE  
C    THIS SUBROUTINE CONVERTS DATES BETWEEN 1900 AND 2099 TO A    
C     COMMON SEQUENTIAL DAY.  THE ZERO POINT IS JANUARY 1, 1900   
C     I.E., 1 IS JANUARY 1,1900.    
C     
C        ARGUMENT LIST: 
C          INPUT:   IYR  - YEAR (4 DIGITS)
C                   IJD  - JULIAN DAY (3 DIGITS)
C          OUTPUT:  IYDY - DATE CONVERTED TO SEQUENTIAL DAY 
C                          SINCE JANUARY 1, 1900
C     
C  CALLED BY:     UAEXT,GETSDG,GETMIX,SFEXT,GETSFC,MERGE    
C     
C  CALLS TO:      ERROR 
C     
C  VERSION DATE:  17 MARCH 1988     
C     
C=======================================================================
      SUBROUTINE CHROND( IYR,IJD,IYDY )   
C     
        INTEGER NDAY,IYDY,IYR,IJD,ICNT,TEST     
      CHARACTER MESS*40,BLNK40*40   
        DATA BLNK40/'                                        '/   
C *** CHECK FOR INVALID DATA  
C     
*AB  SUBTRACT ZERO POINT (1900)     
        IYR = IYR - 1900
*AB        IF ( IYR.LT.0 .OR. IYR.GE.100 ) THEN 
        IF ( IYR.LT.0 .OR. IYR.GE.200 ) THEN    
         MESS = BLNK40  
         WRITE(MESS,1000) IYR 
         PRINT*,MESS    
         IYDY = -9999   
         RETURN   
        END IF    
C     
        IF ( IJD.LE.0 .OR. IJD.GE.367 ) THEN    
         MESS = BLNK40  
         WRITE(MESS,1010) IJD 
         PRINT*,MESS    
         IYDY = -9999   
         RETURN   
        END IF    
C     
C *** TEST FOR LEAP YEAR AND DETERMINE THE NUMBER OF DAYS IN THE YEAR   
        TEST = MOD(IYR,4)     
        IF( TEST.EQ.0 .OR.IYR.EQ.0) THEN  
         NDAY = 366     
        ELSE
         NDAY = 365     
        END IF    
C     
C *** CHECK FOR A VALID JULIAN DAY  
C     
        IF ( IJD .GT. NDAY ) THEN   
         MESS = BLNK40  
         WRITE(MESS,1030) IJD 
         PRINT*,MESS    
         IYDY = -9999   
         RETURN   
        END IF    
C     
C *** COMPUTE THE CHONOLOGICAL DAY  
C     
        ICNT = (IYR+3)/4
        IYDY = 366*ICNT + 365*(IYR-ICNT) + IJD  
*AB   ADD BACK 1900     
        IYR = IYR + 1900
C     
C     
        RETURN    
C-----------------------------------------------------------------------
C *** FORMAT STATEMENTS.
C     
*AB 1000    FORMAT(' INPUT YEAR (',I4,') NOT IN RANGE OF 00-99')  
C     
*AB 1111    FORMAT(' IYR  ',I4,'  IJD  ',I3,'  IYDY ',I6)   
1000    FORMAT(' INPUT YEAR (',I4,') NOT IN RANGE OF 00-199')     
1010    FORMAT(' JULIAN DAY (',I4,') NOT IN RANGE OF 1-366')
1030    FORMAT(' NOT A LEAP YEAR: INVALID JULIAN DAY',I4)   
        END 
C     
C=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
C  SUBROUTINE ICHRND(IYDY,IYR,IJD)  
C     
C  PURPOSE  
C    THIS SUBROUTINE CONVERTS A COMMON SEQUENTIAL DAY TO A DATE   
C     BETWEEN 1900 AND 1999.  THE ZERO POINT IS 1899, 365, 24,    
C     I.E., DAY 1 IS THE FIRST DAY OF 1900.     
C     
C        ARGUMENT LIST: 
C          INPUT:   IYDY - COMMON SEQUENTIAL DAY
*AB INCREASE YR              OUTPUT:  IYR  - YEAR (2 DIGITS)
C          OUTPUT:  IYR  - YEAR (4  
C                   IJD  - JULIAN DAY (3 DIGITS)
C     
C  CALLED BY:     MERGE 
C     
C  CALLS TO:      ERROR 
C     
C  VERSION DATE:  17 MARCH 1988     
C     
C=======================================================================
      SUBROUTINE ICHRND( IYDY,IYR,IJD )   
C     
      INTEGER IYDY,IYR,IJD    
      CHARACTER MESS*40,BLNK40*40   
        DATA BLNK40/'                                        '/   
C      IYDYO   = OLD SEQUENTIAL DAY 
C      IYDYN   = NEW SEQUENTIAL DAY 
C-----------------------------------------------------------------------
C *** CHECK FOR INVALID DATA  
C     
*AB  SUBTRACT 1900 FROM 4 DIGIT YEAR.     
       IYR = IYR - 1900 
      IF ((IYDY.LT.1) .OR. (IYDY.GT.36525)) THEN
       MESS = BLNK40    
       WRITE(MESS,1000) IYDY  
       PRINT*,MESS
*AB       IYR = -99     
       IYR = -199 
*AB       IJD = -99     
       IJD = -199 
       RETURN     
      ENDIF 
C     
C *** CONVERT COMMON SEQUENTIAL DAY TO 2-DIGIT YEAR AND JULIAN DAY
C     
      IYDYO = IYDY
      IYR = 0     
    1 IF(MOD(IYR,4) .EQ. 0) THEN    
       NDAY = 366 
      ELSE  
       NDAY = 365 
      ENDIF 
C     
      IYDYN = IYDYO - NDAY    
C     
C *** WHEN THE CHRONOLOGICAL DAY BECOMES NEGATIVE, THE YEAR AND JULIAN  
C      DAY HAVE BEEN DETERMINED     
C     
      IF(IYDYN .GT. 0) THEN   
       IYR = IYR + 1    
       IYDYO = IYDYN    
       GO TO 1    
      ELSE IF(IYDYN .LE. 0) THEN    
       IJD = IYDYO
      ENDIF 
C     
C *** CHECK FOR VALID YEAR AND JULIAN DAY 
C     
*AB      IF(IYR .GT. 99) THEN 
      IF(IYR .GT. 199) THEN   
       GO TO 999  
      ELSE  
       IF(MOD(IYR,4) .EQ. 0) THEN   
        NDAY = 366
       ELSE 
        NDAY = 365
       ENDIF
C     
       IF((IJD .LT. 0) .OR. (IJD.GT.NDAY)) THEN 
        GO TO 999 
       ENDIF
      ENDIF 
C     
*AB ADD 1900 BACK TO YEAR, IT WAS SUBTRACTED AT BEGINNING   
         IYR = IYR + 1900     
      RETURN
C     
  999 MESS = BLNK40     
      WRITE(MESS,2000)  
      PRINT*,MESS 
*AB      IYR = -99
      IYR = -199  
*AB      IJD = -99
      IJD = -199  
      RETURN
C     
C-----------------------------------------------------------------------
C *** FORMAT STATEMENTS.
C     
 1000 FORMAT(' CHRON. DAY NOT IN 1-36525 RANGE ',I7)  
 2000 FORMAT(' YEAR/JULIAN DAY NOT VALID; SET TO -99' )     
C     
      END   
*********************************************************************   
*********************************************************************   
*****SOFTWARE DEVELOPED BY TOM WHITEHURST 
*****PRODUCES TAPE AND DISK FORMATS.  DISK FILE CAN BE DOWN 
*****LOADED TO A FOLPPY DISK. 
***** 
*********************************************************************   
       SUBROUTINE MIXRPT(STN,YR,MO,SEA,DA,ATYP,AMD,AWSD,AWSS,     
     &PTYP,PMD,PWSD,PWSS)     
       IMPLICIT INTEGER(A-Z)  
       CHARACTER*5 STN  
       CHARACTER*4 SP   
*AB       CHARACTER*57 FMT,FMS
       CHARACTER*61 FMT,FMS   
*AB       DATA    FMT/'(A5,2I2.2,I1,I2.2,I1,I4.4,I3.3,I3.3,I1,I4.4,I3.3,I3.   
       DATA    FMT/'(A5,I4.4,I2.2,I1,I2.2,I1,I4.4,I3.3,I3.3,I1,I4.4,I3.3,     
     &I3.3,A4)'/  
*AB       DATA    FMS/'(A5,2I2.2,I1,I2.2,I1,I4.4,I3.3,I3.3,I1,I4.4,I3.3,I3.   
       DATA    FMS/'(A5,I4.4,I2.2,I1,I2.2,I1,I4.4,I3.3,I3.3,I1,I4.4,I3.3,     
     &I3.3,A4)'/  
       DATA SP/'    '/  
       FMT=FMS    
       MAMD=AMD   
       MPMD=PMD   
       MAWSD=AWSD 
       MPWSD=PWSD 
       MAWSS=AWSS 
       MPWSS=PWSS 
       IF(ATYP.EQ.4)THEN
         MAMD=0   
*AB          FMT(25:25)='0'   
          FMT(29:29)='0'
         MAWSD=0  
*AB     FMT(30:30)='0'  
          FMT(34:34)='0'
         MAWSS=0  
*AB     FMT(35:35)='0'  
          FMT(39:39)='0'
       ELSEIF(ATYP.EQ.3)THEN  
         MAMD=0   
*AB     FMT(25:25)='0'  
          FMT(29:29)='0'
          IF(AWSD.LT.0)THEN   
*AB        FMT(30:30)='0'     
             FMT(34:34)='0'   
            MAWSD=0     
          ENDIF   
          IF(AWSS.LT.0)THEN   
            MAWSS=0     
*AB        FMT(35:35)='0'     
             FMT(39:39)='0'   
          ENDIF   
       ELSE 
          IF(AMD.LT.0)THEN    
            MAMD=0
*AB        FMT(25:25)='0'     
             FMT(29:29)='0'   
          ENDIF   
          IF(AWSS.LT.0)THEN   
            MAWSS=0     
*AB        FMT(35:35)='0'     
             FMT(39:39)='0'   
          ENDIF   
          IF(AWSD.LT.0)THEN   
*AB       FMT(30:30)='0'
            FMT(34:34)='0'    
           MAWSD=0
          ENDIF   
        ENDIF     
       IF(PTYP.EQ.4)THEN
         MPMD=0   
*AB     FMT(43:43)='0'  
          FMT(47:47)='0'
         MPWSD=0  
*AB     FMT(48:48)='0'  
          FMT(52:52)='0'
         MPWSS=0  
*AB     FMT(53:53)='0'  
          FMT(57:57)='0'
       ELSEIF(PTYP.EQ.3)THEN  
         MPMD=0   
*AB     FMT(43:43)='0'  
          FMT(47:47)='0'
          IF(PWSD.LT.0)THEN   
            MPWSD=0     
*AB        FMT(48:48)='0'     
             FMT(52:52)='0'   
          ENDIF   
          IF(PWSS.LT.0)THEN   
            MPWSS=0     
*AB        FMT(53:53)='0'     
             FMT(57:57)='0'   
          ENDIF   
       ELSE 
          IF(PMD.LT.0)THEN    
           MPMD=0 
*AB       FMT(43:43)='0'
            FMT(47:47)='0'    
         ENDIF    
         IF(PWSS.LT.0)THEN    
           MPWSS=0
*AB       FMT(53:53)='0'
            FMT(57:57)='0'    
         ENDIF    
         IF(PWSD.LT.0)THEN    
           MPWSD=0
*AB       FMT(48:48)='0'
            FMT(52:52)='0'    
         ENDIF    
         ENDIF    
       WRITE(23,FMT)STN,YR,MO,SEA,DA,ATYP,MAMD,MAWSD,MAWSS,PTYP   
     &,MPMD,MPWSD,MPWSS,SP    
       WRITE(24,FMT)STN,YR,MO,SEA,DA,ATYP,MAMD,MAWSD,MAWSS,PTYP   
     &,MPMD,MPWSD,MPWSS,SP    
       RETURN     
       END  
*********************************************************************   
*********************************************************************   
*****SOFTWARE DEVELOPED BY TOM WHITEHURST 
*****THIS ROUTINE IS USED TO SET THE TYPE OF REPORT BASED ON
***** PRECIPITATION.  THE PRECIPATION CODES COME FROM THE TD3280  
***** MANUAL.     
* IF PRECIPITATION OCCURRED AT TWO OR MORE OBSERVATION TIMES WITH 
* THE INTENSITY OF LIGHT, OR AT ONE OBSERVATION WITH THE INTENSITY
*OF MODERATE OR HEAVY DURING THE HOURS 1000-2100 LST FOR THE
*AFTERNOON OR DURING THE HOURS 2200(PREVIOUS DAY )-0900 LST FOR   
*THE MORNING.  THE TYPE IS SET TO 2.
*****IPWX = CURRENT WEATHER   
*****CURHR = HOUR 
*****ICALP,ICPLP = AM AND PM LIGHT PRECIP INDICATORS  
*****ICAHP,ICPHP = AM AND PM HEAVY PRECIP INDICATORS  
*****ICNLP,ICNHP = NIGHT BEFORE LIGHT & HEAVY PRECIP INDICATORS   
*********************************************************************   
       SUBROUTINE SETTYP(IPWX,CURHR,ICALP,ICAHP,ICPLP,ICPHP,
     &ICNLP,ICNHP)
       IMPLICIT INTEGER(A-Z)  
       CHARACTER*2 LTPCP(16),HYPCP(32)    
       CHARACTER*5 CPWX 
       DATA LTPCP/'20','23','26','30','33','36','40','43','46',   
     &'50','53','56','60','63','66','90'/ 

       DATA HYPCP/'21','22','24','25','27','28','31','32','34',   
     &'35','37','38','41','42','44','45','47','48','51','52','54',
     &'55','57','58','61','62','64','65','67','68','91','92'/     
       IF(IPWX.EQ.0)THEN
         RETURN   
        ENDIF     
       WRITE(CPWX,'(I5)')IPWX 

       DO 10 I=1,32     
       IF(CPWX(2:3).EQ.HYPCP(I).OR.CPWX(4:5).EQ.HYPCP(I))THEN     
         IF(CURHR.GE.10.AND.CURHR.LE.21)THEN    
            ICPHP=1     
         ELSEIF(CURHR.GE.22.AND.CURHR.LE.23)THEN
            ICNHP=1     
         ELSE     
            ICAHP=1     
         ENDIF    
         RETURN   
        ENDIF     
 10     CONTINUE  
       DO 12 I=1,16     
       IF(CPWX(2:3).EQ.LTPCP(I).OR.CPWX(4:5).EQ.LTPCP(I))THEN     
         IF(CURHR.GE.10.AND.CURHR.LE.21)THEN    
            ICPLP=ICPLP+1     
         ELSEIF(CURHR.GE.22.AND.CURHR.LE.23)THEN
            ICNLP=ICNLP+1     
         ELSE     
            ICALP=ICALP+1     
         ENDIF    
*        PRINT*,CPWX(2:3)     
*        PRINT*,ICALP,ICPLP   
*        PRINT*,CURHR   
         RETURN   
        ENDIF     
12      CONTINUE  
       RETURN     
       END  
*********************************************************************   
*****SOFTWARE DEVELOPED BY TOM WHITEHURST 
***** SET THE SEASON INDICATOR BASED ON MONTH   
***** 
***** 
*********************************************************************   
       INTEGER FUNCTION SETSEA(MO)  
       INTEGER MO 
       IF(MO.EQ.12.OR.MO.EQ.1.OR.MO.EQ.2)THEN   
       SETSEA=1   
       ELSEIF(MO.EQ.3.OR.MO.EQ.4.OR.MO.EQ.5)THEN
       SETSEA=2   
       ELSEIF(MO.EQ.6.OR.MO.EQ.7.OR.MO.EQ.8)THEN
       SETSEA=3   
       ELSEIF(MO.EQ.9.OR.MO.EQ.10.OR.MO.EQ.11)THEN    
       SETSEA=4   
       ENDIF
       RETURN     
       END  
*********************************************************************   
*********************************************************************   
*********************************************************************   
*****SOFTWARE DEVELOPED BY TOM WHITEHURST 
***** 
*********************************************************************   
      SUBROUTINE UALWND(NLVLS,WNDS,Z,LYRWND,VALUE,FLAG)     
C     
C---  THIS INTERPOLATION ROUTINE LINEARLY INTERPOLATES WIND SPEED  TO   
C       A SPECIFIED LEVEL.  MISSING DATA ARE ALLOWED IN THE ARRAYS.     

C---  TO PERFORM THE INTERPOLATION FOR THE HOLZWORTH ESTIMATE OF THE    
C       MIXING HEIGHTS, THE WNDSS AND Z ARRAY FUNCTIONS ARE INT  -
C       CHANGED   
C     
C     
      REAL WNDS(200),Z(200),INTWND,VALUE,AVGWND,LYRWND
      INTEGER NLVLS,FLAG,TOP,BOTTOM 
      REAL MSSCHK 
C     REAL MISS,MSSCHK,MLOW,MHIGH   
      DATA MSSCHK/-9999./     
C---  LOCATE THE LEVEL OF NON-MISSING HEIGHT AND TEMPERATURE DATA ABOVE 
C       THE INTERPOLATION HEIGHT, BEGINNING WITH THE BOTTOM LEVEL.
C     
      DO 200 I=2,NLVLS  
        IF ((Z(I).GE.VALUE) .AND. (Z(I).NE.MSSCHK) .AND.    
     &      (WNDS(I).NE.MSSCHK)) THEN     
C     
C---  ALL NECESSARY DATA PRESENT ABOVE THE INTERPOLATION LEVEL;   
C      LOCATE THE LEVEL OF NON-MISSING HEIGHT AND TEMPERATURE DATA
C      BELOW THE INTERPOLATION HEIGHT.    
C     
            DO 210 II=I-1,1,-1
              IF ((WNDS(II).NE.MSSCHK) .AND. (Z(II).NE.MSSCHK)) THEN    
C     
C---  ALL REQUIRED DATA ARE AVAILABLE - INTERPOLATE   
C     
                TOP=I   
                BOTTOM=II     
                TOTDIF = Z(TOP) - Z(BOTTOM)     
                PRTDIF = VALUE - Z(BOTTOM)
                FRAC   = PRTDIF/TOTDIF    
                PRMDIF = WNDS(TOP) - WNDS(BOTTOM)     
                INTWND  = WNDS(BOTTOM) + FRAC*PRMDIF  
                FLAG   = 0    
                LYRWND = AVGWND(BOTTOM,WNDS,INTWND)   
                RETURN  
              ELSE
C     
C--- THE LOWER ARRAY HEIGHT OR WIND SPEED  IS MISSING 
C     
                CONTINUE
              ENDIF     
C     
210         CONTINUE    
C     
C--- NO LEVEL OF NON-MISSING DATA WAS LOCATED BELOW THE INTERPOLATION   
C      HEIGHT     
C     
            LYRWND = MSSCHK   
            FLAG = 1    
            RETURN
        ELSE
C     
C--- THE ARRAY HEIGHT IS MISSING OR IS BELOW THE INTERPOLATION HEIGHT   
C     
          CONTINUE
        ENDIF     
C     
200   CONTINUE    
C     
C---  NO LEVEL OF NON-MISSING DATA WAS LOCATED ABOVE THE    
C       INTERPOLATION HEIGHT  
C     
      LYRWND = MSSCHK   
      FLAG = 1    
      RETURN
C     
      END   
C     
*********************************************************************   
*********************************************************************   
*****SOFTWARE DEVELOPED BY TOM WHITEHURST 
***** AVERAGE THE WIND SPEEDS IN THE MIXING LAYER     
***** RETURN THE AVERAGE WIND SPEED 
***** 
*********************************************************************   
      REAL FUNCTION AVGWND(BOTTOM,WNDS,INTWND)  
       INTEGER BOTTOM, NWNDS  
       REAL WNDS,INTWND,SUMWND
       DIMENSION WNDS(BOTTOM) 
       SUMWND=0   
       NWNDS=0    
       DO 10 I=1,BOTTOM 
          IF(WNDS(I).GE.0.AND.WNDS(I).NE.-9999)THEN   
             SUMWND=SUMWND+WNDS(I)  
             NWNDS=NWNDS+1    
          ENDIF   
  10   CONTINUE   
       IF(SUMWND.EQ.0.AND.NWNDS.EQ.0)THEN 
          AVGWND=-9999  
       ELSE 
       SUMWND=SUMWND+INTWND   
       NWNDS=NWNDS+1    
       AVGWND=SUMWND/NWNDS    
       ENDIF
*      PRINT*,AVGWND    
*      PRINT*,INTWND    
       RETURN     
       END  
*********************************************************************   
*********************************************************************   
*****SOFTWARE DEVELOPED BY TOM WHITEHURST 
***** 
******AVERAGE WIND SPEEDS SURFACE TO 4000 FEET  
***** DONE WHEN THE TYPE IS "3" COLD ADVECTION  
*********************************************************************   
      REAL FUNCTION SUM4WD(NLVLS,WNDS,Z)  
       INTEGER NLVLS    
       REAL WNDS,Z,SUMWND     
       DIMENSION WNDS(NLVLS),Z(NLVLS)     
       SUM4WD=0   
       NWNDS=0    
       SUMWND=0   
       DO 10 I=1,NLVLS  
       IF(Z(I).NE.-9999..AND.WNDS(I).NE.-9999.)THEN   
          IF(Z(I).GE.0.AND.Z(I).LE.4000.AND.WNDS(I).GE.0)THEN     
             SUMWND=SUMWND+WNDS(I)  
             NWNDS=NWNDS+1    
          ENDIF   
       ENDIF
10     CONTINUE   
       IF(SUMWND.EQ.0.AND.NWNDS.EQ.0)THEN 
          SUM4WD=-9999  
       ELSE 
          SUM4WD=SUMWND/NWNDS 
       ENDIF
*      PRINT*,SUM4WD    
       RETURN     
       END  
*********************************************************************   
** CREATE THE PRINT FILE
*********************************************************************   
       SUBROUTINE MIXFMT (IYR,IMO,IDA,IATYP,IAMD,IAWSD,IAWSS,     
     &IPTYP,IPMD,IPWSD,IPWSS,EOFFLG,STNS,STNU,STNNS,STNNU)  
       CHARACTER*1 SP   
       CHARACTER*4,CAWSD,CAWSS,CPWSD,CPWSS
       CHARACTER*5 STNS,STNU  
       CHARACTER*30 STNNS,STNNU     
*********************************************************************   
       CHARACTER*64 HRD1,HRD2,HRD3  
       CHARACTER*50 FOOT5     
       CHARACTER*44 TITL1     
       CHARACTER*26 FOOT1,FOOT2,FOOT3,FOOT4,FOOT6     
       CHARACTER*18 TITL2,TITL3     
       LOGICAL EOFFLG   
         DATA SP/' '/   
      DATA KOUNT/0/     
      DATA TITL1/'DAILY MIXING HEIGHTS AND AVERAGE WIND SPEEDS'/  
      DATA TITL2/' U/A DATA STA. NO.'/    
      DATA TITL3/' SFC DATA STA. NO.'/    
*AB      DATA HRD3/'YR  MO  DAY  T   HGT   LYR   SFC             T   HGT
      DATA HRD3/'YEAR  MO  DAY  T   HGT   LYR   SFC           T   HGT   
     &LYR   SFC'/ 
      DATA HRD2/'                   MXG  WIND  WIND               MXG  W
     &IND  WIND'/ 
      DATA HRD1/'                    MORNING                       AFTER
     &NOON     '/ 
      DATA FOOT1/'MIXING HEIGHT IN METERS   '/  
      DATA FOOT2/'WIND SPEED IN METERS/SEC  '/  
      DATA FOOT3/'TYPE  1 = NO PRECIPITATION'/  
      DATA FOOT4/'TYPE  2 = PRECIPITATION   '/  
      DATA FOOT5/'TYPE  3 = SELECTED TEMPERATURE LESS THAN RAOB TEMP'/  
      DATA FOOT6/'TYPE  4 = MISSING         '/  
      KOUNT=KOUNT+1     
      IF(IATYP.EQ.4)THEN
         CAWSD='    '   
         CAWSS='    '   
         MAMD=0   
      ELSE  
          IF(IAMD.LT.0)THEN   
             MAMD=0     
          ELSE    
             MAMD=IAMD  
          ENDIF   
          IF(IAWSD.LT.0)THEN  
             CAWSD='   '
          ELSE    
             AWSD=IAWSD/10.0  
             WRITE(CAWSD,'(F4.1)')AWSD    
          ENDIF   
          IF(IAWSS.LT.0)THEN  
             CAWSS='   '
          ELSE    
             AWSS=IAWSS/10.0  
             WRITE(CAWSS,'(F4.1)')AWSS    
          ENDIF   
          ENDIF   
      IF(IPTYP.EQ.4)THEN
         CPWSD='    '   
         CPWSS='    '   
         MPMD=0   
      ELSE  
          IF(IPMD.LT.0)THEN   
             MPMD=0     
          ELSE    
             MPMD=IPMD  
          ENDIF   
          IF(IPWSD.LT.0)THEN  
             CPWSD='   '
          ELSE    
             PWSD=IPWSD/10.0  
             WRITE(CPWSD,'(F4.1)')PWSD    
          ENDIF   
          IF(IPWSS.LT.0)THEN  
             CPWSS='   '
          ELSE    
             PWSS=IPWSS/10.0  
             WRITE(CPWSS,'(F4.1)')PWSS    
          ENDIF   
      ENDIF 
1     IF(KOUNT.EQ.1)THEN
         LMO=IMO  
         WRITE(21,*)SP  
         WRITE(21,111)TITL1   
111      FORMAT('1',16X,A44)  
         WRITE(21,*)SP  
         WRITE(21,110)TITL2,STNU,STNNU    
         WRITE(21,*)SP  
         WRITE(21,110)TITL3,STNS,STNNS    
         WRITE(21,*)SP  
         WRITE(21,*)SP,HRD1   
         WRITE(21,*)SP,HRD2   
         WRITE(21,*)SP,HRD3   
      ENDIF 
      IF(IMO.EQ.LMO)THEN
      WRITE(21,100)IYR,IMO,IDA,IATYP,MAMD,CAWSD,CAWSS,IPTYP,MPMD, 
     &CPWSD,CPWSS 
      ELSE  
         WRITE(21,*)SP  
      WRITE(21,*)SP,FOOT1     
         WRITE(21,*)SP  
      WRITE(21,*)SP,FOOT2     
         WRITE(21,*)SP  
      WRITE(21,*)SP,FOOT3     
         WRITE(21,*)SP  
      WRITE(21,*)SP,FOOT4     
         WRITE(21,*)SP  
      WRITE(21,*)SP,FOOT5     
         WRITE(21,*)SP  
      WRITE(21,*)SP,FOOT6     
      KOUNT=1     
      GO TO 1     
      ENDIF 
      IF(EOFFLG)THEN    
         WRITE(21,*)SP  
      WRITE(21,*)SP,FOOT1     
         WRITE(21,*)SP  
      WRITE(21,*)SP,FOOT2     
         WRITE(21,*)SP  
      WRITE(21,*)SP,FOOT3     
         WRITE(21,*)SP  
      WRITE(21,*)SP,FOOT4     
         WRITE(21,*)SP  
      WRITE(21,*)SP,FOOT5     
         WRITE(21,*)SP  
      WRITE(21,*)SP,FOOT6     
      ENDIF 
*     
110   FORMAT(A17,5X,A5,5X,A30)
*AB 100   FORMAT(2X,I2.2,2X,I2.2,3X,I2.2,2X,I1,2X,I4.0,2X,A4,2X,A4,13X  
100   FORMAT(2X,I4.4,2X,I2.2,3X,I2.2,2X,I1,2X,I4.0,2X,A4,2X,A4,11X
     &,I1,2X,I4.0,2X,A4,2X,A4)
      RETURN
      END   
