PROGRAM arpsdas,175
!
!##################################################################
!##################################################################
!######                                                      ######
!######                   PROGRAM ARPSDAS                    ######
!######             ARPS Data Analysis System                ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!  PURPOSE:
!
!  This program does three-dimensional data analyses on
!  the ARPS grid.  It is loosely based on the rewritten
!  OLAPS surface analysis.
!
!  The driver establishes sizing of variables
!  and adjustable parameters.
!  It calls a routine to set up the grid, calls the data readers,
!  a routine to analyse the data and another to write it out.
!
!  AUTHOR:
!
!  Keith Brewster, CAPS, March, 1994
!  OLAPS sfc analysis Keith Brewster, CAPS, March, 1994
!
!  MODIFICATION HISTORY:
!  01/17/96  (KB) Named changed to ADAS, dimensions imported from
!            adas.dims and input variables from namelist
!            rather than parameter statements in this code.
!  11/27/96  (KB) Cleaned-up to remove unused arrays.
!  01/25/96  (KB) Added tem variables passed to anxiter to accomodate
!            new version of anxiter to speed-up execution.
!  06/15/97  (KB) Update for version 4.3.0, new arguments to call
!            to INITVAR.
!  07/16/97  (KB) Added nt to calling list of INITVAR and made nt
!            equal to one to save memory.
!  05/26/98  (KB) Modified logic in superadiabatic check. (kb)
!  08/31/98  (KB) Added code for nudging to official version (kb)
!  09/22/98  Jason Levit
!            Updated call in initgrdvar for arps4.4.0, including
!            change of dimension for qcumsrc.
!  12/07/98  Added background-based surface correlations to
!            account for cumulus-parameterization small-scale changes.
!  12/20/98  Changed moisture analysis variable to qv.
!            Error stats are still in RH, converted to qv based on
!            background value.
!  12/30/98  Added processing of qr from retrieval data.
!
!  2000/04/14 (Gene Bassett)
!            Merged the ADAS input file back into the ARPS input file.
!
!  2000/05/09 (Gene Bassett)
!            Converted to F90, creating allocation and adas main
!            subroutines.
!
!  2000/07/28 (Ming Xue)
!            Converted to F90 free format. Did away with adaasamin for 
!            better compilation efficiency. Use ALLOCATABLE instead of
!            POINTER allocation to avoid double memory usage.
!
!  07/23/2001 (K. Brewster)
!            Added dynamic array allocation for the IAU increment arrays.
!            Made corresponding changes to the increment handling 
!            routines to account for the array allocation.
!
!  08/07/2001 (K. Brewster)
!            Added a preliminary calculation of vertical velocity 
!            before the cloud analyis (before cmpclddrv call).
!
!  11/01/2001 (K. Brewster)
!            Renamed radar_ref_3d array to ref_mos_3d since its a mosaic.
!            Removed potentially very large 4-d temporary arrays for 
!            radar by reorganizing logic for radar mosaic creation.
!
!  04/10/2002 (K. Brewster)
!            Added code for soil surface temperature adjustment.
!            Includes allocation and deallocation of arrays temporarily
!            needed for radiation calculations.
!
!  05/16/2002 (K. Brewster)
!            Added BMJ arrays to initgrdvar call to conform with
!            latest updates to that subroutine.  Added deallocate after
!            the call to initgrdvar for some arrays not needed after
!            that call.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INCLUDE 'alloc.inc'

  INCLUDE 'adas.inc'
  INCLUDE 'nudging.inc'
  INCLUDE 'exbc.inc'
  INCLUDE 'radcst.inc'

  REAL :: exbcbuf( 1 ) ! EXBC buffer array (unused)
!
!-----------------------------------------------------------------------
!
!  Arrays defining model grid
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: x     (:)      ! The x-coord. of the physical and
                                      ! computational grid. Defined at u-point.
  REAL, ALLOCATABLE :: y     (:)      ! The y-coord. of the physical and
                                      ! computational grid. Defined at v-point.
  REAL, ALLOCATABLE :: z     (:)      ! The z-coord. of the computational grid.
                                      ! Defined at w-point on the staggered grid.
  REAL, ALLOCATABLE :: zp    (:,:,:)  ! The physical height coordinate defined at
                                      ! w-point of the staggered grid.

  REAL, ALLOCATABLE :: hterain(:,:)   ! The height of the terrain.
  REAL, ALLOCATABLE :: mapfct(:,:,:)  ! Map factors at scalar, u and v points

  REAL, ALLOCATABLE :: j1    (:,:,:)  ! Coordinate transformation Jacobian defined
                                      ! as - d( zp )/d( x ).
  REAL, ALLOCATABLE :: j2    (:,:,:)  ! Coordinate transformation Jacobian defined
                                      ! as - d( zp )/d( y ).
  REAL, ALLOCATABLE :: j3    (:,:,:)  ! Coordinate transformation Jacobian defined
                                      ! as d( zp )/d( z ).
  REAL, ALLOCATABLE :: aj3x  (:,:,:)  ! Coordinate transformation Jacobian defined
                                      ! as d( zp )/d( z ) AVERAGED IN THE X-DIR.
  REAL, ALLOCATABLE :: aj3y  (:,:,:)  ! Coordinate transformation Jacobian defined
                                      ! as d( zp )/d( z ) AVERAGED IN THE Y-DIR.
  REAL, ALLOCATABLE :: aj3z  (:,:,:)  ! Coordinate transformation Jacobian defined
                                      ! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
  REAL, ALLOCATABLE :: j3inv (:,:,:)  ! Inverse of j3
!
!-----------------------------------------------------------------------
!
!  ARPS Time-dependent variables
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: u     (:,:,:)  ! Total u-velocity (m/s)
  REAL, ALLOCATABLE :: v     (:,:,:)  ! Total v-velocity (m/s)
  REAL, ALLOCATABLE :: w     (:,:,:)  ! Total w-velocity (m/s)
  REAL, ALLOCATABLE :: wcont (:,:,:)  ! Contravariant vertical velocity (m/s)
  REAL, ALLOCATABLE :: ptprt (:,:,:)  ! Perturbation potential temperature (K)
  REAL, ALLOCATABLE :: pprt  (:,:,:)  ! Perturbation pressure (Pascal)

  REAL, ALLOCATABLE :: qv    (:,:,:)  ! Water vapor specific humidity (kg/kg)
  REAL, ALLOCATABLE :: qc    (:,:,:)  ! Cloud water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: qr    (:,:,:)  ! Rain water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: qi    (:,:,:)  ! Cloud ice mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: qs    (:,:,:)  ! Snow mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: qh    (:,:,:)  ! Hail mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: tke   (:,:,:)  ! Turbulent Kinetic Energy ((m/s)**2)

  REAL, ALLOCATABLE :: ubar  (:,:,:)  ! Base state u-velocity (m/s)
  REAL, ALLOCATABLE :: vbar  (:,:,:)  ! Base state v-velocity (m/s)
  REAL, ALLOCATABLE :: ptbar (:,:,:)  ! Base state potential temperature (K)
  REAL, ALLOCATABLE :: pbar  (:,:,:)  ! Base state pressure (Pascal).
  REAL, ALLOCATABLE :: rhostr(:,:,:)  ! Base state density rhobar times j3.
  REAL, ALLOCATABLE :: qvbar (:,:,:)  ! Base state water vapor specific
                                      ! humidity(kg/kg)

  REAL, ALLOCATABLE :: udteb (:,:)    ! T-tendency of u at e-boundary (m/s**2)
  REAL, ALLOCATABLE :: udtwb (:,:)    ! T-tendency of u at w-boundary (m/s**2)
  REAL, ALLOCATABLE :: vdtnb (:,:)    ! T-tendency of v at n-boundary (m/s**2)
  REAL, ALLOCATABLE :: vdtsb (:,:)    ! T-tendency of v at s-boundary (m/s**2)

  REAL, ALLOCATABLE :: pdteb (:,:)    ! T-tendency of pprt at e-boundary (PASCAL/s)
  REAL, ALLOCATABLE :: pdtwb (:,:)    ! T-tendency of pprt at w-boundary (PASCAL/s)
  REAL, ALLOCATABLE :: pdtnb (:,:)    ! T-tendency of pprt at n-boundary (PASCAL/s)
  REAL, ALLOCATABLE :: pdtsb (:,:)    ! T-tendency of pprt at s-boundary (PASCAL/s)
  REAL, ALLOCATABLE :: trigs1(:)      ! Array containing pre-computed trig
                                      ! function for use in fft.
  REAL, ALLOCATABLE :: trigs2(:)      ! Array containing pre-computed trig
                                      ! function for use in fft.
  INTEGER, ALLOCATABLE :: ifax1(:)    ! Array containing the factors of nx.
  INTEGER, ALLOCATABLE :: ifax2(:)    ! Array containing the factors of ny.
  REAL, ALLOCATABLE :: vwork1 (:,:)   ! 2-D work array for fftopt=2.
  REAL, ALLOCATABLE :: vwork2 (:,:)   ! 2-D work array for fftopt=2.
  REAL, ALLOCATABLE :: wsave1 (:)     ! Work array for fftopt=2.
  REAL, ALLOCATABLE :: wsave2 (:)     ! Work array for fftopt=2.

  REAL, ALLOCATABLE :: ppi(:,:,:)     ! Exner function.
  REAL, ALLOCATABLE :: csndsq(:,:,:)  ! Speed of sound squared
 
  REAL, ALLOCATABLE :: radfrc(:,:,:)  ! Radiation forcing (K/s)
  REAL, ALLOCATABLE :: radsw(:,:)     ! Solar radiation reaching the surface
  REAL, ALLOCATABLE :: rnflx(:,:)     ! Net absorbed radiation by the surface
!
!-----------------------------------------------------------------------
!
!  ARPS Surface variables:
!
!-----------------------------------------------------------------------
!
  INTEGER, ALLOCATABLE :: soiltyp (:,:,:) ! Soil type
  REAL, ALLOCATABLE ::    stypfrct(:,:,:) ! Soil type fraction
  INTEGER, ALLOCATABLE :: vegtyp (:,:)    ! Vegetation type
  REAL, ALLOCATABLE ::    lai    (:,:)    ! Leaf Area Index
  REAL, ALLOCATABLE ::    roufns (:,:)    ! Surface roughness
  REAL, ALLOCATABLE ::    veg    (:,:)    ! Vegetation fraction

  REAL, ALLOCATABLE :: tsfc(:,:,:)    ! Temperature at ground (K)
                                      ! (in top 1 cm layer)
  REAL, ALLOCATABLE :: qvsfc(:,:,:)   ! Effective qv at sfc.
  REAL, ALLOCATABLE :: tsoil(:,:,:)   ! Deep soil temperature (K)
                                      ! (in deep 1 m layer)
  REAL, ALLOCATABLE :: wetsfc(:,:,:)  ! Surface soil moisture in the top
                                      ! 1 cm layer
  REAL, ALLOCATABLE :: wetdp(:,:,:)   ! Deep soil moisture in the deep
                                      ! 1 m layer
  REAL, ALLOCATABLE :: wetcanp(:,:,:) ! Canopy water amount
  REAL, ALLOCATABLE :: snowdpth(:,:)  ! Snow depth (:)

  REAL, ALLOCATABLE :: ptcumsrc(:,:,:)! Source term in pt-equation due
                                      ! to cumulus parameterization
  REAL, ALLOCATABLE :: qcumsrc(:,:,:,:) ! Source term in Qv equation due
                                      ! to cumulus parameterization

  REAL, ALLOCATABLE :: w0avg(:,:,:)   ! a closing running average vertical
                                      ! velocity in 10min for K-F scheme
  REAL, ALLOCATABLE :: kfraincv(:,:)  ! K-F convective rainfall (:)
  INTEGER, ALLOCATABLE :: nca(:,:)    ! K-F counter for CAPE release
  REAL, ALLOCATABLE :: cldefi(:,:)    ! BMJ cloud efficiency
  REAL, ALLOCATABLE :: xland(:,:)     ! BMJ land mask
                                      !   (1.0 = land, 2.0 = sea)
  REAL, ALLOCATABLE :: bmjraincv(:,:) ! BMJ convective rainfall (cm)


  REAL, ALLOCATABLE :: raing(:,:)     ! Grid supersaturation rain
  REAL, ALLOCATABLE :: rainc(:,:)     ! Cumulus convective rain
  REAL, ALLOCATABLE :: prcrate(:,:,:) ! precipitation rate (kg/(m**2*s))
                                      ! prcrate(:,:,:) = total precipitation rate
                                      ! prcrate(:,:,:) = grid scale precip. rate
                                      ! prcrate(:,:,:) = cumulus precip. rate
                                      ! prcrate(:,:,:) = microphysics precip. rate

  REAL, ALLOCATABLE :: usflx (:,:)    ! Surface flux of u-momentum (kg/(m*s**2))
  REAL, ALLOCATABLE :: vsflx (:,:)    ! Surface flux of v-momentum (kg/(m*s**2))
  REAL, ALLOCATABLE :: ptsflx(:,:)    ! Surface heat flux (K*kg/(m*s**2))
  REAL, ALLOCATABLE :: qvsflx(:,:)    ! Surface moisture flux (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
!  Arrays for analysis increment updating (a type of nudging) output
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: uincr(:,:,:)
  REAL, ALLOCATABLE :: vincr(:,:,:)
  REAL, ALLOCATABLE :: wincr(:,:,:)
  REAL, ALLOCATABLE :: pincr(:,:,:)
  REAL, ALLOCATABLE :: ptincr(:,:,:)
  REAL, ALLOCATABLE :: qvincr(:,:,:)
  REAL, ALLOCATABLE :: qcincr(:,:,:)
  REAL, ALLOCATABLE :: qrincr(:,:,:)
  REAL, ALLOCATABLE :: qiincr(:,:,:)
  REAL, ALLOCATABLE :: qsincr(:,:,:)
  REAL, ALLOCATABLE :: qhincr(:,:,:)
!
!-----------------------------------------------------------------------
!
!  Analysis scratch space
!
!-----------------------------------------------------------------------
!
  INTEGER, ALLOCATABLE :: knt(:,:)
  REAL, ALLOCATABLE :: rngsqi(:)
  REAL, ALLOCATABLE :: wgtsum(:,:)
  REAL, ALLOCATABLE :: zsum(:,:)
  REAL, ALLOCATABLE :: sqsum(:,:)
!
!-----------------------------------------------------------------------
!
!  Additional grid variables
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: xs(:)
  REAL, ALLOCATABLE :: ys(:)
  REAL, ALLOCATABLE :: zs(:,:,:)
  REAL, ALLOCATABLE :: latgr(:,:)
  REAL, ALLOCATABLE :: longr(:,:)
!
!-----------------------------------------------------------------------
!
!  Temporary arrays
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: tem1  (:,:,:)     ! Temporary work array.
  REAL, ALLOCATABLE :: tem2  (:,:,:)     ! Temporary work array.
  REAL, ALLOCATABLE :: tem3  (:,:,:)     ! Temporary work array.
  REAL, ALLOCATABLE :: tem4  (:,:,:)     ! Temporary work array.
  REAL, ALLOCATABLE :: tem5  (:,:,:)     ! Temporary work array.
  REAL, ALLOCATABLE :: tem6  (:,:,:)     ! Temporary work array.
  REAL, ALLOCATABLE :: tem7  (:,:,:)     ! Temporary work array.
  REAL, ALLOCATABLE :: tem8  (:,:,:)     ! Temporary work array.
  REAL, ALLOCATABLE :: tem9  (:,:,:)     ! Temporary work array.
  REAL, ALLOCATABLE :: tem10 (:,:,:)     ! Temporary work array.
  INTEGER, ALLOCATABLE :: ibegin(:)
  INTEGER, ALLOCATABLE :: iend(:)
!
!-----------------------------------------------------------------------
!
!  Radiation variables only needed for adjtsfc
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: rsirbm(:,:)
  REAL, ALLOCATABLE :: rsirdf(:,:)
  REAL, ALLOCATABLE :: rsuvbm(:,:)
  REAL, ALLOCATABLE :: rsuvdf(:,:)
  REAL, ALLOCATABLE :: cosz(:,:)
  REAL, ALLOCATABLE :: cosss(:,:)
  REAL, ALLOCATABLE :: fdirir(:,:)
  REAL, ALLOCATABLE :: fdifir(:,:)
  REAL, ALLOCATABLE :: fdirpar(:,:)
  REAL, ALLOCATABLE :: fdifpar(:,:)
  REAL, ALLOCATABLE :: st4(:,:)
  REAL, ALLOCATABLE :: radbuf(:)

  REAL, ALLOCATABLE :: tem11 (:,:,:)     ! Temporary work array.
  REAL, ALLOCATABLE :: tem12 (:,:,:)     ! Temporary work array.
  REAL, ALLOCATABLE :: tem13 (:,:,:)     ! Temporary work array.
  REAL, ALLOCATABLE :: tem14 (:,:,:)     ! Temporary work array.
  REAL, ALLOCATABLE :: tem15 (:,:,:)     ! Temporary work array.
  REAL, ALLOCATABLE :: tem16 (:,:,:)     ! Temporary work array.

!
!-----------------------------------------------------------------------
!
!  Cloud analysis variables other than model's ones:
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: ref_mos_3d (:,:,:)
                     ! 3D remapped radar reflectivity field as input
                     ! may be modified by sate. VIS as output
  REAL, ALLOCATABLE :: w_cld (:,:,:)      ! vertical velocity (m/s) in clouds
  INTEGER(kind=selected_int_kind(4)), ALLOCATABLE :: cldpcp_type_3d(:,:,:)
                                       ! cloud/precip type field

  INTEGER, ALLOCATABLE :: icing_index_3d(:,:,:) ! icing severity index
  LOGICAL, ALLOCATABLE :: l_mask_pcptype(:,:) ! analyze precip type 
!
!-----------------------------------------------------------------------
!
!  Analysis variables
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: anx(:,:,:,:)
!
!-----------------------------------------------------------------------
!
!  Cross-category correlations
!
!-----------------------------------------------------------------------
!
  REAL,    ALLOCATABLE :: xcor(:,:)
  INTEGER, ALLOCATABLE :: icatg(:,:)
!
!-----------------------------------------------------------------------
!
!  ARPS dimensions:
!
!  nx, ny, nz: Dimensions of analysis grid.
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx       ! Number of grid points in the x-direction
  INTEGER :: ny       ! Number of grid points in the y-direction
  INTEGER :: nz       ! Number of grid points in the z-direction

  INTEGER :: nxndg    ! Number of x grid points for IAU
  INTEGER :: nyndg    ! Number of y grid points for IAU
  INTEGER :: nzndg    ! Number of z grid points for IAU
!
!-----------------------------------------------------------------------
!
!  Soil types.
!
!-----------------------------------------------------------------------
!
  INTEGER :: nstyps            ! Number of soil types
!
!-----------------------------------------------------------------------
!
!  Misc
!
!-----------------------------------------------------------------------
!
  INTEGER :: nt                ! Number of time levels of data
  INTEGER :: ncat              ! Number of correlation categories

!-----------------------------------------------------------------------
!
!  Analysis parameters
!
!-----------------------------------------------------------------------
!
  REAL :: rhmin
  INTEGER :: iqspr,klim
  PARAMETER (rhmin=0.05,  & ! rh safety net value to prevent neg qv
             iqspr=3,     & ! Use qobs of pstn to combine x,y,elev
             klim= 1)       ! Min of one other sfc station for QC
!
!-----------------------------------------------------------------------
!
!  ARPS cloud cover analysis variables:
!
!-----------------------------------------------------------------------
!
! Remapping parameters of the remapped radar data grid
!
  INTEGER :: strhopt_rad                  ! streching option
  INTEGER :: mapproj_rad                  ! map projection indicator
  REAL :: dx_rad,dy_rad,dz_rad,dzmin_rad  ! grid spcngs
  REAL :: ctrlat_rad,ctrlon_rad           ! central lat and lon
  REAL :: tlat1_rad,tlat2_rad,tlon_rad    ! true lat and lon
  REAL :: scale_rad                       ! map scale factor
!
!-----------------------------------------------------------------------
!
!  Common block that stores remapping parameters for the radar
!  data file.
!
!-----------------------------------------------------------------------
!
  COMMON /remapfactrs_rad/ strhopt_rad,mapproj_rad
  COMMON /remapfactrs_rad2/ dx_rad,dy_rad,dz_rad,dzmin_rad,             &
           ctrlat_rad,ctrlon_rad,                                       &
           tlat1_rad,tlat2_rad,tlon_rad,scale_rad
!
!-----------------------------------------------------------------------
!
!  Non-dimensional smoothing coefficient
!
!-----------------------------------------------------------------------
!
  INTEGER, PARAMETER :: tsfcadj=1
!
!-----------------------------------------------------------------------
!
!  ARPS include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
  INCLUDE 'bndry.inc'
  INCLUDE 'grid.inc'        ! Grid parameters
  INCLUDE 'adjust.inc'
!
!-----------------------------------------------------------------------
!
!  Surface Station variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: isrcsng(mx_sng)
  INTEGER :: icatsng(mx_sng)
  REAL :: latsng(mx_sng,ntime),lonsng(mx_sng,ntime)
  REAL :: hgtsng(mx_sng,ntime)
  REAL :: xsng(mx_sng),ysng(mx_sng)
  INTEGER :: timesng(mx_sng,ntime)
  CHARACTER (LEN=5) :: stnsng(mx_sng,ntime)
!
!-----------------------------------------------------------------------
!
!  Surface (single-level) read-in observation variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=8) :: wx(mx_sng,ntime)
  CHARACTER (LEN=8) :: csrcsng(mx_sng,ntime)
  CHARACTER (LEN=1) :: store_emv(mx_sng,5,ntime)
  CHARACTER (LEN=4) :: store_amt(mx_sng,5,ntime)
  INTEGER :: kloud(mx_sng,ntime),idp3(mx_sng,ntime)
  REAL :: store_hgt(mx_sng,5,ntime)
  REAL :: obrdsng(mx_sng,nvar_sng,ntime)
  REAL :: obsng(nvar_anx,mx_sng)
  REAL :: odifsng(nvar_anx,mx_sng)
  REAL :: corsng(nvar_anx,mx_sng)
  REAL :: oanxsng(nvar_anx,mx_sng)
  REAL :: thesng(mx_sng)
  INTEGER :: ival(mx_sng)
!
!-----------------------------------------------------------------------
!
!  Upper Air Station variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: isrcua(mx_ua)
  REAL :: xua(mx_ua),yua(mx_ua)
  REAL :: elevua(mx_ua)
  REAL :: hgtua(nz_ua,mx_ua)
  INTEGER :: nlevsua(mx_ua)
  CHARACTER (LEN=5) :: stnua(mx_ua)
!
!-----------------------------------------------------------------------
!
!  Upper-air observation variables
!
!-----------------------------------------------------------------------
!
  REAL :: obsua(nvar_anx,nz_ua,mx_ua)
  REAL :: odifua(nvar_anx,nz_ua,mx_ua)
  REAL :: corua(nvar_anx,nz_ua,mx_ua)
  REAL :: oanxua(nvar_anx,nz_ua,mx_ua)
  REAL :: theua(nz_ua,mx_ua)
!
!-----------------------------------------------------------------------
!
!  Radar site variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=5) :: stnrad(mx_rad)
  INTEGER :: isrcrad(0:mx_rad)
  REAL :: latrad(mx_rad),lonrad(mx_rad)
  REAL :: elvrad(mx_rad)
!
!-----------------------------------------------------------------------
!
!  Radar observation variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: irad(mx_colrad)
  INTEGER :: nlevrad(mx_colrad)
  REAL :: latradc(mx_colrad),lonradc(mx_colrad)
  REAL :: xradc(mx_colrad),yradc(mx_colrad)
  REAL :: distrad(mx_colrad),uazmrad(mx_colrad),vazmrad(mx_colrad)
  REAL :: hgtradc(nz_rdr,mx_colrad)
  REAL :: theradc(nz_rdr,mx_colrad)
  REAL :: obsrad(nvar_radin,nz_rdr,mx_colrad)
  REAL :: odifrad(nvar_rad,nz_rdr,mx_colrad)
  REAL :: corrad(nvar_rad,nz_rdr,mx_colrad)
  REAL :: oanxrad(nvar_rad,nz_rdr,mx_colrad)
!
!-----------------------------------------------------------------------
!
!  Retrieval radar variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=5) :: stnret(mx_ret)
  INTEGER :: isrcret(0:mx_ret)
  REAL :: latret(mx_ret),lonret(mx_ret)
  REAL :: elvret(mx_ret)
!
!-----------------------------------------------------------------------
!
!  Retrieval observation variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iret(mx_colret)
  INTEGER :: nlevret(mx_colret)
  REAL :: latretc(mx_colret),lonretc(mx_colret)
  REAL :: xretc(mx_colret),yretc(mx_colret)
!  real distret(mx_colret),uazmret(mx_colret),vazmret(mx_colret)
  REAL :: hgtretc(nz_ret,mx_colret)
  REAL :: theretc(nz_ret,mx_colret)
  REAL :: obsret(nvar_anx,nz_ret,mx_colret)
  REAL :: qrret(nz_ret,mx_colret)
  REAL :: odifret(nvar_anx,nz_ret,mx_colret)
  REAL :: corret(nvar_anx,nz_ret,mx_colret)
  REAL :: oanxret(nvar_anx,nz_ret,mx_colret)
!
!-----------------------------------------------------------------------
!
!  Namen de variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=6) :: var_sfc(nvar_sng)
  DATA var_sfc/ 't     ','td    ','dd    ','ff    ',                    &
                'urot  ','vrot  ','pstn  ','pmsl  ',                    &
                'alt   ','ceil  ','low   ','cover ',                    &
                'solar ','vis'/
  CHARACTER (LEN=6) :: var_ua(nvar_ua)
  DATA var_ua / 'press ',                                               &
                'temp  ','dewpt ','dd    ','ff    '/
  CHARACTER (LEN=6) :: var_anx(nvar_anx)
  DATA var_anx/ 'u     ','v     ',                                      &
                'press ','theta ','qv    '/
!
!-----------------------------------------------------------------------
!
!  Source-dependent parameters
!  Qsrc is the expected square error it is used for
!  setting the QC threshold (qclim) and for relative
!  weighting of observations.
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=80) :: srcback
  REAL :: hqback(nz_tab)
  INTEGER :: nlvqback
  REAL :: qback(nvar_anx,nz_tab)

  CHARACTER (LEN=80) :: srcsng_full(nsrc_sng)
  REAL :: qsrcsng(nvar_anx,nsrc_sng)

  CHARACTER (LEN=80) :: srcua_full(nsrc_ua)
  INTEGER :: nlvuatab(nsrc_ua)
  REAL :: huaqsrc(nz_tab,nsrc_ua)
  REAL :: qsrcua(nvar_anx,nz_tab,nsrc_ua)

  CHARACTER (LEN=80) :: srcrad_full(nsrc_rad)
  REAL :: qsrcrad(nvar_radin,nsrc_rad)

  CHARACTER (LEN=80) :: srcret_full(nsrc_ret)
  INTEGER :: nlvrttab(nsrc_ret)
  REAL :: hrtqsrc(nz_tab,nsrc_ret)
  REAL :: qsrcret(nvar_anx,nz_tab,nsrc_ret)
!
!-----------------------------------------------------------------------
!
!  Quality Control Parameters and Variables
!
!-----------------------------------------------------------------------
!
  REAL :: rmiss
  PARAMETER (rmiss=-99.)
!
  INTEGER :: qualrdsng(mx_sng,nvar_sng,ntime)
  INTEGER :: qualsng(nvar_anx,mx_sng)
  REAL :: qobsng(nvar_anx,mx_sng)
  REAL :: qcmulsng(nsrc_sng)
  REAL :: qcthrsng(nvar_anx,nsrc_sng)
!
  INTEGER :: qualua(nvar_anx,nz_ua,mx_ua)
  REAL :: qobsua(nvar_anx,nz_ua,mx_ua)
  REAL :: qcmulua(nsrc_ua)
  REAL :: qcthrua(nvar_anx,nsrc_ua)
!
  INTEGER :: qualrad(nvar_rad,nz_rdr,mx_colrad)
  REAL :: qobsrad(nvar_rad,nz_rdr,mx_colrad)
  REAL :: qcmulrad(nsrc_rad)
  REAL :: qcthrrad(nvar_radin,nsrc_rad)
!
  INTEGER :: qualret(nvar_anx,nz_ret,mx_colret)
  REAL :: qobsret(nvar_anx,nz_ret,mx_colret)
  REAL :: qcmulret(nsrc_ret)
  REAL :: qcthrret(nvar_anx,nsrc_ret)
!
!-----------------------------------------------------------------------
!
!  Climin is the minimum possible value of each observed variable.
!  Climax is the maximum possible value of each observed variable.
!  Used for surface data only.
!
!-----------------------------------------------------------------------
!
  REAL :: climin(nvar_sng),climax(nvar_sng)
  DATA climin /  -50.,    -50.,    0.,    0.,                           &
                -100.,   -100.,  700.,  880.,                           &
                 880.,      0.,    0.,    0.,                           &
                   0.,      0./
  DATA climax /  125.,    125.,  360.,  100.,                           &
                 100.,    100., 1100., 1090.,                           &
                1090.,  30000.,30000.,    1.,                           &
                1500.,    150./
!
!-----------------------------------------------------------------------
!
!  Filenames and such
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=12) :: suffix
  CHARACTER (LEN=128) :: froot
  CHARACTER (LEN=132) :: qclist
  CHARACTER (LEN=132) :: qcsave
  CHARACTER (LEN=132) :: stats
  INTEGER :: istatus,jstatus
  INTEGER :: nobsng,nobsrd(ntime)
  INTEGER :: i4timanx
  INTEGER :: lenfnm,maxsuf,lensuf,dotloc,mxroot,lenroot
  INTEGER :: iqclist,iqcsave,iwstat

  CHARACTER (LEN=132) :: status_file
  INTEGER n_used_sng(nsrc_sng), n_used_ua(nsrc_sng)
  INTEGER jsta, klev, ngoodlev
!
!-----------------------------------------------------------------------
!
!  Physical constants
!
!-----------------------------------------------------------------------
!
  REAL :: kts2ms,mb2pa,f2crat
  PARAMETER (kts2ms=0.514444,                                           &
             mb2pa=100.,                                                &
             f2crat=(5./9.))
!
!-----------------------------------------------------------------------
!
!  ptlapse is minimum potential temperature lapse rate to use in
!  superadiabatic check.  Used to ensure initial static stability.
!
!-----------------------------------------------------------------------

  REAL :: ptlapse
  PARAMETER (ptlapse = (1.0/4000.) )   ! deg K / m
!
!-----------------------------------------------------------------------
!
!  Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_qvsat

!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=8) :: rdsource
  CHARACTER (LEN=80) :: basdmpfn
  CHARACTER (LEN=80) :: header
  INTEGER :: i,j,k,ios,ktab,isrc,jsrc,ivar,ifile
  INTEGER :: grdbas,lbasdmpf,rbufsz
  INTEGER :: nobsua,ncolrad,ncolret

  REAL :: tgrid,qvmin,qvsat,qsrcmax,rngsq,sprkm
  REAL :: temp,tvbar,qvprt,qtot,pconst,pr1,pr2
!
  INTEGER :: ixrad(mx_rad),jyrad(mx_rad)

  INTEGER :: retqropt
  PARAMETER (retqropt=1)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  nt = 1    ! Number of time levels of data
  ncat = 4  ! Number of correlation categories

!-----------------------------------------------------------------------
!  Set grid mode to non-nesting and grid index, mgrid, to 1.
!-----------------------------------------------------------------------

  nestgrd=0
  mgrid=1
 
!-----------------------------------------------------------------------
!  read in ARPS namelists
!-----------------------------------------------------------------------
 
  CALL initpara(nx,ny,nz,nstyps)

  IF (lbcopt == 2) THEN
    WRITE (*,*) "INITADAS: resetting lbcopt to 1 ",                     &
                "& lateral bc's to 4"
    lbcopt = 1
    ebc = 4
    wbc = 4
    nbc = 4
    sbc = 4
  END IF

!-----------------------------------------------------------------------
!  Read in adas namelist parameters
!-----------------------------------------------------------------------
!
  CALL initadas

  IF(incrdmp > 0 ) THEN
    nxndg=nx
    nyndg=ny
    nzndg=nz
  ELSE
    nxndg=1
    nyndg=1
    nzndg=1
  END IF
!
!-----------------------------------------------------------------------
!
!  Allocate adas arrays
!
!-----------------------------------------------------------------------
!
  ALLOCATE(  CALL check_alloc_status(istatus, "adas:x")
  ALLOCATE(  CALL check_alloc_status(istatus, "adas:y")
  ALLOCATE(  CALL check_alloc_status(istatus, "adas:z")

  ALLOCATE(hterain(nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:hterain")
  ALLOCATE(mapfct (nx,ny,8),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:mapfct")

  ALLOCATE(zp  (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:zp")
  ALLOCATE(j1  (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:j1")
  ALLOCATE(j2  (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:j2")
  ALLOCATE(j3  (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:j3")
  ALLOCATE(aj3x(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:aj3x")
  ALLOCATE(aj3y(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:aj3y")
  ALLOCATE(aj3z(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:aj3z")
  ALLOCATE(j3inv(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:j3inv")

  ALLOCATE(u    (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:u")
  ALLOCATE(v    (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:v")
  ALLOCATE(w    (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:w")
  ALLOCATE(wcont(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:wcont")
  ALLOCATE(ptprt(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:ptprt")
  ALLOCATE(pprt (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:pprt")

  ALLOCATE(qv   (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qv")
  ALLOCATE(qc   (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qc")
  ALLOCATE(qr   (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qr")
  ALLOCATE(qi   (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qi")
  ALLOCATE(qs   (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qs")
  ALLOCATE(qh   (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qh")
  ALLOCATE(tke  (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:tke")

  ALLOCATE(ubar (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:ubar")
  ALLOCATE(vbar (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:vbar")
  ALLOCATE(ptbar(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:ptbar")
  ALLOCATE(pbar (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:pbar")
  ALLOCATE(rhostr(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:rhostr")
  ALLOCATE(qvbar(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qvbar")

  ALLOCATE(udteb(ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:udteb")
  ALLOCATE(udtwb(ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:udtbw")
  ALLOCATE(vdtnb(nx,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:vdtnb")
  ALLOCATE(vdtsb(nx,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:vdtsb")

  ALLOCATE(pdteb(ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:pdteb")
  ALLOCATE(pdtwb(ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:pdtwb")
  ALLOCATE(pdtnb(nx,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:pdtnb")
  ALLOCATE(pdtsb(nx,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:pdtsb")

  ALLOCATE(trigs1(3*(nx-1)/2+1),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:trigs1")
  ALLOCATE(trigs2(3*(ny-1)/2+1),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:trigs2")
  ALLOCATE(ifax1(13),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:ifax1")
  ALLOCATE(ifax2(13),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:ifax2")
   
  ALLOCATE(vwork1(nx+1,ny+1),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:vwork1")
  ALLOCATE(vwork2(ny,nx+1),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:vwork2")

  ALLOCATE(wsave1(3*(ny-1)+15),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:wsave1")
  ALLOCATE(wsave2(3*(nx-1)+15),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:wsave2")

  ALLOCATE(ppi   (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:ppi")
  ALLOCATE(csndsq(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:csndsq")

  ALLOCATE(radfrc(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:radfrc")
  ALLOCATE(radsw (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:radsw")
  ALLOCATE(rnflx (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:rnflx")

  ALLOCATE(soiltyp(nx,ny,nstyps),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:soiltyp")
  ALLOCATE(stypfrct(nx,ny,nstyps),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:stypfrct")
  ALLOCATE(vegtyp (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:vegtyp")
  ALLOCATE(lai    (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:lai")
  ALLOCATE(roufns (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:roufns")
  ALLOCATE(veg    (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:veg")

  ALLOCATE(tsfc   (nx,ny,0:nstyps),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:tsfc")
  ALLOCATE(qvsfc  (nx,ny,0:nstyps),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qvsfc")
  ALLOCATE(tsoil  (nx,ny,0:nstyps),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:tsoil")
  ALLOCATE(wetsfc (nx,ny,0:nstyps),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:wetsfc")
  ALLOCATE(wetdp  (nx,ny,0:nstyps),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:wetdp")
  ALLOCATE(wetcanp(nx,ny,0:nstyps),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:wetcanp")
  ALLOCATE(snowdpth(nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:snowdpth")

  ALLOCATE(ptcumsrc(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:ptcumsrc")
  ALLOCATE(qcumsrc (nx,ny,nz,5),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qcumsrc")
  ALLOCATE(w0avg   (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:w0avg")
  ALLOCATE(nca    (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:nca")
  ALLOCATE(kfraincv (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:kfraincv")
  ALLOCATE(cldefi (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:cldefi")
  ALLOCATE(xland  (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:xland")
  ALLOCATE(bmjraincv (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:bmjraincv")
  ALLOCATE(raing  (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:raing")
  ALLOCATE(rainc  (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:rainc")
  ALLOCATE(prcrate(nx,ny,4),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:prcrate")

  ALLOCATE(usflx (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:usflx")
  ALLOCATE(vsflx (nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:vsflx")
  ALLOCATE(ptsflx(nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:ptsflx")
  ALLOCATE(qvsflx(nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qvsflx")

  ALLOCATE(uincr(nxndg,nyndg,nzndg),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:uincr")
  ALLOCATE(vincr(nxndg,nyndg,nzndg),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:vincr")
  ALLOCATE(wincr(nxndg,nyndg,nzndg),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:wincr")
  ALLOCATE(pincr(nxndg,nyndg,nzndg),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:pincr")
  ALLOCATE(ptincr(nxndg,nyndg,nzndg),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:ptincr")
  ALLOCATE(qvincr(nxndg,nyndg,nzndg),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qvincr")
  ALLOCATE(qcincr(nxndg,nyndg,nzndg),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qcincr")
  ALLOCATE(qrincr(nxndg,nyndg,nzndg),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qrincr")
  ALLOCATE(qiincr(nxndg,nyndg,nzndg),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qiincr")
  ALLOCATE(qsincr(nxndg,nyndg,nzndg),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qsincr")
  ALLOCATE(qhincr(nxndg,nyndg,nzndg),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:qhincr")

  ALLOCATE(knt(nvar_anx,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:knt")
  ALLOCATE(rngsqi(nvar_anx),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:rngsqi")
  ALLOCATE(wgtsum(nvar_anx,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:wgtsum")
  ALLOCATE(zsum(nvar_anx,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:zsum")
  ALLOCATE(sqsum(nvar_anx,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:sqsum")

  ALLOCATE(xs(nx),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:xs")
  ALLOCATE(ys(ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:ys")
  ALLOCATE(zs(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:zs")
  ALLOCATE(latgr(nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:latgr")
  ALLOCATE(longr(nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:longr")

  ALLOCATE(tem1(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:tem1")
  ALLOCATE(tem2(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:tem2")
  ALLOCATE(tem3(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:tem3")
  ALLOCATE(tem4(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:tem4")
  ALLOCATE(tem5(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:tem5")
  ALLOCATE(tem6(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:tem6")
  ALLOCATE(tem7(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:tem7")
  ALLOCATE(tem8(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:tem8")
  ALLOCATE(tem9(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:tem9")
  ALLOCATE(tem10(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:tem10")
  ALLOCATE(ibegin(ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:ibegin")
  ALLOCATE(iend(ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:iend")

  ALLOCATE(anx(nx,ny,nz,nvar_anx),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:anx")

  ALLOCATE(xcor(ncat,ncat),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:xcor")
  ALLOCATE(icatg(nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, "adas:icatg")

  WRITE (*,*) "Finished allocating arrays"
 
!-----------------------------------------------------------------------
! Initialize the allocated arrays to zero.
!-----------------------------------------------------------------------

  x=0.0
  y=0.0
  z=0.0

  hterain=0.0
  mapfct =0.0

  zp  =0.0
  j1  =0.0
  j2  =0.0
  j3  =0.0
  aj3x=0.0
  aj3y=0.0
  aj3z=0.0
  j3inv=0.0

  u    =0.0
  v    =0.0
  w    =0.0
  wcont=0.0
  ptprt=0.0
  pprt =0.0

  qv   =0.0
  qc   =0.0
  qr   =0.0
  qi   =0.0
  qs   =0.0
  qh   =0.0
  tke  =0.0

  ubar =0.0
  vbar =0.0
  ptbar=0.0
  pbar =0.0
  rhostr=0.0
  qvbar=0.0

  udteb=0.0
  udtwb=0.0
  vdtnb=0.0
  vdtsb=0.0

  pdteb=0.0
  pdtwb=0.0
  pdtnb=0.0
  pdtsb=0.0

  trigs1=0.0
  trigs2=0.0
  ifax1=0.0
  ifax2=0.0
   
  vwork1=0.0
  vwork2=0.0

  wsave1=0.0
  wsave2=0.0

  ppi   =0.0
  csndsq=0.0

  radfrc=0.0
  radsw =0.0
  rnflx =0.0

  soiltyp=0.0
  stypfrct=0.0
  vegtyp =0.0
  lai    =0.0
  roufns =0.0
  veg    =0.0

  tsfc   =0.0
  qvsfc  =0.0
  tsoil  =0.0
  wetsfc =0.0
  wetdp  =0.0
  wetcanp=0.0
  snowdpth=0.0

  ptcumsrc=0.0
  qcumsrc =0.0
  w0avg   =0.0
  nca    =0.0
  kfraincv =0.0
  cldefi = 0.0
  xland = 0.0
  bmjraincv = 0.0
  raing  =0.0
  rainc  =0.0
  prcrate=0.0

  usflx =0.0
  vsflx =0.0
  ptsflx=0.0
  qvsflx=0.0

  knt=0.0
  rngsqi=0.0
  wgtsum=0.0
  zsum=0.0
  sqsum=0.0

  xs=0.0
  ys=0.0
  zs=0.0
  latgr=0.0
  longr=0.0

  tem1=0.0
  tem2=0.0
  tem3=0.0
  tem4=0.0
  tem5=0.0
  tem6=0.0
  tem7=0.0
  tem8=0.0
  tem9=0.0
  tem10=0.0
!wdt-williams update
  ibegin=0
  iend=0

  anx=0.0

  xcor=0.0
  icatg=0.0
!
!-----------------------------------------------------------------------
!
!  Set expected _squared_ background error for each variable
!  This will depend on the particular background field used,
!  season, age of forecast, etc.
!
!-----------------------------------------------------------------------
!
  PRINT *, 'Reading ', TRIM(backerrfil)
  OPEN(4,FILE=trim(backerrfil),STATUS='old')
  READ(4,'(a80)') srcback
  READ(4,'(a80)') header
  WRITE(6,'(//,a,/a)') 'Background Std Error for',srcback
  WRITE(6,'(1x,a)')                                                     &
      '   k    hgt(m)  u(m/s)  v(m/s) pres(mb) temp(K)  RH(%)'
  DO ktab=1,nz_tab
    READ(4,*,END=2) hqback(ktab),                                       &
                qback(3,ktab),                                          &
                qback(4,ktab),                                          &
                qback(5,ktab),                                          &
                qback(1,ktab),                                          &
                qback(2,ktab)
    WRITE(6,'(1x,i4,f10.2,5f8.2)') ktab,hqback(ktab),                   &
                               (qback(ivar,ktab),ivar=1,5)
    qback(3,ktab)=100.*qback(3,ktab)
    qback(5,ktab)=0.01*qback(5,ktab)
  END DO
  2 nlvqback=ktab-1
  CLOSE(4)
!
  isrc=1
  DO jsrc=1,nsrc_sng
    IF(srcsng(jsrc) /= 'NULL') THEN
      PRINT *, 'Reading ', TRIM(sngerrfil(jsrc))
      OPEN(4,ERR=3,FILE=trim(sngerrfil(jsrc)),STATUS='old')
      READ(4,'(a8)',ERR=3,END=3) rdsource
      IF(rdsource /= srcsng(jsrc)) THEN
        WRITE(6,'(a,i4,a,a,a,a,a)')                                     &
            ' Mismatch of source names',jsrc,                           &
            ' read ',rdsource,' expected ',srcsng(jsrc)
        STOP
      END IF
      READ(4,'(a80)',ERR=3,END=3) srcsng_full(isrc)
      READ(4,*,ERR=3,END=3) qcmulsng(isrc)
      READ(4,'(a80)',ERR=3,END=3) header
      READ(4,*,ERR=3,END=3)                                             &
             qsrcsng(3,isrc),                                           &
             qsrcsng(4,isrc),                                           &
             qsrcsng(5,isrc),                                           &
             qsrcsng(1,isrc),                                           &
             qsrcsng(2,isrc)
      WRITE(6,'(//,a,a,/a)') 'Single-level std error for ',             &
                            srcsng(isrc),srcsng_full(isrc)
      WRITE(6,'(a,f5.1)') ' QC multiplier: ',qcmulsng(isrc)
      WRITE(6,'(1x,a)')                                                 &
          '   u (m/s)  v (m/s)  pres(mb) temp(K) RH(%)'
      WRITE(6,'(1x,5f8.2)') (qsrcsng(ivar,isrc),ivar=1,5)
      qsrcsng(3,isrc)=100.*qsrcsng(3,isrc)
      qsrcsng(5,isrc)=0.01*qsrcsng(5,isrc)
      CLOSE(4)
      isrc=isrc+1
    END IF
  END DO

3 CONTINUE
!
  isrc=1
  DO jsrc=1,nsrc_ua
    IF(srcua(jsrc) /= 'NULL') THEN
      PRINT *, 'Reading ', TRIM(uaerrfil(jsrc))
      OPEN(4,ERR=6,FILE=trim(uaerrfil(jsrc)),STATUS='old')
      READ(4,'(a8)',ERR=6,END=6) rdsource
      IF(rdsource /= srcua(jsrc)) THEN
        WRITE(6,'(a,i4,a,a,a,a,a)')                                     &
            ' Mismatch of source names',jsrc,                           &
            ' read ',rdsource,' expected ',srcua(jsrc)
        STOP
      END IF
      READ(4,'(a80)',ERR=6,END=6) srcua_full(isrc)
      READ(4,*,ERR=6,END=6) qcmulua(isrc)
      READ(4,'(a80)',ERR=6,END=6) header
      WRITE(6,'(//,a,a,/a)') 'UA data std error for ',                  &
                            srcua(isrc),srcua_full(isrc)
      WRITE(6,'(a,f5.1)') ' QC multiplier: ',qcmulua(isrc)
      WRITE(6,'(1x,a)')                                                 &
          '   k    hgt(m)  u(m/s)  v(m/s) pres(mb) temp(K)  RH(%)'
      DO ktab=1,nz_tab
        READ(4,*,ERR=5,END=5) huaqsrc(ktab,isrc),                       &
                qsrcua(3,ktab,isrc),                                    &
                qsrcua(4,ktab,isrc),                                    &
                qsrcua(5,ktab,isrc),                                    &
                qsrcua(1,ktab,isrc),                                    &
                qsrcua(2,ktab,isrc)
        WRITE(6,'(1x,i4,f10.2,5f8.2)') ktab,huaqsrc(ktab,isrc),         &
               (qsrcua(ivar,ktab,isrc),ivar=1,5)
        qsrcua(3,ktab,isrc)=100.*qsrcua(3,ktab,isrc)
        qsrcua(5,ktab,isrc)=0.01*qsrcua(5,ktab,isrc)
      END DO
      5     nlvuatab(isrc)=ktab-1
      CLOSE(4)
      isrc=isrc+1
    END IF
  END DO

6 CONTINUE
!
  isrc=1
  DO jsrc=1,nsrc_rad
    IF(srcrad(jsrc) /= 'NULL') THEN
      PRINT *, 'Reading ', TRIM(raderrfil(jsrc))
      OPEN(4,ERR=7,FILE=trim(raderrfil(jsrc)),STATUS='old')
      READ(4,'(a8)',ERR=7,END=7) rdsource
      IF(rdsource /= srcrad(jsrc)) THEN
        WRITE(6,'(a,i4,a,a,a,a,a)')                                     &
            ' Mismatch of source names',jsrc,                           &
            ' read ',rdsource,' expected ',srcrad(jsrc)
        STOP
      END IF
      READ(4,'(a80)',ERR=7,END=7) srcrad_full(isrc)
      READ(4,*,ERR=7,END=7) qcmulrad(isrc)
      READ(4,'(a80)',ERR=7,END=7) header
      READ(4,*,ERR=7,END=7)                                             &
             qsrcrad(1,isrc),                                           &
             qsrcrad(2,isrc),                                           &
             qsrcrad(3,isrc)
      WRITE(6,'(//,a,a,/a)') 'Radar data std error for ',               &
                          srcrad(isrc),srcrad_full(isrc)
      WRITE(6,'(a,f5.1)') ' QC multiplier: ',qcmulrad(isrc)
      WRITE(6,'(1x,a)')                                                 &
          '    ref(dBz) Vrad(m/s) SpWid(m/s)'
      WRITE(6,'(1x,4f8.2)')                                             &
          (qsrcrad(ivar,isrc),ivar=1,nvar_radin)
      CLOSE(4)
      isrc=isrc+1
    END IF
  END DO

7 CONTINUE
!
  isrc=1
  DO jsrc=1,nsrc_ret
    IF(srcret(jsrc) /= 'NULL') THEN
      PRINT *, 'Reading ', TRIM(reterrfil(jsrc))
      OPEN(4,FILE=trim(reterrfil(jsrc)),ERR=10,STATUS='old')
      READ(4,'(a8)',ERR=10,END=10) rdsource
      IF(rdsource /= srcret(jsrc)) THEN
        WRITE(6,'(a,i4,a,a,a,a,a)')                                     &
            ' Mismatch of source names',jsrc,                           &
            ' read ',rdsource,' expected ',srcret(jsrc)
        STOP
      END IF
      READ(4,'(a80)',ERR=10,END=10) srcret_full(isrc)
      READ(4,*,ERR=10,END=10) qcmulret(isrc)
      READ(4,'(a80)',ERR=10,END=10) header
      WRITE(6,'(//,a,a,/a)') 'Retrieval std error for ',                &
                          srcret(isrc),srcret_full(isrc)
      WRITE(6,'(a,f5.1)') ' QC multiplier: ',qcmulrad(isrc)
      WRITE(6,'(1x,a)')                                                 &
          '   k    hgt(m)  u(m/s)  v(m/s) pres(mb) temp(K)  RH(%)'
      DO ktab=1,nz_tab
        READ(4,*,END=9,ERR=9) hrtqsrc(ktab,isrc),                       &
                qsrcret(3,ktab,isrc),                                   &
                qsrcret(4,ktab,isrc),                                   &
                qsrcret(5,ktab,isrc),                                   &
                qsrcret(1,ktab,isrc),                                   &
                qsrcret(2,ktab,isrc)
        WRITE(6,'(1x,i4,f10.2,5f8.2)') ktab,hrtqsrc(ktab,isrc),         &
                      (qsrcret(ivar,ktab,isrc),ivar=1,5)
        qsrcret(3,ktab,isrc)=100.*qsrcret(3,ktab,isrc)
        qsrcret(5,ktab,isrc)=0.01*qsrcret(5,ktab,isrc)
      END DO
      9     nlvrttab(isrc)=ktab-1
      CLOSE(4)
      isrc=isrc+1
    END IF
  END DO

10 CONTINUE
!
!-----------------------------------------------------------------------
!
!  Change standard error to standard error variance by squaring
!  Calculate quality control thresholds
!
!-----------------------------------------------------------------------
!
  DO ktab=1,nlvqback
    DO ivar=1,nvar_anx
      qback(ivar,ktab)=qback(ivar,ktab)*qback(ivar,ktab)
    END DO
  END DO
!
  DO isrc=1,nsrc_sng
    DO ivar=1,nvar_anx
      qsrcsng(ivar,isrc)=qsrcsng(ivar,isrc)*qsrcsng(ivar,isrc)
      qcthrsng(ivar,isrc)=qcmulsng(isrc)*                               &
                       SQRT(qback(ivar,1)+qsrcsng(ivar,isrc))
    END DO
  END DO
!
  DO isrc=1,nsrc_ua
    DO ivar=1,nvar_anx
      qsrcmax=0.
      DO ktab=1,nlvuatab(isrc)
        qsrcua(ivar,ktab,isrc) =                                        &
            qsrcua(ivar,ktab,isrc)*qsrcua(ivar,ktab,isrc)
        qsrcmax=AMAX1(qsrcmax,qsrcua(ivar,ktab,isrc))
      END DO
      qcthrua(ivar,isrc)=qcmulua(isrc)*                                 &
                         SQRT(qback(ivar,1)+qsrcmax)
    END DO
  END DO
!
  DO isrc=1,nsrc_rad
    DO ivar=1,nvar_radin
      qsrcrad(ivar,isrc) =                                              &
          qsrcrad(ivar,isrc)*qsrcrad(ivar,isrc)
      qcthrrad(ivar,isrc)=qcmulrad(isrc)*                               &
                          SQRT(qback(ivar,1)+qsrcrad(ivar,isrc))
    END DO
  END DO
!
  DO isrc=1,nsrc_ret
    DO ivar=1,nvar_anx
      qsrcmax=0.
      DO ktab=1,nlvrttab(isrc)
        qsrcret(ivar,ktab,isrc) =                                       &
            qsrcret(ivar,ktab,isrc)*qsrcret(ivar,ktab,isrc)
        qsrcmax=AMAX1(qsrcmax,qsrcret(ivar,ktab,isrc))
      END DO
      qcthrret(ivar,isrc)=qcmulret(isrc)*                               &
                          SQRT(qback(ivar,1)+qsrcmax)
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Set cross-correlations between numerical categories.
!  Roughly 1=clear,2=some evidence of outflow,3=precip,4=conv precip
!  This could be read in as a table.
!
!-----------------------------------------------------------------------
!
  DO j=1,ncat
    DO i=1,ncat
      xcor(i,j)=1.0-(IABS(i-j))*0.25
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Initialize grids, forming first guess fields based on
!  the model initial option specified in the ARPS init file.
!
!-----------------------------------------------------------------------
!
  CALL initgrdvar(nx,ny,nz,1,nstyps,1,                                  &
                  x,y,z,zp,hterain,mapfct,                              &
                  j1,j2,j3,aj3x,aj3y,aj3z,j3inv,                        &
                  u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,         &
                  udteb, udtwb, vdtnb, vdtsb,                           &
                  pdteb,pdtwb ,pdtnb ,pdtsb,                            &
                  trigs1,trigs2,ifax1,ifax2,                            &
                  wsave1,wsave2,vwork1,vwork2,                          &
                  ubar,vbar,ptbar,pbar,tem10,tem10,                     &
                  rhostr,tem10,qvbar,ppi,csndsq,                        &
                  soiltyp,stypfrct,vegtyp,lai,roufns,veg,               &
                  tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc,       &
                  ptcumsrc,qcumsrc,w0avg,nca,kfraincv,                  &
                  cldefi,xland,bmjraincv,                               &
                  raing,rainc,prcrate,exbcbuf,                          &
                  radfrc,radsw,rnflx,                                   &
                  usflx,vsflx,ptsflx,qvsflx,                            &
                  tem1,tem2,tem3,tem4,tem5,                             &
                  tem6,tem7,tem8,tem9)
!
!-----------------------------------------------------------------------
!
!  Deallocate some arrays needed only for initgrdvar
!
!-----------------------------------------------------------------------
!
  DEALLOCATE(ptcumsrc)
  DEALLOCATE(qcumsrc)
  DEALLOCATE(w0avg)
  DEALLOCATE(kfraincv)
  DEALLOCATE(nca)
  DEALLOCATE(cldefi)
  DEALLOCATE(xland)
!
!-----------------------------------------------------------------------
!
!  Set location of scalar fields.
!
!-----------------------------------------------------------------------
!
  CALL xytoll(nx,ny,x,y,latgr,longr)

  DO i=1,nx-1
    xs(i)=0.5*(x(i)+x(i+1))
  END DO
  xs(nx)=2.0*xs(nx-1)-xs(nx-2)
  DO j=1,ny-1
    ys(j)=0.5*(y(j)+y(j+1))
  END DO
  ys(ny)=2.0*ys(ny-1)-ys(ny-2)
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        zs(i,j,k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
      END DO
    END DO
  END DO
  DO i=1,nx-1
    DO j=1,ny-1
      zs(i,j,nz)=2.*zs(i,j,nz-1)-zs(i,j,nz-2)
    END DO
  END DO

  CALL ctim2abss( year,month,day,hour,minute,second, i4timanx)
!
!-----------------------------------------------------------------------
!
!  Identify the background field correlation category for
!  each 2-d point.   This is based on precip rate, cumulus
!  parameterization switch, and surface relative humidity.
!
!-----------------------------------------------------------------------
!
  CALL setcat(nx,ny,nz,ccatopt,zs,                                      &
              ptprt,pprt,qv,qc,qr,qi,qs,qh,                             &
              ptbar,pbar,                                               &
              prcrate,icatg)
!
!-----------------------------------------------------------------------
!
!  Load "background fields" for analysis
!  Note, analysis is done at scalar points
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        anx(i,j,k,1)=0.5*(u(i,j,k)+u(i+1,j,k))
        anx(i,j,k,2)=0.5*(v(i,j,k)+v(i,j+1,k))
        anx(i,j,k,3)=pbar(i,j,k)+pprt(i,j,k)
        anx(i,j,k,4)=ptbar(i,j,k)+ptprt(i,j,k)
        anx(i,j,k,5)=qv(i,j,k)
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Save background fields on staggered grid in the increment
!  arrays.
!
!-----------------------------------------------------------------------
!
  IF(incrdmp > 0) THEN
    CALL incrsave(nx,ny,nz,nxndg,nyndg,nzndg,                           &
                  u,v,w,pprt,ptprt,qv,                                  &
                  qc,qr,qi,qs,qh,                                       &
                  uincr,vincr,wincr,pincr,ptincr,qvincr,                &
                  qcincr,qrincr,qiincr,qsincr,qhincr,                   &
                  istatus)
  END IF
!
!-----------------------------------------------------------------------
!
!  Filename creation
!
!-----------------------------------------------------------------------
!
  WRITE(stats,'(a,a)') runname(1:lfnkey),'.lst'
  CALL getunit(iwstat)
  PRINT *, 'Writing ', TRIM(stats)
  OPEN(iwstat,IOSTAT=ios,FILE=trim(stats),STATUS='unknown')
  IF(ios /= 0) iwstat=0

  maxsuf=LEN(suffix)
  mxroot=LEN(froot)
  nobsng=0
  DO ifile=1,nsngfil
    PRINT *, 'Processing file ', ifile, ' ', TRIM(sngfname(ifile))
    lenfnm=LEN(sngfname(ifile))
    CALL strlnth(sngfname(ifile),lenfnm)
    CALL exsufx(sngfname(ifile),lenfnm,suffix,maxsuf,dotloc,lensuf)
    IF(lensuf == 3. .AND. suffix(1:3) == 'lso') THEN

      CALL exfroot(sngfname(ifile),lenfnm,froot,mxroot,lenroot)
      WRITE(qclist,'(a,a)') froot(1:lenroot),'.sqc'
      WRITE(qcsave,'(a,a)') froot(1:lenroot),'.lsq'
!
!-----------------------------------------------------------------------
!
!  Open the files for listing QC info
!  To suppress listing set the unit numbers to zero
!
!-----------------------------------------------------------------------
!
      CALL getunit(iqclist)
      PRINT *, 'Opening qclist: ', TRIM(qclist)
      OPEN(iqclist,IOSTAT=ios,FILE=trim(qclist),STATUS='unknown')
      IF(ios /= 0) iqclist=0
      CALL getunit(iqcsave)
      PRINT *, 'Opening qclist: ', TRIM(qcsave)
      OPEN(iqcsave,IOSTAT=ios,FILE=trim(qcsave),STATUS='unknown')
      IF(ios /= 0) iqcsave=0
!
!-----------------------------------------------------------------------
!
!  Read surface data, QC and convert units.
!
!-----------------------------------------------------------------------
!
      sprkm=0.001*sprdist
      rngsq=sfcqcrng*sfcqcrng
      PRINT *, 'Calling prepsfcobs'
      CALL prepsfcobs(ntime,mx_sng,                                     &
          nvar_sng,nvar_anx,nsrc_sng,nobsng,                            &
          sngfname(ifile),sngtmchk(ifile),blackfil,                     &
          var_sfc,var_anx,srcsng,qsrcsng,                               &
          rmiss,iqspr,sprkm,nobsrd,timesng,                             &
          stnsng,csrcsng,isrcsng,latsng,lonsng,hgtsng,xsng,ysng,        &
          wx,kloud,idp3,store_emv,store_hgt,store_amt,                  &
          obrdsng,obsng,qualrdsng,qobsng,qualsng,                       &
          ival,climin,climax,                                           &
          rngsq,klim,wlim,qcthrsng,                                     &
          knt,wgtsum,zsum,sqsum,iqclist,iqcsave,jstatus)

      IF(iqclist /= 0) THEN
        CLOSE(iqclist)
        CALL retunit(iqclist)
      END IF
      IF(iqcsave /= 0) THEN
        CLOSE(iqcsave)
        CALL retunit(iqcsave)
      END IF

    ELSE
!
!-----------------------------------------------------------------------
!
!  Read other single-level data and convert units.
!
!-----------------------------------------------------------------------
!
      PRINT *, 'Calling prepsglobs'
      CALL prepsglobs(mx_sng,ntime,srcsng,nsrc_sng,                     &
              sngfname(ifile),stnsng,latsng,lonsng,xsng,ysng,           &
              hgtsng,obsng,qobsng,qualsng,isrcsng,qsrcsng,              &
              csrcsng,nx,ny,nz,nvar_anx,anx,xs,ys,zp,                   &
              tem1,tem2,tem3,tem4,tem5,tem6,                            &
              rmiss,nobsng,istatus)

    END IF
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate initial obs differences for single-level data
!
!-----------------------------------------------------------------------
!
  PRINT *, 'Calling grdtosng'
  CALL grdtosng(nx,ny,nz,nz_tab,mx_sng,nvar_anx,nobsng,                 &
                    xs,ys,zp,icatg,anx,qback,hqback,nlvqback,           &
                    tem1,tem2,tem3,tem4,tem5,tem6,                      &
                    stnsng,isrcsng,icatsng,hgtsng,xsng,ysng,            &
                    obsng,qobsng,qualsng,                               &
                    odifsng,oanxsng,thesng)
!
!-----------------------------------------------------------------------
!
!  Read upper-air data, QC and convert units
!
!-----------------------------------------------------------------------
!
  PRINT *, 'Calling prepuaobs'
  CALL prepuaobs(nvar_anx,nz_ua,mx_ua,nz_tab,nsrc_ua,                   &
                    mx_ua_file,nuafil,uafname,srcua,                    &
                    stnua,elevua,xua,yua,hgtua,obsua,                   &
                    qsrcua,huaqsrc,nlvuatab,                            &
                    qobsua,qualua,isrcua,nlevsua,                       &
                    rmiss,nobsua,istatus)
!
!-----------------------------------------------------------------------
!
!  Calculate initial obs differences for upper-air data
!
!-----------------------------------------------------------------------
!
  PRINT *, 'Calling grdtoua'
  CALL grdtoua(nx,ny,nz,nz_tab,nz_ua,mx_ua,nvar_anx,nobsua,             &
              xs,ys,zp,anx,qback,hqback,nlvqback,                       &
              tem1,tem2,tem3,tem4,tem5,tem6,                            &
              stnua,isrcua,elevua,xua,yua,hgtua,                        &
              obsua,qobsua,qualua,nlevsua,                              &
              odifua,oanxua,theua)
!
!-----------------------------------------------------------------------
!
!  Read radar data, unfold and convert units
!
!-----------------------------------------------------------------------
!
  PRINT *, 'Calling prepradar'
  CALL prepradar(nx,ny,nz,nz_tab,nvar_anx,nvar_radin,                   &
             nvar_rad,mx_rad,nsrc_rad,nz_rdr,mx_colrad,mx_pass,         &
             raduvobs,radrhobs,radistride,radkstride,                   &
             iuserad,npass,refrh,rhradobs,                              &
             xs,ys,zs,hterain,anx,qback,hqback,nlvqback,                &
             nradfil,radfname,srcrad,isrcrad,qsrcrad,qcthrrad,          &
             stnrad,latrad,lonrad,elvrad,                               &
             latradc,lonradc,xradc,yradc,irad,nlevrad,                  &
             distrad,uazmrad,vazmrad,hgtradc,theradc,                   &
             obsrad,oanxrad,odifrad,qobsrad,qualrad,                    &
             ncolrad,istatus,tem1,tem2,tem3,tem4,tem5,tem6)

!
!-----------------------------------------------------------------------
!
!  Read retrieval data.
!
!-----------------------------------------------------------------------
!
  PRINT *, 'Calling prepretr'
  CALL prepretr(nx,ny,nz,nvar_anx,                                      &
                nz_ret,mx_ret,mx_colret,nz_tab,nsrc_ret,                &
                nretfil,retfname,                                       &
                isrcret,srcret,nlvrttab,qsrcret,hrtqsrc,                &
                stnret,latret,lonret,elvret,                            &
                latretc,lonretc,xretc,yretc,iret,nlevret,               &
                hgtretc,obsret,qrret,qobsret,qualret,                   &
                rmiss,ncolret,tem1,istatus)

  PRINT *, 'Calling grdtoret'
  CALL grdtoret(nx,ny,nz,nz_tab,                                        &
                nz_ret,mx_ret,mx_colret,nvar_anx,ncolret,               &
                xs,ys,zp,anx,qback,hqback,nlvqback,                     &
                tem1,tem2,tem3,tem4,tem5,tem6,                          &
                stnret,iret,xretc,yretc,hgtretc,                        &
                obsret,qobsret,qualret,nlevret,                         &
                odifret,oanxret,theretc)
!
!-----------------------------------------------------------------------
!
!  Quality-control observation differences
!
!-----------------------------------------------------------------------
!
  PRINT *, 'Calling qcdiff'
  CALL qcdiff(nvar_anx,nvar_rad,nvar_radin,mx_sng,nsrc_sng,             &
              nz_ua,mx_ua,nsrc_ua,                                      &
              nz_rdr,mx_rad,mx_colrad,nsrc_rad,                         &
              nz_ret,mx_ret,mx_colret,nsrc_ret,                         &
              nobsng,nobsua,ncolrad,ncolret,var_anx,                    &
              stnsng,isrcsng,hgtsng,obsng,oanxsng,odifsng,              &
              qcthrsng,qualsng,                                         &
              stnua,isrcua,hgtua,obsua,oanxua,odifua,                   &
              qcthrua,qualua,nlevsua,                                   &
              stnrad,irad,isrcrad,hgtradc,obsrad,odifrad,               &
              qcthrrad,qualrad,nlevrad,                                 &
              stnret,iret,isrcret,hgtretc,                              &
              obsret,oanxret,odifret,                                   &
              qcthrret,qualret,nlevret,                                 &
              wgtsum,zsum)
!
!-----------------------------------------------------------------------
!
!  Analyze data
!
!-----------------------------------------------------------------------
!
!wdt
  CALL FLUSH (06)
  PRINT *, 'Calling anxiter'
  CALL anxiter(nx,ny,nz,                                                &
          nvar_anx,nvar_radin,nvar_rad,nz_ua,nz_rdr,nz_ret,             &
          mx_sng,mx_ua,mx_rad,mx_colrad,mx_ret,mx_colret,               &
          nsrc_sng,nsrc_ua,nsrc_rad,nsrc_ret,ncat,                      &
          mx_pass,npass,iwstat,xs,ys,zs,icatg,xcor,var_anx,             &
          xsng,ysng,hgtsng,thesng,                                      &
          obsng,odifsng,qobsng,qualsng,isrcsng,icatsng,nobsng,          &
          xua,yua,hgtua,theua,                                          &
          obsua,odifua,qobsua,qualua,isrcua,nlevsua,nobsua,             &
          elvrad,xradc,yradc,                                           &
          distrad,uazmrad,vazmrad,hgtradc,theradc,                      &
          obsrad,odifrad,qobsrad,qualrad,                               &
          irad,isrcrad,nlevrad,ncolrad,                                 &
          xretc,yretc,hgtretc,theretc,                                  &
          obsret,odifret,qobsret,qualret,                               &
          iret,isrcret,nlevret,ncolret,                                 &
          srcsng,srcua,srcrad,srcret,                                   &
          ianxtyp,iusesng,iuseua,iuserad,iuseret,                       &
          xyrange,kpvrsq,wlim,zrange,zwlim,                             &
          thrng,rngsqi,knt,wgtsum,zsum,                                 &
          corsng,corua,corrad,corret,                                   &
          oanxsng,oanxua,oanxrad,oanxret,                               &
          anx,tem1,tem2,tem3,ibegin,iend,istatus)
!
!wdt
!-----------------------------------------------------------------------
!
!	Write to ARPSControl status file
!
!-----------------------------------------------------------------------
! Only u (index 1) is checked. Eventually should check all
! Only 1 level of UA data is checked.
!
! 10   Good
! 100  Good
! 200  Good - superob
! -9   Outside domain
! -99  Initializing value
! -199 Rejected in quality test
!
!  > 0    : used++ : in file, in domain, good point
!  = -199 : rej++  : in file, in domain, bad point
!  = -9   :        : in file, outside domain
!  = 0    :        : not in file
!
! Good data has isrc* > 0 and qual* > 0.

! single-point: SA (sao), MESO (okmeso)

  WRITE (status_file, '(2A)') runname(1:lfnkey),'.adasstat'
  PRINT *, 'Creating ', TRIM(status_file)
  OPEN (UNIT=09, FILE=status_file)

  n_used_sng = 0
  DO jsta=1,nobsng
    IF (isrcsng(jsta) .GT. 0 .AND. qualsng(1,jsta) .GT. 0) THEN
      n_used_sng(isrcsng(jsta)) = n_used_sng(isrcsng(jsta)) + 1
    ELSE
    END IF
  END DO

  DO jsrc=1,nsrc_sng
    IF (srcsng(jsrc)(:2) .EQ. 'SA') THEN
      WRITE (09,9901) '$sao_n_used = ', n_used_sng(jsrc)
    ELSE IF (srcsng(jsrc)(:4) .EQ. 'MESO') THEN
      WRITE (09,9901) '$meso_n_used = ', n_used_sng(jsrc)
    ELSE IF (srcsng(jsrc)(:5) .EQ. 'MDCRS') THEN
      WRITE (09,9901) '$mdcrs_n_used = ', n_used_sng(jsrc)
!wdt update 2002-03-05 from kwthomas -- added buoy
    ELSE IF (srcsng(jsrc)(:4) .EQ. 'BUOY') THEN
      WRITE (09,9901) '$buoy_n_used = ', n_used_sng(jsrc)
    ELSE IF (srcsng(jsrc)(:4) .NE. 'NULL') THEN
      WRITE (09,'(3A)') '# Unable to process srcsng: ', &
			srcsng(jsrc), stnsng(jsta,1)
    END IF
  END DO

! upper-air: NWS RAOB (raob), WPDN PRO (pro)

  n_used_ua = 0
  DO jsta=1,nobsua
    IF ((isrcua(jsta) <= 0) .or. (isrcua(jsta) > nsrc_ua)) THEN
      WRITE (09,*) '# Unable to process: irsrcua=',  &
         isrcua(jsta), ' ', TRIM(stnua(jsta))
    ELSE IF (srcua(isrcua(jsta)) .EQ. 'NWS RAOB') THEN
      WRITE (09,9902) 'raob', TRIM(stnua(jsta)), nlevsua(jsta)
      n_used_ua(isrcua(jsta)) = n_used_ua(isrcua(jsta)) + 1
    ELSE IF (srcua(isrcua(jsta)) .EQ. 'WPDN PRO') THEN
      WRITE (09,9902) 'pro', TRIM(stnua(jsta)), nlevsua(jsta)
      n_used_ua(isrcua(jsta)) = n_used_ua(isrcua(jsta)) + 1
    ELSE IF (srcua(isrcua(jsta))(:4) .NE. 'NULL') THEN
      WRITE (09,'(3A)') '# Unable to process: ',  &
         srcua(isrcua(jsta)), TRIM(stnua(jsta))
    ENDIF
  END DO

!depricated:
!  DO jsta=1,nobsua
!! criterion: 3 good u obs out of the first 5 levels
!    IF (isrcua(jsta) .GT. 0) THEN
!      ngoodlev = 0
!
!      DO klev=1,5	!nlevsua(jsta)
!	IF (qualua(1,klev,jsta) .GT. 0) ngoodlev = ngoodlev + 1
!      END DO
!
!      IF (ngoodlev .GE. 3) THEN
!	IF (srcua(isrcua(jsta)) .EQ. 'NWS RAOB') THEN
!	  WRITE (09,9902) 'raob', TRIM(stnua(jsta)), nlevsua(jsta)
!	ELSE IF (srcua(isrcua(jsta)) .EQ. 'WPDN PRO') THEN
!	  WRITE (09,9902) 'pro', TRIM(stnua(jsta)), nlevsua(jsta)
!	ELSE IF (srcua(jsrc)(:4) .NE. 'NULL') THEN
!	  WRITE (09,'(3A)') '# Unable to process: ', &
!	                     srcua(jsrc), TRIM(stnua(jsta))
!	END IF
!	n_used_ua(isrcua(jsta)) = n_used_ua(isrcua(jsta)) + 1
!      END IF
!
!    END IF ! good isrcua
!  END DO

  DO jsrc=1,nsrc_ua
    IF (srcua(jsrc)(:8) .EQ. 'NWS RAOB') THEN
      WRITE (09,9901) '$raob_n_used = ', n_used_ua(jsrc)
    ELSE IF (srcua(jsrc)(:8) .EQ. 'WPDN PRO') THEN
      WRITE (09,9901) '$pro_n_used  = ', n_used_ua(jsrc)
    END IF
  END DO

  9902 FORMAT ('${', A, '}{stn_levels}{"', A, '"} = ', I6, ';')
  9901 FORMAT (A,I6,';')
! wdt 
  CALL FLUSH (06)
  CLOSE (09)

!-----------------------------------------------------------------------
!
!  Use reflectivity to establish cloud water
!
!-----------------------------------------------------------------------
!
  IF( radcldopt > 0 ) CALL radmcro(nx,ny,nz,                            &
                 mx_rad,nsrc_rad,nvar_radin,nz_rdr,mx_colrad,mx_pass,   &
                 radistride,radkstride,iuserad,npass,                   &
                 nradfil,radfname,srcrad,isrcrad,stnrad,                &
                 latrad,lonrad,elvrad,latradc,lonradc,irad,             &
                 nlevrad,xradc,yradc,hgtradc,obsrad,ncolrad,            &
                 radqvopt,radqcopt,radqropt,radptopt,                   &
                 refsat,rhrad,                                          &
                 refcld,cldrad,ceilopt,ceilmin,dzfill,                  &
                 refrain,radsetrat,radreflim,radptgain,                 &
                 xs,ys,zs,zp,j3,                                        &
                 anx(1,1,1,3),anx(1,1,1,4),anx(1,1,1,5),                &
                 ptbar,qvbar,rhostr,                                    &
                 qc,qr,qi,qs,qh,                                        &
                 tem1,tem2,tem3,tem4,tem5,tem6,tem7)
!
!-----------------------------------------------------------------------
!
!  Transfer analysed fields to ARPS variables for writing.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        pprt(i,j,k)=anx(i,j,k,3)-pbar(i,j,k)
        ptprt(i,j,k)=anx(i,j,k,4)-ptbar(i,j,k)
        tgrid =anx(i,j,k,4) * (anx(i,j,k,3)/p0) ** rddcp
        qvsat=f_qvsat( anx(i,j,k,3),tgrid)
        qvmin=rhmin*qvsat
        qv(i,j,k)=MAX(qvmin,MIN(anx(i,j,k,5),qvsat))
      END DO
    END DO
  END DO

  IF(spradopt == 1) THEN

!-----------------------------------------------------------------------
!
!  Check for superadiabatic layers in each column.
!  For spradopt=1,
!  Where superadiabatic levels are found, make the next level down
!  cooler, applying the minimum potential temperature lapse rate,
!  ptlapse.
!
!-----------------------------------------------------------------------
!
    DO j=1,ny-1
      DO i=1,nx-1
        ptprt(i,j,nz-1)=anx(i,j,nz-1,4)-ptbar(i,j,nz-1)
        DO k=nz-2,2,-1
          anx(i,j,k,4)=MIN(anx(i,j,k,4),                                &
                (anx(i,j,k+1,4)-ptlapse*(zs(i,j,k+1)-zs(i,j,k))))
          ptprt(i,j,k)=anx(i,j,k,4)-ptbar(i,j,k)
        END DO
      END DO
    END DO

  ELSE IF(spradopt == 2) THEN
!
!-----------------------------------------------------------------------
!
!  For spradopt=2,
!  Where superadiabatic levels are found, make the next level up
!  warmer, applying the minimum potential temperature lapse rate,
!  ptlapse.
!
!-----------------------------------------------------------------------
!
    DO j=1,ny-1
      DO i=1,nx-1
        ptprt(i,j,2)=anx(i,j,2,4)-ptbar(i,j,2)
        DO k=3,nz-1
          anx(i,j,k,4)=MAX(anx(i,j,k,4),                                &
                (anx(i,j,k-1,4)+ptlapse*(zs(i,j,k)-zs(i,j,k-1))))
          ptprt(i,j,k)=anx(i,j,k,4)-ptbar(i,j,k)
        END DO
      END DO
    END DO

  END IF
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=2,nx-1
        u(i,j,k)=0.5*(anx(i,j,k,1)+anx(i-1,j,k,1))
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      u(1,j,k)=u(2,j,k)
      u(nx,j,k)=u(nx-1,j,k)
    END DO
  END DO
!
  DO k=1,nz-1
    DO j=2,ny-1
      DO i=1,nx-1
        v(i,j,k)=0.5*(anx(i,j,k,2)+anx(i,j-1,k,2))
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO i=1,nx-1
      v(i,1,k)=v(i,2,k)
      v(i,ny,k)=v(i,ny-1,k)
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Adjust wind fields to minimize 3-d divergence
!  This is the first run-through for "w" in order to obtain
!  a good estimate of w for the cloud analysis.  Note here wndadj
!  is set to "2" for integration of 2-d divergence for "w" and
!  the cloud_w opt is set to 0.
!
!-----------------------------------------------------------------------
!
  CALL adjuvw( nx,ny,nz,u,v,w,                                          &
               wcont,ptprt,ptbar,                                       &
               zp,j1,j2,j3,aj3z,mapfct,rhostr,w_cld,                    &
               2,obropt,obrzero,0,                                      &
               tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)
!
!-----------------------------------------------------------------------
!
!  Enforce vertical w boundary conditions
!
!-----------------------------------------------------------------------
!
  CALL rhouvw(nx,ny,nz,rhostr,tem1,tem2,tem3)
  CALL vbcw(nx,ny,nz,w,wcont,tbc,bbc,u,v,                               &
              rhostr,tem1,tem2,tem3,                                    &
              j1,j2,j3)
!
!-----------------------------------------------------------------------
!
!  Free-up some memory.
!
!-----------------------------------------------------------------------
!
  DEALLOCATE(anx)
!
!-----------------------------------------------------------------------
!
!  Call cloud analysis.
!
!-----------------------------------------------------------------------
!
  IF( cloudopt > 0 ) THEN
    ALLOCATE(ref_mos_3d(nx,ny,nz),STAT=istatus)
    CALL check_alloc_status(istatus, "adas:ref_mos_3d")
    ref_mos_3d=0.

    ALLOCATE(w_cld(nx,ny,nz),STAT=istatus)
    CALL check_alloc_status(istatus, "adas:w_cld")
    w_cld=0.
    ALLOCATE(cldpcp_type_3d(nx,ny,nz),STAT=istatus)
    CALL check_alloc_status(istatus, "adas:cldpcp_type_3d")
    cldpcp_type_3d=0
    ALLOCATE(icing_index_3d(nx,ny,nz),STAT=istatus)
    CALL check_alloc_status(istatus, "adas:icing_index_3d")
    icing_index_3d=0
    ALLOCATE(l_mask_pcptype(nx,ny),STAT=istatus)
    CALL check_alloc_status(istatus, "adas:l_mask_pcptype")
    l_mask_pcptype=.false.

    CALL cmpclddrv(nx,ny,nz,i4timanx,                  &
              xs,ys,zs,j3,hterain,latgr,longr,                          &
              nobsng,stnsng,csrcsng,xsng,ysng,                          &
              timesng,latsng,lonsng,hgtsng,                             &
              kloud,store_amt,store_hgt,                                &
              stnrad,isrcrad,latrad,lonrad,elvrad,ixrad,jyrad,          &
              pprt,ptprt,qv,qc,qr,qi,qs,qh,w,                           &
              pbar,ptbar,qvbar,rhostr,                                  &
              ref_mos_3d,w_cld,cldpcp_type_3d,                          &
              icing_index_3d,l_mask_pcptype,                            &
              tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9,tem10)
  END IF
!
  IF(retqropt > 0 .AND. ncolret > 0) CALL retmcro(nx,ny,nz,             &
                 mx_ret,nsrc_ret,nvar_anx,nz_ret,mx_colret,             &
                 srcret,isrcret,iret,nlevret,                           &
                 xretc,yretc,hgtretc,qrret,ncolret,                     &
                 dzfill,                                                &
                 xs,ys,zs,qr)
!
!-----------------------------------------------------------------------
!
!  Create rhobar from rhostr
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem6(i,j,k)=rhostr(i,j,k)/j3(i,j,k)
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate and store the sound wave speed squared in csndsq.
!  Use original definition of sound speed.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        csndsq(i,j,k)= cpdcv*pbar(i,j,k)/tem6(i,j,k)
      END DO
    END DO
  END DO
!
  IF(hydradj == 1 .OR. hydradj == 2) THEN
    pconst=0.5*g*dz
!
!-----------------------------------------------------------------------
!
!  Create thermal buoyancy at each scalar point,
!  which is stored in tem2
!
!-----------------------------------------------------------------------
!
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          qvprt=qv(i,j,k)-qvbar(i,j,k)
          qtot=qc(i,j,k)+qr(i,j,k)+                                     &
               qi(i,j,k)+qs(i,j,k)+qh(i,j,k)
          tem2(i,j,k)=j3(i,j,k)*tem6(i,j,k)*                            &
                    g*( (ptprt(i,j,k)/ptbar(i,j,k)) +                   &
                        (qvprt/(rddrv+qvbar(i,j,k))) -                  &
                        ((qvprt+qtot)/(1.+qvbar(i,j,k))) )
        END DO
      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  Average thermal buoyancy to w points
!  which is stored in tem3
!
!-----------------------------------------------------------------------
!
    CALL avgsw(tem2,nx,ny,nz,1,nx,1,ny, tem3)

    IF(hydradj == 1) THEN

      DO i=1,nx
        DO j=1,ny
          tem1(i,j,1)=pprt(i,j,1)
          DO k=2,nz-2
            tem1(i,j,k)=                                                &
                  ( (1. - (pconst*j3(i,j,k-1)/csndsq(i,j,k-1)) )*       &
                      tem1(i,j,k-1) +                                   &
                      dz*tem3(i,j,k) ) /                                &
                      (1. + (pconst*j3(i,j,k)/csndsq(i,j,k)) )
            IF(j == 26 .AND. MOD(i,10) == 0) THEN
              IF(k == nz-2) PRINT *, '            Point i= ',i, '  j=26'
              PRINT *, ' k,pprt,tem1,diff =',                           &
                  k,pprt(i,j,k),tem1(i,j,k),                            &
                  (tem1(i,j,k)-pprt(i,j,k))
            END IF
            pprt(i,j,k)=tem1(i,j,k)
          END DO
          pprt(i,j,nz-1)=pprt(i,j,nz-2)
        END DO
      END DO

    ELSE IF(hydradj == 2) THEN

      DO i=1,nx
        DO j=1,ny
          tem1(i,j,nz-1)=pprt(i,j,nz-1)
          DO k=nz-2,2,-1
            tem1(i,j,k)=                                                &
                  ( (1.+ (pconst*j3(i,j,k+1)/csndsq(i,j,k+1)) )*        &
                      tem1(i,j,k+1) -                                   &
                      dz*tem3(i,j,k+1) ) /                              &
                      (1.- (pconst*j3(i,j,k  )/csndsq(i,j,k  )) )
            IF(j == 26 .AND. MOD(i,10) == 0) THEN
              IF(k == nz-2) PRINT *, '            Point i= ',i, '  j=26'
              PRINT *, ' k,pprt,tem1,diff =',                           &
                  k,pprt(i,j,k),tem1(i,j,k),                            &
                  (tem1(i,j,k)-pprt(i,j,k))
            END IF
            pprt(i,j,k)=tem1(i,j,k)
          END DO
          pprt(i,j,1)=pprt(i,j,2)
        END DO
      END DO
    END IF
  ELSE IF (hydradj == 3) THEN
!
!-----------------------------------------------------------------------
!
!  Calculate total pressure, store in tem1.
!  Calculate virtual temperature, store in tem2.
!
!-----------------------------------------------------------------------
!
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          tem1(i,j,k) = pbar(i,j,k)+pprt(i,j,k)
          temp = (ptbar(i,j,k)+ptprt(i,j,k))*                           &
                 ((tem1(i,j,k)/p0)**rddcp)
          tem2(i,j,k) = temp*(1.0+rvdrd*qv(i,j,k))/                     &
                             (1.0+qv(i,j,k))
        END DO
      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  Integrate hydrostatic equation, begining at k=2
!
!-----------------------------------------------------------------------
!
    pconst=-g/rd
    DO k=3,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tvbar=0.5*(tem2(i,j,k)+tem2(i,j,k-1))
          tem1(i,j,k)=tem1(i,j,k-1)*                                    &
                     EXP(pconst*(zs(i,j,k)-zs(i,j,k-1))/tvbar)
          pprt(i,j,k)=tem1(i,j,k)-pbar(i,j,k)
        END DO
      END DO
    END DO
    DO j=1,ny-1
      DO i=1,nx-1
        tvbar=tem2(i,j,2)
        tem1(i,j,1)=tem1(i,j,2)*                                        &
                   EXP(pconst*(zs(i,j,1)-zs(i,j,2))/tvbar)
        pprt(i,j,1)=tem1(i,j,1)-pbar(i,j,1)
      END DO
    END DO
  END IF
!
!-----------------------------------------------------------------------
!
!  Adjust wind fields to minimize 3-d divergence
!
!-----------------------------------------------------------------------
!
  CALL adjuvw( nx,ny,nz,u,v,w,                                          &
               wcont,ptprt,ptbar,                                       &
               zp,j1,j2,j3,aj3z,mapfct,rhostr,w_cld,                    &
               wndadj,obropt,obrzero,cldwopt,                           &
               tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)
!
!-----------------------------------------------------------------------
!
!  Enforce vertical w boundary conditions
!
!-----------------------------------------------------------------------
!
  CALL rhouvw(nx,ny,nz,rhostr,tem1,tem2,tem3)
  CALL vbcw(nx,ny,nz,w,wcont,tbc,bbc,u,v,                               &
              rhostr,tem1,tem2,tem3,                                    &
              j1,j2,j3)
!
!-----------------------------------------------------------------------
!
!  Enforce boundary conditions
!  Use zero gradient for lateral bc, and rigid for top,bottom.
!
!-----------------------------------------------------------------------
!
  ebc = 3
  wbc = 3
  nbc = 3
  sbc = 3
  tbc = 1
  bbc = 1
!
  DO k=1,nz
    DO j=1,ny
      pdteb(j,k)=0.
      pdtwb(j,k)=0.
    END DO
    DO i=1,nx
      pdtnb(i,k)=0.
      pdtsb(i,k)=0.
    END DO
  END DO
!
  CALL bcu(nx,ny,nz,0., u,                                              &
           pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)

  CALL bcv(nx,ny,nz,0., v,                                              &
           pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)

  CALL lbcw(nx,ny,nz,0.0,w,wcont,                                       &
            pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc)

  CALL bcsclr(nx,ny,nz,0.,ptprt,ptprt,ptprt,                            &
            pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)

  CALL bcsclr(nx,ny,nz,0., qv,qv,qv,                                    &
           pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)

  CALL bcsclr(nx,ny,nz,0., qc,qc,qc,                                    &
           pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)

  CALL bcsclr(nx,ny,nz,0., qr,qr,qr,                                    &
           pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)

  CALL bcsclr(nx,ny,nz,0., qi,qi,qi,                                    &
           pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)

  CALL bcsclr(nx,ny,nz,0., qs,qs,qs,                                    &
           pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)

  CALL bcsclr(nx,ny,nz,0., qh,qh,qh,                                    &
           pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)

!
!-----------------------------------------------------------------------
!
!  Apply extrapolated gradient bc to pressure for geostrophic
!  type balance to winds along lateral boundaries.  Consistent
!  with horizontal winds constant.
!
!-----------------------------------------------------------------------
!
  DO k=2,nz-2
    DO j=1,ny-1
      pprt(1,j,k)=2.*pprt(2,j,k)-pprt(3,j,k)
    END DO
  END DO

  DO k=2,nz-2
    DO j=1,ny-1
      pprt(nx-1,j,k)=2.*pprt(nx-2,j,k)-pprt(nx-3,j,k)
    END DO
  END DO

  DO k=2,nz-2
    DO i=1,nx-1
      pprt(i,ny-1,k)=2.*pprt(i,ny-2,k)-pprt(i,ny-3,k)
    END DO
  END DO

  DO k=2,nz-2
    DO i=1,nx-1
      pprt(i,1,k)=2.*pprt(i,2,k)-pprt(i,3,k)
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  For top and bottom bc, apply hydrostatic pressure equation to
!  total pressure, then subtract pbar.
!
!-----------------------------------------------------------------------
!
  pconst=-g/rd
  DO j=1,ny-1
    DO i=1,nx-1
      pr2 = pbar(i,j,2)+pprt(i,j,2)
      temp = (ptbar(i,j,2)+ptprt(i,j,2))*                               &
             ((pr2/p0)**rddcp)
      tvbar = temp*(1.0+rvdrd*qv(i,j,2))/                               &
                         (1.0+qv(i,j,2))
      pr1=pr2*EXP(pconst*(zs(i,j,1)-zs(i,j,2))/tvbar)
      pprt(i,j,1)=pr1-pbar(i,j,1)

      pr2 = pbar(i,j,nz-2)+pprt(i,j,nz-2)
      temp = (ptbar(i,j,nz-2)+ptprt(i,j,nz-2))*                         &
             ((pr2/p0)**rddcp)
      tvbar= temp*(1.0+rvdrd*qv(i,j,nz-2))/                             &
                         (1.0+qv(i,j,nz-2))
      pr1=pr2*EXP(pconst*(zs(i,j,nz-1)-zs(i,j,nz-2))/tvbar)
      pprt(i,j,nz-1)=pr1-pbar(i,j,nz-1)
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Adjust surface (skin) temperature based on radiation and
!  current air temperature.
!
!-----------------------------------------------------------------------
!
  IF(tsfcadj > 0 ) THEN

    IF (radopt == 2) THEN
      rbufsz = n2d_radiat*nx*ny + n3d_radiat*nx*ny*nz
    ELSE
      rbufsz = 1
    END IF

    ALLOCATE(rsirbm(nx,ny))
    rsirbm = 0.
    ALLOCATE(rsirdf(nx,ny))
    rsirdf = 0.
    ALLOCATE(rsuvbm(nx,ny))
    rsuvbm = 0.
    ALLOCATE(rsuvdf(nx,ny))
    rsuvdf = 0.

    ALLOCATE(cosz(nx,ny))
    cosz = 0.
    ALLOCATE(cosss(nx,ny))
    cosss = 0.


    ALLOCATE(fdirir(nx,ny))
    fdirir = 0.
    ALLOCATE(fdifir(nx,ny))
    fdifir = 0.
    ALLOCATE(fdirpar(nx,ny))
    fdirpar= 0.
    ALLOCATE(fdifpar(nx,ny))
    fdifpar= 0.
    ALLOCATE(st4(nx,ny))
    st4 = 0.

    ALLOCATE(tem11(nx,ny,nz))
    tem11 = 0.
    ALLOCATE(tem12(nx,ny,nz))
    tem12 = 0.
    ALLOCATE(tem13(nx,ny,nz))
    tem13 = 0.
    ALLOCATE(tem14(nx,ny,nz))
    tem14 = 0.
    ALLOCATE(tem15(nx,ny,nz))
    tem15 = 0.
    ALLOCATE(tem16(nx,ny,nz))
    tem16 = 0.
    ALLOCATE(radbuf(rbufsz))
    radbuf = 0.

    CALL adjtsfc(nx,ny,nz,rbufsz,                                       &
               ptprt,pprt,qv,qc,qr,qi,qs,qh,                            &
               ptbar,pbar,ppi,rhostr,                                   &
               x,y,z,zp,j3inv,                                          &
               soiltyp, tsfc, wetsfc,snowdpth,                          &
               radfrc, radsw, rnflx,                                    &
               rsirbm,rsirdf,rsuvbm,rsuvdf,                             &
               cosz, cosss,                                             &
               fdirir,fdifir,fdirpar,fdifpar,st4,                       &
               tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,                 &
               tem9,tem10,tem11,tem12,tem13,tem14,tem15,tem16,          &
               radbuf)

    DEALLOCATE(rsirbm)
    DEALLOCATE(rsirdf)
    DEALLOCATE(rsuvbm)
    DEALLOCATE(rsuvdf)

    DEALLOCATE(cosz)
    DEALLOCATE(cosss)

    DEALLOCATE(fdirir)
    DEALLOCATE(fdifir)
    DEALLOCATE(fdirpar)
    DEALLOCATE(fdifpar)
    DEALLOCATE(st4)

    DEALLOCATE(tem11)
    DEALLOCATE(tem12)
    DEALLOCATE(tem13)
    DEALLOCATE(tem14)
    DEALLOCATE(tem15)
    DEALLOCATE(tem16)
    DEALLOCATE(radbuf)

  END IF





!
!-----------------------------------------------------------------------
!
!  Calculate analysis increments and write to file, if desired.
!
!-----------------------------------------------------------------------
!

  IF(incrdmp > 0) THEN
    CALL incrcalc(nx,ny,nz,nxndg,nyndg,nzndg,                           &
                  u,v,w,pprt,ptprt,qv,                                  &
                  qc,qr,qi,qs,qh,                                       &
                  uincr,vincr,wincr,pincr,ptincr,qvincr,                &
                  qcincr,qrincr,qiincr,qsincr,qhincr,                   &
                  istatus)

    CALL incrdump(nxndg,nyndg,nzndg,incrdmp,incdmpf,                    &
                  uincr,vincr,wincr,pincr,ptincr,qvincr,                &
                  qcincr,qrincr,qiincr,qsincr,qhincr,                   &
                  uincdmp,vincdmp,wincdmp,                              &
                  pincdmp,ptincdmp,qvincdmp,                            &
                  qcincdmp,qrincdmp,qiincdmp,qsincdmp,qhincdmp,         &
                  istatus)

  END IF
!
!-----------------------------------------------------------------------
!
!  Data dump of the model grid and base state arrays:
!
!-----------------------------------------------------------------------
!
  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        tem1(i,j,k)=0.
      END DO
    END DO
  END DO
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem2(i,j,k)=rhostr(i,j,k)/j3(i,j,k)
      END DO
    END DO
  END DO
!
  IF( hdmpfmt == 5 ) GO TO 700
  IF( hdmpfmt == 9 ) GO TO 700
!
!-----------------------------------------------------------------------
!
!  Find a unique name basdmpfn(1:lbasdmpf) for the grid and
!  base state array dump file
!
!-----------------------------------------------------------------------
!
  CALL gtbasfn(runname(1:lfnkey),dirname,ldirnam,hdmpfmt,               &
               mgrid,nestgrd,basdmpfn,lbasdmpf)

  WRITE(6,'(1x,a,a)')                                                   &
       'Data dump of grid and base state arrays into file ',            &
        basdmpfn(1:lbasdmpf)

  CALL setcornerll( nx,ny,x,y )   ! set the corner lat/lon

  grdbas = 1                ! Dump out grid and base state arrays only

  CALL dtadump(nx,ny,nz,nstyps,                                         &
               hdmpfmt,nchdmp,basdmpfn(1:lbasdmpf),                     &
               grdbas,filcmprs,                                         &
               u,v,w,ptprt,pprt,qv,                                     &
               qc,qr,qi,qs,qh,                                          &
               tem1,tem1,tem1,                                          &
               ubar,vbar,tem1,ptbar,pbar,tem2,qvbar,                    &
               x,y,z,zp,hterain,j1,j2,j3,                               &
               soiltyp,stypfrct,vegtyp,lai,roufns,veg,                  &
               tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,                &
               raing,rainc,prcrate,                                     &
               radfrc,radsw,rnflx,                                      &
               usflx,vsflx,ptsflx,qvsflx,                               &
               tem3,tem4,tem5)
!
!-----------------------------------------------------------------------
!
!  Find a unique name hdmpfn(1:ldmpf) for history dump data set
!  at time 'curtim'.
!
!-----------------------------------------------------------------------
!
  700 CONTINUE

  CALL gtdmpfn(runname(1:lfnkey),dirname,                               &
               ldirnam,curtim,hdmpfmt,                                  &
               mgrid,nestgrd, hdmpfn, ldmpf)

  WRITE(6,'(1x,a,a)') 'History data dump in file ',hdmpfn(1:ldmpf)


  grdbas = 0                ! No base state or grid array is dumped.

  CALL dtadump(nx,ny,nz,nstyps,                                         &
               hdmpfmt,nchdmp,hdmpfn(1:ldmpf),                          &
               grdbas,filcmprs,                                         &
               u,v,w,ptprt,pprt,qv,                                     &
               qc,qr,qi,qs,qh,                                          &
               tem1,tem1,tem1,                                          &
               ubar,vbar,tem1,ptbar,pbar,tem2,qvbar,                    &
               x,y,z,zp,hterain, j1,j2,j3,                              &
               soiltyp,stypfrct,vegtyp,lai,roufns,veg,                  &
               tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,                &
               raing,rainc,prcrate,                                     &
               radfrc,radsw,rnflx,                                      &
               usflx,vsflx,ptsflx,qvsflx,                               &
               tem3,tem4,tem5)

END PROGRAM arpsdas  


!

SUBROUTINE exsufx(fname,lenfnm,suffix,maxsuf,dotloc,lensuf) 2
  IMPLICIT NONE
  INTEGER :: lenfnm
  INTEGER :: maxsuf
  CHARACTER (LEN=1) :: fname(lenfnm)
  CHARACTER (LEN=1) :: suffix(maxsuf)
  INTEGER :: lensuf
  INTEGER :: dotloc
!
  INTEGER :: i
!
  DO i=lenfnm,1,-1
    IF(fname(i) == '.') GO TO 200
  END DO
  dotloc=0
  lensuf=0
  DO i=1,maxsuf
    suffix(i)=' '
  END DO
  RETURN
  200 CONTINUE
  dotloc=i
  lensuf=MIN(maxsuf,lenfnm-i)
  DO i=1,lensuf
    suffix(i)=fname(dotloc+i)
  END DO
  RETURN
END SUBROUTINE exsufx
!

SUBROUTINE exfroot(fname,lenfnm,froot,mxroot,lenroot) 1
  IMPLICIT NONE
  INTEGER :: lenfnm
  INTEGER :: mxroot
  CHARACTER (LEN=1) :: fname(lenfnm)
  CHARACTER (LEN=1) :: froot(mxroot)
  INTEGER :: lenroot
!
  INTEGER :: i
  INTEGER :: slashloc
  INTEGER :: dotloc
!
  DO i=lenfnm,1,-1
    IF(fname(i) == '/') EXIT
  END DO
!  101 CONTINUE
  slashloc=i
  DO i=slashloc,lenfnm
    IF(fname(i) == '.') EXIT
  END DO
!  201 CONTINUE
  dotloc=i
  lenroot=(dotloc-slashloc)-1
  DO i=1,lenroot
    froot(i)=fname(slashloc+i)
  END DO
  RETURN
END SUBROUTINE exfroot
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  SUBROUTINE SETCAT                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!

SUBROUTINE setcat(nx,ny,nz,ccatopt,zs,                                  & 1,2
           ptprt,pprt,qv,qc,qr,qi,qs,qh,                                &
           ptbar,pbar,                                                  &
           prcrate,icatg)
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Assign categories to grid locations according to the background
!  state to account for decorrelation across areas with active
!  parameterized convection and those without.
!
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  8/7/98
!
!  MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!  OUTPUT:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  INTEGER :: ccatopt
  REAL :: zs    (nx,ny,nz)  ! The physical height coordinate defined at
                            ! w-point of the staggered grid.
  REAL :: ptprt (nx,ny,nz)  ! Perturbation potential temperature (K)
  REAL :: pprt  (nx,ny,nz)  ! Perturbation pressure (Pascal)
  REAL :: qv    (nx,ny,nz)  ! Water vapor specific humidity (kg/kg)
  REAL :: qc    (nx,ny,nz)  ! Cloud water mixing ratio (kg/kg)
  REAL :: qr    (nx,ny,nz)  ! Rain water mixing ratio (kg/kg)
  REAL :: qi    (nx,ny,nz)  ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs    (nx,ny,nz)  ! Snow mixing ratio (kg/kg)
  REAL :: qh    (nx,ny,nz)  ! Hail mixing ratio (kg/kg)

  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
  REAL :: prcrate(nx,ny,4)     ! precipitation rate (kg/(m**2*s))
                               ! prcrate(1,1,1) = total precipitation rate
                               ! prcrate(1,1,2) = grid scale precip. rate
                               ! prcrate(1,1,3) = cumulus precip. rate
                               ! prcrate(1,1,4) = microphysics precip. rate
  INTEGER :: icatg(nx,ny)
!
  REAL :: prctrw,bltop,rhrkul,thevrkul
  PARAMETER (prctrw=1.0E-04,    & !  (kg/(m**2*s))
         bltop=1000.,   & !  m
         rhrkul=0.8,                                                    &
             thevrkul=9.0)  !  K**2
!
  INTEGER :: i,j
!  integer knt,k
  INTEGER :: knt1,knt2a,knt2b,knt3,knt4
  REAL :: pr,temp,qvsat,rh
!  real rhsum,rhmean
!  real the,thesum,thesq,themean,thevar
!
!-----------------------------------------------------------------------
!
!  Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_qvsat

!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat

!
!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
!
  INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
  knt1=0
  knt2a=0
  knt2b=0
  knt3=0
  knt4=0
!
!-----------------------------------------------------------------------
!
!  Initialize with default of no precip or precip influence
!
!-----------------------------------------------------------------------
!
  DO j=1,ny-1
    DO i=1,nx-1
      icatg(i,j)=1
    END DO
  END DO
!
  IF(ccatopt == 1) THEN
    DO j=1,ny-1
      DO i=1,nx-1
!
!-----------------------------------------------------------------------
!
!  Is convective precip turned on or heavy rain occuring?
!
!-----------------------------------------------------------------------
!
        IF (prcrate(i,j,3) > 0. .OR. prcrate(i,j,1) > prctrw) THEN
          knt4=knt4+1
          icatg(i,j)=4
!      print *, ' set i,j:',i,j,' cat 4',' prcrate(1):',
!    :             prcrate(i,j,1),' prcrate(3):',prcrate(i,j,3)
!
!-----------------------------------------------------------------------
!
!  Is it currently raining?
!
!-----------------------------------------------------------------------
!
        ELSE IF ( prcrate(i,j,1) > 0. ) THEN
          knt3=knt3+1
          icatg(i,j)=3
!      print *, ' set i,j:',i,j,' cat 3',' prcrate(1):',
!    :             prcrate(i,j,1),' prcrate(3):',prcrate(i,j,3)
!
!-----------------------------------------------------------------------
!
!  Evidence of rain-cooled air?
!  Lapse rate close to a moist adiabatic lapse rate in the lowest 2-km,
!  and a high relative humidity in that layer.  Or a high relative
!  humidity at the first level above ground.
!
!-----------------------------------------------------------------------
!
        ELSE
!
!-----------------------------------------------------------------------
!
!  For the layer 0-2km above first grid level (k=2).
!  Calculate the mean relative humidity.
!  Calculate the standard deviation of saturated equivalent
!  potential temperature.
!
!-----------------------------------------------------------------------
!
!      rhsum=0.
!      thesum=0.
!      thesq=0.
!      knt=0
!      DO 40 k=2,nz-1
!        IF(k.gt.2 .and. (zs(i,j,k)-zs(i,j,2)).gt.bltop) GO TO 41
!        knt=knt+1
!        pr=pprt(i,j,k)+pbar(i,j,k)
!        temp=(ptprt(i,j,k)+ptbar(i,j,k)) * (pr/p0)**rddcp
!        qvsat=f_qvsat( pr,temp )
!        rh=min(1.,(qv(i,j,k)/qvsat))
!        the=F_PT2PTE( pr,(ptprt(i,j,k)+ptbar(i,j,k)),qvsat)
!        rhsum=rhsum+rh
!        thesum=thesum+the
!        thesq=thesq+(the*the)
!  40     CONTINUE
!  41     CONTINUE
!      IF(knt.gt.0) THEN
!        rhmean=rhsum/float(knt)
!        themean=thesum/float(knt)
!      ELSE
!        rhmean=0.
!        themean=0.
!      END IF
!      IF(knt.gt.1) THEN
!        thevar=(thesq-((thesum*thesum)/float(knt)))/float(knt-1)
!      ELSE
!        thevar=9999.
!      END IF
!
!      IF (rhmean.gt.rhrkul .and. thevar.lt.thevrkul) THEN
!        knt2a=knt2a+1
!        icatg(i,j)=2
!        print *, ' set i,j:',i,j,' cat 2',' prcrate(1):',
!    :             prcrate(i,j,1),' prcrate(3):',prcrate(i,j,3)
!      ELSE
!
          pr=pprt(i,j,2)+pbar(i,j,2)
          temp=(ptprt(i,j,2)+ptbar(i,j,2)) * (pr/p0)**rddcp
          qvsat=f_qvsat( pr,temp )
          rh=qv(i,j,2)/qvsat
          IF(rh > 0.8) THEN
            knt2b=knt2b+1
            icatg(i,j)=2
!          print *, ' set i,j:',i,j,' cat 2*',' RH(2) = ',rh
          END IF

!      END IF
!      IF (mod(i,20).eq.0 .and. mod(j,10).eq.0)  THEN
!        print *, i,j,' rhmean=',rhmean,
!    :               ' thevar=',thevar,' icat=',icatg(i,j)
!        print *, ' '
!      END IF

        END IF
      END DO
    END DO
    knt1=((nx-1)*(ny-1))-knt2a-knt2b-knt3-knt4
    WRITE(6,'(a)') ' Precip correlation category assignment counts'
    WRITE(6,'(a,i6)') '     cat  1: ',knt1
    WRITE(6,'(a,i6)') '     cat 2a: ',knt2a
    WRITE(6,'(a,i6)') '     cat 2b: ',knt2b
    WRITE(6,'(a,i6)') '     cat  3: ',knt3
    WRITE(6,'(a,i6)') '     cat  4: ',knt4
  END IF
!
  CALL iedgfill(icatg,1,nx,1,nx-1, 1,ny,1,ny-1,1,1,1,1)
!
  RETURN
END SUBROUTINE setcat