      SUBROUTINE ZSTATS(DA,SI,OBUF,SI2,J,KLEV)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    ZSTATS      COMPUTE AND PRINT STATS ON HEIGHTS
C   PRGMMR: DIMEGO           ORG: W/NMC22    DATE: 86-03-31
C
C ABSTRACT: COMPUTES STATISTICS ON THE COMPARISON OF ORIGINAL Z
C   PROFILES FROM FIRST-GUESS WITH NEW PROFILES OF Z. ENTRY POINT
C   DSTATS PRINTS OUT A TABLE OF STATISTICS FOR THE COMPARISON AT
C   A FIXED SET OF 16 PRESSURES. BOTH INPUT PROFILES ARE IN SIGMA
C   (GUESS & RAFS) AND MUST BE INTERPOLATED TO THE 16 PRESSURES.
C   QUADRATIC IN LN P USED ALA THE LFM PROCEDURE.
C
C PROGRAM HISTORY LOG:
C   86-03-31  DIMEGO
C   88-06-30  DIMEGO      HALF PRECISION / PL1 / FTN200
C
C USAGE:    CALL ZSTATS(DA,SI,OBUF,SI2,J,KLEV)
C   INPUT ARGUMENT LIST:
C     DA       - MULTI-VARIABLE ARRAY OF GUESS STRIPS AT KLEV LEVELS
C              - FIRST     VARIABLE IS GEOPOTENTIAL AT INTERFACES
C              - SECOND    VARIABLE IS TEMPERATURE AT INTERFACES
C              - THIRD     VARIABLE IS U COMPONENT OF WIND
C              - FOURTH    VARIABLE IS V COMPONENT OF WIND
C              - FIFTH     VARIABLE IS SPECIFIC HUMIDITY
C     SI       - INTERFACE SIGMAS FOR THE GUESS VERTICAL STRUCTURE
C     OBUF     - MULTI-VARIABLE ARRAY OF GUESS STRIPS AFTER THEY HAVE
C              - BEEN INTERPOLATED TO RAFS SIGMA STRUCTURE
C     SI2      - INTERFACE SIGMAS FOR THE RAFS VERTICAL STRUCTURE
C     J        - INDEX OF THE STRIP BEING PROCESSED
C     KLEV     - NUMBER OF LAYERS IN THE GUESS VERTICAL STRUCTURE
C
C   OUTPUT FILES:
C     UNIT6    - DIAGNOSTIC PRINTOUT OF STATISTICS THROUGH
C              - CALL TO DSTATS    NO ARGUMENTS
C
C REMARKS:   NO CAVEATS, OTHER HELPFUL HINTS NOR INFORMATION
C
C ATTRIBUTES:
C   LANGUAGE: STANDARD FORTRAN
C   MACHINE:
C
C$$$
      IMPLICIT REAL (A-H,O-Z)
        save
      INCLUDE "parmanl"
      DIMENSION DA( IMAX , KMAX ,6),OBUF( IMAX , LMAX ,6)
      DIMENSION PIN( KMAX ),TIN( KMAX ),SI( KMAXP1 ),SI2( LMAXP1 )
      DIMENSION POUT( LMAXP1 ),OUTZ( LMAXP1 ),TOUT( LMAXP1 )
      DIMENSION MAXI(16),MAXJ(16)
      DIMENSION STATZ(8,16)
      DIMENSION IPMAND(16)
      DATA STATZ/128*0.0/
      DATA MAXI/16*0/,MAXJ/16*0/
      DATA IPMAND/1000,925,850,700,600,500,400,300,
     +             250,200,150,100, 70, 50, 30, 20/
      R = 287.05
      G = 9.8
      ROG = R / G
      DO 15 I=1, IMAX
      PST1 = DA(I,KLEV-1,5)
      DO 7 KK=1,KLEV
C-------
C     KM = MAX0(1,KK-1)
C     KL = KM + 1
C------- KM & KL WERE USED BELOW WHEN IT WAS NECESSARY
C------- TO INTERPOLATE THE MID-POINT TEMPERATURES TO
C------- THE INTERFACES WHERE WE HAVE THE GEOPOTENTIAL
C------- WITH BOTH T & Z AT THE SAME LEVEL, THE QUADRATIC
C------- INTERPOLATION REQUIRES ONLY THE TWO BRACKETTING
C------- LEVELS. INTERFACE TEMPERATURES ARE NOW PASSEDTHRU DA.
C-------
      PIN(KK) = PST1 * SI(KK)
      TIN(KK) = DA(I,KK,2)
C------- THE AFOREMENTIONED  INTERPOLATION OF TEMPERATURE
C     TIN(KK) = DA(I,KM,2) + (DA(I,KL,2) - DA(I,KM,2)) *
C    + LOG(PIN(KK)/DA(I,KM,6))/LOG(DA(I,KL,6)/DA(I,KM,6))
C-------
    7 CONTINUE
      POUT(1) = OBUF(I,LMAX-1,5)
      OUTZ(1) = OBUF(I,LMAX  ,5)
      TOUT(1) = OBUF(I,LMAX-2,5)
      DO 8 KK=2,LMAXP1
      POUT(KK) = OBUF(I,KK-1,6)
      OUTZ(KK) = OBUF(I,KK-1,1)
      TOUT(KK) = OBUF(I,KK-1,2)
    8 CONTINUE
      DO 14 K=1,16
      P = IPMAND(K)
      PLNP = LOG(P)
      DO 1 KK=1,KLEV
      KL = KK
      IF(PIN(KK).LE.P) GO TO 2
    1 CONTINUE
      GO TO 14
    2 CONTINUE
      IF( KL.LE.1 )  GO TO 14
      KM = KL - 1
      ALNPKL = LOG(PIN(KL))
      ALNPKM = LOG(PIN(KM))
      DZ1= DA(I,KM,1)
      DZ = DA(I,KL,1)
      ALNP3 = 0.5 * (ALNPKM + ALNPKL)
      DLNP = ALNPKL - ALNPKM
      RLNP = 1.0 / DLNP
      AZ = (DZ - DZ1) * RLNP
      BT = ROG * (TIN(KL) - TIN(KM)) * RLNP
      Z3 = 0.5*(DZ+DZ1) + 0.125*BT*DLNP*DLNP
      ALNP = (PLNP - ALNP3)
      ZIN = Z3 + (AZ - 0.5*BT*ALNP)*ALNP
      DO 4 KK=1,LMAXP1
      KL = KK
      IF(POUT(KK).LE.P) GO TO 5
    4 CONTINUE
      GO TO 14
    5 CONTINUE
      IF( KL.LE.1 )  GO TO 14
      KM = KL - 1
      ALNPKL = LOG(POUT(KL))
      ALNPKM = LOG(POUT(KM))
      DZ1= OUTZ(KM)
      DZ = OUTZ(KL)
      ALNP3 = 0.5 * (ALNPKM + ALNPKL)
      DLNP = ALNPKL - ALNPKM
      RLNP = 1.0 / DLNP
      AZ = (DZ - DZ1) * RLNP
      BT = ROG * (TOUT(KL) - TOUT(KM)) * RLNP
      Z3 = 0.5*(DZ+DZ1) + 0.125*BT*DLNP*DLNP
      ALNP = (PLNP - ALNP3)
      ZOUT= Z3 + (AZ - 0.5*BT*ALNP)*ALNP
      ZDIF = ZOUT - ZIN
      STATZ(8,K) = STATZ(8,K) + 1.0
      STATZ(1,K) = STATZ(1,K) + ZIN
      STATZ(2,K) = STATZ(2,K) + ZOUT
      STATZ(3,K) = STATZ(3,K) + ZDIF
      STATZ(4,K) = STATZ(4,K) + ZIN*ZIN
      STATZ(5,K) = STATZ(5,K) + ZOUT*ZOUT
      STATZ(6,K) = STATZ(6,K) + ZDIF*ZDIF
      IF(ABS(ZDIF).LE.STATZ(7,K)) GO TO 10
      STATZ(7,K) = ABS(ZDIF)
      MAXI(K) = I
      MAXJ(K) = J
   10 CONTINUE
   14 CONTINUE
   15 CONTINUE
      RETURN
         ENTRY DSTATS
      WRITE(6,20)
   20 FORMAT('1       STATISTICAL ANALYSIS BY PRESSURE ',//,
     + '  N        KP',10X,'BIAS',27X,'RMS',15X,'MAXDIF',//)
      DO 25 K=1,16
      COUNT = MAX(1.0 E+0,STATZ(8,K))
      COUNT = 1.0 / COUNT
      DO 25 I=1,3
      STATZ(I,K) = STATZ(I,K) * COUNT
      STATZ(I+3,K) = SQRT(STATZ(I+3,K)*COUNT)
   25 CONTINUE
      DO 35 KK=1,16
      K=16-KK+1
      KNT = MAX0(1,INT(STATZ(8,K)))
      WRITE(6,30)KNT,IPMAND(K),(STATZ(I,K),I=1,7),MAXI(K),MAXJ(K)
   30 FORMAT(I7,I6,6F10.3,F10.3,2I7)
   35 CONTINUE
      RETURN
      END
