CMAX: [Version 1.0 (6/22/93)]
CMAX: ----------------------------------------------------
CMAX: Transformation of ZUVANL from zuvanl.f
CMAX: ----------------------------------------------------
      SUBROUTINE zuvanl(fcerr,z,u,v,t,psfc,z1000,cov,aquan,asdae,ze,
     &    zelat,zelon,nze,data,nobs)
      
      INTEGER MaxObservPoints
      PARAMETER (MaxObservPoints = 150000)
C     DRIVER FOR THE UPPER AIR  Z - U - V  ANALYSIS
c
c     Maximal observation points
c
c
c
      COMMON / BADLIST / NUMBADS, BADLIST (5000), BADLEV (5000)
      COMMON / scan / maxobs, rscan, onestage, MaxData
      COMMON / pmh / pmand (18)
C     ANALYSIS CONTROL PARAMETERS
      LOGICAL STRAT, SLP, FDCOS, FSOAR, FLOW, SIGDAT
      COMMON / CONTRL / IM, JNP, NB, NDALT, NHMS, NLAY, NLAYM1, NLAYP1, 
     &NYMD, STRAT, IDIMEN, SLP, FDCOS, FSOAR, FLOW, SIGDAT, MAXREP
      COMMON / DATAPTRS / IKX, ILAT, ILON, IPRES, ITYPE, IDELT, IFCER
      COMMON / MLAY / MLEV, MLEVM1, MMIX, MMIXM1, LS, LSM1
      LOGICAL MR
      COMMON / MRFLAG / MR
      COMMON / sparam / lsolve, ns, nr, nv, nh (13104), nq, kq (3), volc
     &entre_x (13104), ngv, volcentre_y (13104), glatup (13104), glonup 
     &(13104)
      COMMON / CONST / SIGE (20), DLAT, DLON, PI, PI180, PTOP, PINT, GI,
     & OMG2RE, RADEI, RADEI2, GRADEI, GRADEI2, G2, POLEFAC, RTD, POLELAT
!hpf$ distribute aquan(block,block,block)
!hpf$ distribute asdae(block,block)
!hpf$ distribute badlev(block)
!hpf$ distribute badlist(block)
!hpf$ distribute data(block)
!hpf$ distribute fcerr(block,block)
!hpf$ distribute glatup(block)
!hpf$ distribute glonup(block)
!hpf$ distribute kq(block)
!hpf$ distribute nh(block)
!hpf$ distribute pmand(block)
!hpf$ distribute psfc(block)
!hpf$ distribute sige(block)
!hpf$ distribute t(block,block)
!hpf$ distribute u(block,block)
!hpf$ distribute ucount(block)
!hpf$ distribute uinc_mean(block)
!hpf$ distribute uinc_stdev(block)
!hpf$ distribute uquan_mean(block)
!hpf$ distribute uquan_stdev(block)
!hpf$ distribute v(block,block)
!hpf$ distribute vcount(block)
!hpf$ distribute vinc_mean(block)
!hpf$ distribute vinc_stdev(block)
!hpf$ distribute volcentre_x(block)
!hpf$ distribute volcentre_y(block)
!hpf$ distribute vquan_mean(block)
!hpf$ distribute vquan_stdev(block)
!hpf$ distribute z(block,block)
!hpf$ distribute z1000(block)
!hpf$ distribute zcount(block)
!hpf$ distribute ze(block)
!hpf$ distribute zelat(block)
!hpf$ distribute zelon(block)
!hpf$ distribute zinc_mean(block)
!hpf$ distribute zinc_stdev(block)
!hpf$ distribute zquan_mean(block)
!hpf$ distribute zquan_stdev(block)
      DATA fill / 1.e10 /
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      INTEGER NZE
      LOGICAL first
      DATA first / .true. /
      SAVE first
      DIMENSION AQUAN (IM * JNP, MLEV + 1, nq)
      DIMENSION ASDAE (IM * JNP, MLEV + 1)
      REAL data (MaxObservPoints * 10)
      EXTERNAL COV
      DIMENSION FCERR (IM * JNP, MLEV)
      INTEGER I
      INTEGER IMJM
      INTEGER LEV
      INTEGER NNQ
      REAL fill
      LOGICAL onestage
      REAL psfc (IM * JNP)
      DIMENSION T (IM * JNP, MLEV)
      DIMENSION U (IM * JNP, MLEV)
      DIMENSION V (IM * JNP, MLEV)
      DIMENSION Z (IM * JNP, MLEV)
      DIMENSION ZE (20000)
      DIMENSION ZELAT (20000)
      DIMENSION ZELON (20000)
      DIMENSION Z1000 (IM * JNP)
      CHARACTER * 60 message

C Statistical arrays.

      real zquan_mean( mlev ), uquan_mean( mlev ), vquan_mean( mlev ),
     &   zinc_mean( mlev ), uinc_mean( mlev ), vinc_mean( mlev ),
     &   zquan_stdev( mlev ), uquan_stdev( mlev ), vquan_stdev( mlev ),
     &   zinc_stdev( mlev ), uinc_stdev( mlev ), vinc_stdev( mlev )
      integer zcount( mlev ), ucount( mlev ), vcount( mlev )

C     INITIALIZE
      INTEGER Ti1
      imjm = im * jnp
      mr = .FALSE.
      WRITE( 6,900 ) 
      message(1:60) = ' '

      zcount( 1:mlev ) = 0
      ucount( 1:mlev ) = 0
      vcount( 1:mlev ) = 0
      zquan_stdev( 1:mlev ) = 0.
      uquan_stdev( 1:mlev ) = 0.
      vquan_stdev( 1:mlev ) = 0.
      zinc_stdev( 1:mlev ) = 0.
      uinc_stdev( 1:mlev ) = 0.
      vinc_stdev( 1:mlev ) = 0.
      zquan_mean( 1:mlev ) = 0.
      uquan_mean( 1:mlev ) = 0.
      vquan_mean( 1:mlev ) = 0.
      zinc_mean( 1:mlev ) = 0.
      uinc_mean( 1:mlev ) = 0.
      vinc_mean( 1:mlev ) = 0.

      IF (nobs .EQ. 0) THEN
        message = ' ZUVANL: No observations present. '
        CALL flash(3,message,1)
        RETURN
      ENDIF
      IF (first) THEN
        CALL obscorrs(pmand,mlev)
        first = .FALSE.
      ENDIF
C     INTERPOLATE FORECAST ERROR TO DATA POINTS
      ifcer = idelt + 3 * nobs
      CALL vterpv_v1(data,ilat,maxobservpoints * 10,data,ilon,
     &    maxobservpoints * 10,data,ipres,maxobservpoints * 10,data,
     &    ifcer,maxobservpoints * 10,nobs,fcerr)
C     SET INDICES FOR ZUV ANALYSIS
      kq(1) = 2
      kq(2) = 3
      kq(3) = 4
C     INITIALIZE THE WORK SPACE
      DO 7006 lev = 1 , mlev
        DO 7000 i = 1 , imjm
          asdae(i,lev) = fcerr(i,lev)
  100     CONTINUE
 7000   CONTINUE
        DO 7004 nnq = 1 , nq
          DO 7002 i = 1 , imjm
            aquan(i,lev,nnq) = 0.0
  105       CONTINUE
 7002     CONTINUE
          CONTINUE
 7004   CONTINUE
        CONTINUE
 7006 CONTINUE
C     FORM AND SOLVE DATA-GRIDPOINT COVARIANCE MATRIXES
c     CALL ZEITBEG('Z SOLVE ')
      CALL solvedrv(1,mlev,fcerr,cov,aquan,asdae,data,nobs)
c     CALL ZEITEND
      IF (numbads .NE. 0) THEN
        WRITE( 6,901 ) numbads
        numbads = 0
      ENDIF
C     ADD ANALYSIS INCREMENTS AND ANALYSIS ERRORS

      DO 200 LEV = 1, MLEV
         DO 190 I = 1, IMJM
            if ( z(i,lev) .ne. fill ) then
              z(i,lev) = aquan(i,lev,1) + z(i,lev)
              zcount(lev) = zcount(lev)+ 1
              zquan_mean(lev) = zquan_mean( lev ) + z(i,lev)
              zinc_mean(lev) = zinc_mean(lev) + aquan(i,lev,1)
            end if
            if ( u(i,lev) .ne. fill ) then
              u(i,lev) = aquan(i,lev,2) + u(i,lev)
              ucount(lev) = ucount(lev) + 1
              uquan_mean(lev) = uquan_mean( lev ) + u(i,lev)
              uinc_mean(lev) = uinc_mean(lev) + aquan(i,lev,2)
            end if
            if ( v(i,lev) .ne. fill ) then
              v(i,lev) = aquan(i,lev,3) + v(i,lev)
              vcount(lev) = vcount(lev) + 1
              vquan_mean(lev) = vquan_mean( lev ) + v(i,lev)
              vinc_mean(lev) = vinc_mean(lev) + aquan(i,lev,3)
            end if
            FCERR(I,LEV) = ASDAE(I,LEV)
  190    CONTINUE
  200 CONTINUE
 
      do k = 1,mlev
       if ( zcount(k) .ne. 0 ) then
         zquan_mean( k ) = zquan_mean( k ) / float( zcount(k) )
         zinc_mean( k )  = zinc_mean( k ) / float( zcount(k) )
       end if
       if ( ucount(k) .ne. 0 ) then
         uquan_mean( k ) = uquan_mean( k ) / float( ucount(k) )
         uinc_mean( k )  = uinc_mean( k ) / float( ucount(k) )
       end if
       if ( vcount(k) .ne. 0 ) then
         vquan_mean( k ) = vquan_mean( k ) / float( vcount(k) )
         vinc_mean( k ) = vinc_mean( k ) / float( vcount(k) )
       end if
      end do
         
      do k = 1,mlev
          do i = 1,imjm
            zquan_stdev(k) = zquan_stdev(k) + ( z(i,k) -
     &                       zquan_mean(k)  )**2
            uquan_stdev(k) = uquan_stdev(k) + ( u(i,k) -
     &                       uquan_mean(k)  )**2
            vquan_stdev(k) = vquan_stdev(k) + ( v(i,k) -
     &                       vquan_mean(k)  )**2
            zinc_stdev(k) = zinc_stdev(k) + ( aquan(i,k,1) -
     &                       zinc_mean(k)  )**2
            uinc_stdev(k) = uinc_stdev(k) + ( aquan(i,k,2) -
     &                       uinc_mean(k)  )**2
            vinc_stdev(k) = vinc_stdev(k) + ( aquan(i,k,3) -
     &                       vinc_mean(k)  )**2
          end do
          if(zcount(k).ne. 0 ) zquan_stdev(k) = sqrt( zquan_stdev(k) /
     &                     float( zcount(k) ) )
          if (ucount(k).ne.0 ) uquan_stdev(k) = sqrt( uquan_stdev(k) /
     &                     float( ucount(k) ) )
          if (vcount(k).ne.0 ) vquan_stdev(k) = sqrt( vquan_stdev(k) /
     &                     float( vcount(k) ) )
          if(zcount(k).ne. 0 ) zinc_stdev(k) = sqrt( zinc_stdev(k) /
     &                     float( zcount(k) ) )
          if (ucount(k).ne.0 ) uinc_stdev(k) = sqrt( uinc_stdev(k) /
     &                     float( ucount(k) ) )
          if (vcount(k).ne.0 ) vinc_stdev(k) = sqrt( vinc_stdev(k) /
     &                     float( vcount(k) ) )
      end do

C Report the results.
      write(6,'(//," HEIGHT ANALYSIS REPORT ",2(2x,i8),//,14x,"Incr",
     &     "ements",12x,"Quantities",/,2x,"Level",3x,"Mean",7x,"Std",
     &     " Dev",3x,"Mean",9x,"Std Dev",t59,"Comments",/,2x,"-----",
     &     3x,18("-"),3x,20("-"),t59,"--------",/)') nymd, nhms

      do i = 1,mlev

        write(6,'(2x,f5.0,3x,f8.2,2x,f8.2,3x,f9.2,2x,f9.2)')
     &        pmand(i), zinc_mean(i), zinc_stdev(i), zquan_mean(i),
     &        zquan_stdev(i)

      end do

      write(6,'(//," U-WIND ANALYSIS REPORT ",2(2x,i8),//,14x,"Incr",
     &     "ements",12x,"Quantities",/,2x,"Level",3x,"Mean",7x,"Std",
     &     " Dev",3x,"Mean",9x,"Std Dev",t59,"Comments",/,2x,"-----",
     &     3x,18("-"),3x,20("-"),t59,"--------",/)') nymd, nhms

      do i = 1,mlev


        write(6,'(2x,f5.0,3x,f8.2,2x,f8.2,3x,f9.2,2x,f9.2)')
     &        pmand(i), uinc_mean(i), uinc_stdev(i), uquan_mean(i),
     &        uquan_stdev(i)

      end do

      write(6,'(//," V-WIND ANALYSIS REPORT ",2(2x,i8),//,14x,"Incr",
     &     "ements",12x,"Quantities",/,2x,"Level",3x,"Mean",7x,"Std",
     &     " Dev",3x,"Mean",9x,"Std Dev",t59,"Comments",/,2x,"-----",
     &     3x,18("-"),3x,20("-"),t59,"--------",/)') nymd, nhms

      do i = 1,mlev

        write(6,'(2x,f5.0,3x,f8.2,2x,f8.2,3x,f9.2,2x,f9.2)')
     &        pmand(i), vinc_mean(i), vinc_stdev(i), vquan_mean(i),
     &        vquan_stdev(i)

      end do
      RETURN
  900 FORMAT(' ZUVANL:  *** Starting upper air H-U-V analysis ***')
  901 FORMAT('0ZUVANL:  *** NO RESOLVING DONE -  NUMBADS = ',I5)
      END


