      PROGRAM ARPSPRT
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                   PROGRAM ARPSPRT                    ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    university of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
c
c#######################################################################
c
c     PURPOSE:
c
c     Program to examine ARPS grids that have been saved as
c     history files.  This gives print-outs of the variables to a
c     file or screen.
c
c     It shares with the model the include files 'dims.inc' for 
c     definition of dimensions and domain size  and  'globcst.inc'
c     for storage parameters.
c
c     It reads in a history file produced by ARPS 3.0 in a user 
c     specified format.
c
c     Parameters grdin,basin,mstin,icein,trbin are read in from the
c     data file itself, therefore are determined internally.
c     Arrays that are not read in retains their initial zero values.
c     These parameters are passed among subroutines through
c     a common block defined in 'indtflg.inc'. 
c
c#######################################################################
c
c     AUTHOR: Keith Brewster  OU School of Meteorology. April 1992
c
c     MODIFICATION HISTORY:
c      14 May 1992  (KB) changed from arps2.5 to arps3.0
c      03 Aug 1992  (KB) updated to account for changes in arps3.0
c
c      8/27/1992 (M. Xue)
c      To call dtaread to read new data format.
c
c      8/30/1992 (K. Brewster)
c      Moved label 101 so that grid file name is not reentered for
c      second and subsequent data reads.
c
c#######################################################################
c
c     DATA ARRAYS READ IN:
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (km)
c       zp       z coordinate of grid points in computational space (m)
c
c       uprt     x component of perturbation velocity (m/s)
c       vprt     y component of perturbation velocity (m/s)
c       wprt     vertical component of perturbation velocity in Cartesian 
c                coordinates (m/s).
c
c       ptprt    perturbation potential temperature (K)
c       pprt     perturbation pressure (Pascal)
c
c       qvprt    perturbation water vapor mixing ratio (kg/kg)
c       qc       Cloud water mixing ratio (kg/kg)
c       qr       Rainwater mixing ratio (kg/kg)
c       qi       Cloud ice mixing ratio (kg/kg)
c       qs       Snow mixing ratio (kg/kg)
c       qh       Hail mixing ratio (kg/kg)
c
c       ubar     Base state x velocity component (m/s)
c       vbar     Base state y velocity component (m/s)
c       wbar     Base state z velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state air density (kg/m**3)
c       qvbar    Base state water vapor mixing ratio (kg/kg)
c
c     CALCULATED DATA ARRAYS:
c
c       u        x component of velocity (m/s)
c       v        y component of velocity (m/s)
c       w        z component of velocity (m/s)
c       pt       potential temperature (K)
c       qv       water vapor mixing ratio (kg/kg)
c    
c
c     WORK ARRAYS:
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c
c     (These arrays are defined and used locally (i.e. inside this
c      subroutine), they may also be passed into routines called by
c      this one. Exiting the call to this subroutine, these temporary
c      work arrays may be used for other purposes, and therefore their
c      contents may be overwritten. Please examine the usage of work 
c      arrays before you make any change to the code.)
c
c#######################################################################
c
c     Variable Declarations:
c
c#######################################################################
c
      implicit none
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'dims.inc'
      include 'globcst.inc'
      include 'indtflg.inc'
c
c#######################################################################
c
c     Arrays to be read in:
c
c#######################################################################
c
      real x     (nx,ny,nz)
      real y     (nx,ny,nz)
      real z     (nx,ny,nz)
      real zp    (nx,ny,nz)
      real j1    (nx,ny,nz)
      real j2    (nx,ny,nz)
      real j3    (nx,ny,nz)
 
      real u     (nx,ny,nz) 
      real v     (nx,ny,nz)
      real w     (nx,ny,nz)
      real pt    (nx,ny,nz)

      real qv    (nx,ny,nz)
      real qc    (nx,ny,nz)
      real qr    (nx,ny,nz)
      real qi    (nx,ny,nz)
      real qs    (nx,ny,nz)
      real qh    (nx,ny,nz)
      real km    (nx,ny,nz)

      real ubar  (nx,ny,nz)
      real vbar  (nx,ny,nz)
      real wbar  (nx,ny,nz)
      real ptbar (nx,ny,nz)
      real rhobar(nx,ny,nz)
      real pbar  (nx,ny,nz)
      real qvbar (nx,ny,nz)
C
      real uprt  (nx,ny,nz)
      real vprt  (nx,ny,nz)
      real wprt  (nx,ny,nz)
      real ptprt (nx,ny,nz)
      real pprt  (nx,ny,nz)
      real rhoprt(nx,ny,nz)
      real qvprt (nx,ny,nz)
c
c#######################################################################
c
c     Other data variables
c
c#######################################################################
c
      real time
c
c#######################################################################
c
c     Temporary working arrays:
c
c#######################################################################
c
      real tem1(nx,ny,nz)
      real tem2(nx,ny,nz)
c
c#######################################################################
c
c     User request stuff
c
c#######################################################################
c
      character*80 inline
      integer maxpar
      parameter (maxpar=10)
      character*2 params(maxpar)
      integer islice,jslice,kslice
      integer ireply,nparams,ip
c
c#######################################################################
c
c     Misc internal variables
c
c#######################################################################
c
      integer i,j,k,ireturn,mode
c
      character grdbasfn*80,filename*80
      character title*25 
      integer ngchan,nchanl,lenfil,lengbf
c
      integer iplot(3)
      integer ib(3),ie(3),jb(3),je(3),kb(3),ke(3)
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
c#######################################################################
c
c     A few initializations
c
c#######################################################################

      DO 10 mode=1,3
        ib(mode)=2
        ie(mode)=nx-2
        jb(mode)=2
        je(mode)=ny-2
        kb(mode)=2
        ke(mode)=nz-2
 10   CONTINUE

c#######################################################################
c
c     Fill all arrays with zero, which is the default value of the
c     data arrays.
c
c#######################################################################
c
      CALL flzero(x    ,nxyz)
      CALL flzero(y    ,nxyz)
      CALL flzero(z    ,nxyz)
      CALL flzero(zp   ,nxyz)
      CALL flzero(j1   ,nxyz)
      CALL flzero(j2   ,nxyz)
      CALL flzero(j3   ,nxyz)

      CALL flzero(u    ,nxyz)
      CALL flzero(v    ,nxyz)
      CALL flzero(w    ,nxyz)
      CALL flzero(pt   ,nxyz)
      CALL flzero(qv   ,nxyz)
      CALL flzero(qc   ,nxyz)
      CALL flzero(qr   ,nxyz)
      CALL flzero(qi   ,nxyz)
      CALL flzero(qs   ,nxyz)
      CALL flzero(qh   ,nxyz)
      CALL flzero(km   ,nxyz)

      CALL flzero(uprt ,nxyz)
      CALL flzero(vprt ,nxyz)
      CALL flzero(wprt ,nxyz)
      CALL flzero(ptprt,nxyz)
      CALL flzero(pprt ,nxyz)
      CALL flzero(qvprt,nxyz)

      CALL flzero(ubar ,nxyz)
      CALL flzero(vbar ,nxyz)
      CALL flzero(wbar ,nxyz)
      CALL flzero(ptbar,nxyz)
      CALL flzero(pbar ,nxyz)
      CALL flzero(rhobar,nxyz)
      CALL flzero(qvbar,nxyz)

      CALL flzero(tem1 ,nxyz)
      CALL flzero(tem2 ,nxyz)


      write(6,'(/11(/5x,a)/)')

     :'###############################################################',
     :'###############################################################',
     :'###                                                         ###',
     :'###                  Welcome to ARPSPRT                     ###',
     :'###      This program reads in the history dump data        ###',
     :'###      sets generated by ARPS, and prints the array       ###',
     :'###      contents as 2-D arrays tables at user specified    ###',
     :'###      slices.                                            ###',
     :'###                                                         ###',
     :'###############################################################',
     :'###############################################################'


      write(6,'(a)')
     :     ' Please set the input data format flag value 1/2/3/4 ' 
      read(5,*) hdmpfmt 
c
c#######################################################################
c
c     Get the name of the input data set.
c
c#######################################################################
c

      write(6,'(/a/a/a)')
     :' Please give the name of the file containing',
     :'             the grid and base state information.',
     :' This is not used if these fields are stored in each data file.'
      read(5,*) grdbasfn

      lengbf=len(grdbasfn)
      CALL strlnth( grdbasfn, lengbf)

101   CONTINUE

      write(6,'(/a)')' Please give the name of the input data file.'
      read(5,*) filename

      lenfil=len(filename)
      CALL strlnth( filename, lenfil)
      write(6,'(/a,a)')' The data set name is ', filename(1:lenfil) 

c
c#######################################################################
c
c     Read all input data arrays
c
c#######################################################################
c
      CALL dtaread(nx,ny,nz, hdmpfmt,grdbasfn(1:lengbf),lengbf,
     :            filename(1:lenfil),lenfil,time,
     :            x,y,z,zp, uprt ,vprt ,wprt ,ptprt, pprt ,
     :            qvprt, qc, qr, qi, qs, qh, km,
     :            ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,
     :            ireturn, tem1)

      curtim = time
c
c#######################################################################
c 
c     ireturn = 0 for a successful read
c
c#######################################################################
c
      IF( ireturn.eq.0 ) THEN   ! successful read

110     CONTINUE
C
c#######################################################################
c
c     Determine variables to plot
c
c#######################################################################
c
        write(6,'(/a,/a,/a,/a/)')
     :' Do you wish to let the program determine the rest of the',
     :' control parameters? If so, the program will show 2-D' ,
     :' x-y, x-z and y-z slices of most fields through the center',
     :' of model domain. Type in 1 or 0 for yes or no.'
        read(5,*) ireply
C
        IF( ireply.ne. 1) THEN
 
          write(6,'(/a,/a/)')
     :  ' Tell the program what variables to show ,',
     :  ' by using the following two character IDs '
          write(6,801) ' u,v,w wind components:         uu,vv,ww'
          write(6,801) ' u,v,w wind base state:         ub,vb,wb'
          write(6,801) ' u,v,w perturbations            up,vp,wp'
          write(6,801) ' potential temp, pressure:      pt,pr'
          write(6,801) ' base pt, pressure:             tb,pb'
          write(6,801) ' theta, pressure perturbations  tp,pp'
          write(6,801) ' Mixing ratio vapor,rain cloud  qv,qr,qc'
          write(6,801) ' Mixing ratio ice,snow,hail     qi,qs,qh'
          write(6,801) ' qv base state,perturbations    qb,qp'
          write(6,801) ' '

          read(5,840) inline
 840      format(a80)
          CALL strcap(inline,inline,80)
          write(6,'(a,/a80)') '  inline: ',inline
          CALL parsln(inline,params,80,maxpar,nparams)

          write(6,'(a,i4,a,20a)') 
     :     ' read ',nparams,' params: ', (params(ip),ip=1,nparams)
          write(6,801) ' '
          write(6,801) ' Now tell the program which slices you want.'
          write(6,801) ' Do you want x-y to plot the slices? (1 or 0)'
          read(5,*) iplot(1)
          write(6,801) ' Do you want x-z to plot the slices? (1 or 0)'
          read(5,*) iplot(2)
          write(6,801) ' Do you want y-z to plot the slices? (1 or 0)'
          read(5,*) iplot(3)

          write(6,801) ' '
          write(6,801)
     :' Now input the i, j or k index for the slice to be selected.'
          write(6,801) 
     :' if i<0, i=(nx-1)/1+1 is assumed. Same is for j and k.'
          write(6,801) ' '
          write(6,'(a,a,i3,a)') 
     :' Input the index k for the x-y slice to be printed.',
     :' (1=< k =<',nz-1,')'
          read(5,*) kslice
          write(6,'(a,a,i3,a)') 
     :' Input the index j for the x-z slice to be printed.',
     :' (1=< j =<',ny-1,')'
          read(5,*) jslice
          write(6,'(a,a,i3,a)') 
     :' Input the index i for the y-z slice to be printed.',
     :' (1=< i =<',nx-1,')'
          read(5,*) islice

          IF( kslice.le.0) kslice = (nz-2)/2+1
          IF( jslice.le.0) jslice = (ny-2)/2+1
          IF( islice.le.0) islice = (nx-2)/2+1
        ELSE
          nparams=5
          params(1)='uu'
          params(2)='vv'
          params(3)='ww'
          params(4)='TH'
          params(5)='PR'
c
c#######################################################################
c
c     Set the printing parameter for each mode
c
c#######################################################################
c
          iplot(1)=1       ! control for plotting x-y slice
          iplot(2)=1       ! control for plotting x-z slice
          iplot(3)=0       ! control for plotting y-z slice

c          kslice = (nz-2)/2+1
          kslice = 3
          jslice = (ny-2)/2+1
c          islice = (nx-2)/2+1
          islice = 8
        END IF
c
c#######################################################################
c
c     Transfer the islice info into the ib,ie vectors
c     Likewise for j and k
c
c#######################################################################
c
        kb(1)=kslice
        ke(1)=kslice

        jb(2)=jslice
        je(2)=jslice

        ib(3)=islice
        ie(3)=islice
c
c#######################################################################
c
c     Calculate total fields from that for base state and perturbations
c
c#######################################################################
c
        DO 90 k=1,nz
        DO 90 j=1,ny
        DO 90 i=1,nx
          u(i,j,k)=uprt(i,j,k)+ubar(i,j,k)
          v(i,j,k)=vprt(i,j,k)+vbar(i,j,k)
          w(i,j,k)=wprt(i,j,k)+wbar(i,j,k)
          pt(i,j,k)=ptprt(i,j,k)+ptbar(i,j,k)
          qv(i,j,k)=qvprt(i,j,k)+qvbar(i,j,k)
 90     CONTINUE
c
c#######################################################################
c
c     Loop through three plotting modes
c     producing x-y, x-z and y-z slices of 2-d plotting each time.
c     If iplot(mode)=0, plotting is skipped. 
c
c#######################################################################
c
        DO 130 mode = 1, 3
          write(6,'(a,i3)') 'mode =',mode
          IF( iplot(mode).NE.0 ) THEN
            DO 120 ip=1,nparams
              write(6,'(a,i4,a,a2)')
     :        '  ip= ',ip,'  param= ',params(ip)
              IF (params(ip) .EQ. 'UU') THEN
                write(title,'(a6,1x,a)') runname(1:6),'u'
                call wrigar( u,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'VV') THEN
                write(title,'(a6,1x,a)') runname(1:6),'v'
                call wrigar( v,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'WW') THEN
                write(title,'(a6,1x,a)') runname(1:6),'w'
                call wrigar( w,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'PT') THEN
                write(title,'(a6,1x,a)') runname(1:6),'pt'
                call wrigar( pt,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'PP') THEN
                write(title,'(a6,1x,a)') runname(1:6),'pprt ' 
                call wrigar( pprt ,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'QV') THEN
                write(title,'(a6,1x,a)') runname(1:6),'qv'    
                call wrigar( qv,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'QC') THEN
                write(title,'(a6,1x,a)') runname(1:6),'qc'    
                call wrigar( qc,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'QR') THEN
                write(title,'(a6,1x,a)') runname(1:6),'qr'    
                call wrigar( qr,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'QI') THEN
                write(title,'(a6,1x,a)') runname(1:6),'qi'    
                call wrigar( qi,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'QS') THEN
                write(title,'(a6,1x,a)') runname(1:6),'qs'    
                call wrigar( qs,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'QH') THEN
                write(title,'(a6,1x,a)') runname(1:6),'qh'    
                call wrigar( qh,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)
C
C  perturbation quantities
C

              ELSE IF (params(ip) .EQ. 'UP') THEN
                write(title,'(a6,1x,a)') runname(1:6),'uprt '
                call wrigar( uprt ,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'VP') THEN
                write(title,'(a6,1x,a)') runname(1:6),'vprt '
                call wrigar( vprt ,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'WP') THEN
                write(title,'(a6,1x,a)') runname(1:6),'wprt '
                call wrigar( wprt ,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'TP') THEN
                write(title,'(a6,1x,a)') runname(1:6),'ptprt'
                call wrigar( ptprt,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'QP') THEN
                write(title,'(a6,1x,a)') runname(1:6),'qvprt '
                call wrigar( qvprt ,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'UB') THEN
                write(title,'(a6,1x,a)') runname(1:6),'u base' 
                call wrigar( ubar,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'VB') THEN
                write(title,'(a6,1x,a)') runname(1:6),'v base' 
                call wrigar( vbar,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'WB') THEN
                write(title,'(a6,1x,a)') runname(1:6),'w base' 
                call wrigar( wbar,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'TB') THEN
                write(title,'(a6,1x,a)') runname(1:6),'thbase' 
                call wrigar( ptbar,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'PB') THEN
                write(title,'(a6,1x,a)') runname(1:6),'pbar' 
                call wrigar( pbar,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              ELSE IF (params(ip) .EQ. 'QB') THEN
                write(title,'(a6,1x,a)') runname(1:6),'qvbase' 
                call wrigar( qvbar,
     :               1,nx,1,ny,1,nz,ib(mode),ie(mode),
     :               jb(mode),je(mode),kb(mode),ke(mode),
     :               title,0.,mode)

              END IF
 120        CONTINUE   ! parameter loop
          END IF  ! iplot(mode).NE.0
 130    CONTINUE   ! mode loop

        write(6,801)
     :  ' Do you want to look at another field in this',
     :  ' data set? (0 or 1)'
        read(5,*) ireply 

        IF( ireply.eq.1 ) GO TO 110

      END IF ! successful read

      write(6,801)
     :  ' Do you want to process another data set'
      read(5,*) ireply 

      IF( ireply.eq.1 ) GO TO 101

 801  format(a)
      STOP
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE FLZERO                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE FLZERO(a,n)

c
c#######################################################################
c
c     PURPOSE:
c
c     Fill an entire array with zeros.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/1991.
c
c     MODIFICATION HISTORY:
c
c     5/02/92 (M. Xue)
c     Added full documentation.
C
c     5/03/92 (M. Xue)
c     Further documentation.
C
c#######################################################################
c
c     INPUT :
c
c       n        The dimension of 1-D array a.
c
c     OUTPUT:
c
c       a        1-D array to be filled with zeros.
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, n

      real a(n)

c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      DO 100 i=1,n
        a(i)=0.0
100   CONTINUE

      RETURN
      END
