    module driver
!
!  SYNOPSIS
!    use driver
!
!  DESCRIPTION
!  This file defines the main program and "driver" subroutines to run
!  the data assimilation code.
!
!  COMMAND-LINE OPTIONS SUPPORTED IN THIS VERSION
!  -h         : print explanation of command line and exit
!  -g         : indicates a global field; useful for debugging
!  -b  file   : file containing background state (default: 'egrid')
!  -a  file   : file containing analysis vectors (default: 'egrid.out')
!  -o  file   : file containing the observations (default: 'yobs')
!  -d  file   : file to which we print diagnostics
!  -m  file   : file containing analysis ensemble mean (default: none)
!  -n  nsol  : number of ensemble solutions (default 10)
!  -N  file   : reads namelist input from file.  Values in this file
!               overwrite those set by the command line arguments that
!               appear to the left of -N and are overwritten by command line
!               arguments that appear to the right of -N.
!  -c  check  : type of consistency check (if any) to apply to the observations.
!               CHECK may be one of "whitaker" or "obs"/"observation" or "none"
!               according as Whitaker's criterion or an observation-based
!               criterion is to be used.  The default is "none".
!  -e  errthresh : if the data assimilation cannot be completed on
!               more than ERRTHRESH of the patches, then judge the run
!               a failure.
!  -q  id     : sequence id (if any; usually the date in yyyymmddhh format)
!  -t  truegrid_file : file containing the "true" atmospheric state vector
!  -v  lev    : verbosity level: 0 (quiet), 1 (default: basic info)
!               2 (shows progress), 3 (chatty)
!  -W         : apply Whitaker's internal data consistency check
!  -x  xwidth : number of longitude points ON EACH SIDE of center (default 3)
!  -y  ywidth : number of latitude points ON EACH SIDE of center (default 3)
!               The default yields a 7x7 patch (3*2+1 points per side)
!  -z  file   : file containing vertical levels as function of latitude
!
!  REQUIRED PACKAGES
!  io, precision, ferror, command_line_support - in support directory
!  brian, diagnostics, grid_model, grid_obs, grid_io - in current directory
!
!  ROUTINES IMPLEMENTED
!
!   subroutine readin(partial_grid)
!   subroutine mainloop(retcode)
!   subroutine update_errcount(diag,errcount,ok)
!   subroutine parseline(partial_grid)
!   subroutine print_usage
!   subroutine initialize(partial)
!   program data_assim
!
!  MODIFICATION HISTORY
!  09/27/05 - Added facilities to vary the vertical coordinate of observations
!             based on the surface pressure of individual ensemble solutions.
!  08/08/05 - Added facilities for verbose dumps of background and obs in
!             selected patches.
!  07/11/05 - Made all internal computations in energy units; garbage values
!             are now regarded as a non-fatal error (we return the background
!             instead).
!  06/22/05 - Added enhanced variance inflation (brian.f90 and covar_setup.f90).
!  06/11/05 - Returns an error status if any garbage values are generated
!             along the way, irrespective of any other error threshold.
!  03/29/05 - Modified namelist to accept observation consistency check option
!  01/27/05 - Bug fixes to 01/21/05 version; moved main program to a separate 
!             source file; added -T option so temperature observations can 
!             be specified in either Celsius or Kelvin
!  01/21/05 - Changes to handle real observations.
!  12/14/04 - bug fixes and changes to grid_model for more efficient array
!             handling in pre- and postanalysis.  (Cache-unfriendly parts are
!             relegated here for overall cluster speed.)
!  10/11/04 - Changes in called routines to handle "real" data.
!  09/29/04 - modified internal arrays to have NSOL as the first dimension.
!             Only P, T, U, V are assimilated.
!  09/18/04 - Removed original assimilation scheme and replaced with Brian's.
!             Return matrix E that gives optimal linear combination of
!             the perturbations at the center for eventual use in handling
!             unassimilated variables and surface boundary conditions.
!  08/01/04 - Added code to assim.f90 to determine the extent to which the
!             analysis ensemble agrees with the computed analysis covariance
!             matrix.  Also added -B flag to omit Brian's method.
!  07/28/04 - Small changes to logfile printouts.
!  06/26/04 - Bug fixes and optimization for Brian's method to compute only
!             the center point of the patch.
!  06/17/04 - Changes to support both Brian's method and the original one,
!             including diagnostics for both.
!  06/14/04 - Added option to read namelist input.
!  06/08/04 - New version to implement Brian's assimilation scheme.
!             Incorporates some previous changes to precompute the ensemble
!             means.
!  03/19/04 - bug fix to correct patch initialization.
!  12/05/03 - added option to vary random perturbations, covariance inflation
!             and method by latitude and vertical level.
!  11/10/03 - added option for random pertubations at surface and to use
!             enhanced variance inflation at surface and simple inflation
!             elsewhere.
!  11/02/03 - added extra diagnostics to detect NaNs and Infs on input.
!  10/31/03 - version to implement vertical localization.  Also added
!             -e errthresh flag to set failure threshold.
!  08/31/03 - accommodate assimilation by latitudinal strips instead of
!             longitudinal ones.  Removed support for averaging.
!             Added support for centering patches at the poles.
!  08/12/03 - added measures of ensemble spread to diagnostics package
!             (extra.f90).  We compute and print the standard deviation of
!             selected atmospheric vector components at the center of each 
!             patch.
!  07/28/03 - modified solution 3 to use Eq. (66) in revised MWR paper instead
!             of Eqs. (79)-(80).  This version appears to eliminate the
!             catastrophic cancelation that occurs with the other formulation.
!  07/23/03 - added option for conventional variance inflation instead of
!             the enhanced scheme (which remains the default)
!  05/19/03 - revised diagnostics for case of missing observations
!  05/15/03 - added extra diagnostics
!  03/03/03 - more modifications to diagnostics output
!  02/25/03 - added -d option and extra diagnostics on patches
!  02/08/03 - miscellaneous enhancements and bug fixes.
!  01/28/03 - "parallelized" version that can assimilate in longitudinal
!             strips on independent processors.  Includes new "grid manager"
!             interface for handling data transmission to and from the
!             model grid.
!  12/09/02 - bug fixes: ensemble size set to k+1 if kprime is not specified;
!             prints xwin, ywin and patch size on log file.
!  11/21/02 - first production version by Eric; incorporated Solution 3
!             and averaging
!-----------------------------------------------------------------------------
!  DRIVER - supporting driver routines.
!  For convenience, we maintain the list of input and output files here.
!  These names are used and manipulated only by the routines in this
!  file.
!
    use assimilation
    use diagnostics
    use command_line_support
    implicit none
!
!  Counts of reasons that a given assimilation might fail.
!  See assim_vars and assim for more details.
!
    character(*),parameter::VERSION='09/27/05: varying z in observations'
!
    type error_reason  ! see assim_vars
       integer::attempted     ! number of times the method was attempted
       integer::no_obs=0      ! no observations in patch
       integer::symsqrt_fail=0   ! can't compute symmetric sqrt
       integer::no_pa_ev=0       ! can't do eigen decomp of P^a
       integer::neg_cov_ev=0  ! P^a has negative eigenvalues
       integer::pa_sing=0     ! P^a is singular
       integer::enh_fail=0    ! enhanced variance inflation failed
       integer::garbage=0     ! method generates garbage values in analysis
       integer::brian=0       ! Brian's method failed
       integer::other=0       ! any other reason
    end type error_reason
!
!  File containing the background model grid
!
    character(MAX_FILE_NAME_LEN),save,private::bfile='egrid'
!
!  File to be created that contains the analysis model grid
!
    character(MAX_FILE_NAME_LEN),save,private::afile='egrid.out'
!
!  File containing the observations
!
    character(MAX_FILE_NAME_LEN),save,private::ofile='yobs'
!
!  File containing the vertical levels as function of latitude.
!  If no such file is present (the default), then we don't do
!  vertical localization or covariance inflation.
!
    character(MAX_FILE_NAME_LEN),save,private::zfile=' '
!
!  File (if any) to be created that contains the analysis mean
!
    character(MAX_FILE_NAME_LEN),save,private::mfile=' '
!
!  File (if any) to be created that contains the diagnostics
!
    character(MAX_FILE_NAME_LEN),save,private::dfile='diag'
!
!  File (if any) containing the "true" nature grid
!
    character(MAX_FILE_NAME_LEN),save,private::tfile=' '
!
!  Debugging options: process just one or a few model grid points.
!  idxgrid(:,1) contains longitudes, (:,2) latitudes, (:,3) levels.
!  Normally, (:,2) is set in mainloop and the remaining values are
!  unchanged.  For debugging, however, values can be set via
!  NAMELIST using the -N option; see PARSELINE, below.
!
    integer,save,private::idxgrid(2,3)
    data idxgrid/1,LONF,0,0,1,LEVS/
!
!  If a fraction greater than ERRTHRESH of the total patches has
!  a numerical or other problem such that the assimilation cannot
!  be completed, then judge the run a failure.
!
    real,save,private::errthresh=0.0
!
    contains
!-----------------------------------------------------------------------------
    subroutine readin(partial_grid)
!  READIN - read in background grid and observations; allocate storage (once
!  per program invocation).  Exit if the indicated files cannot be opened
!  for writing.
!
    logical,intent(in)::partial_grid  ! true if assimilating partial global grid
!
!  Local variables
!
    integer::ix,iy,south,north,kept,nobs,nsol
    real::t0,t1
    character(4),save::gridhdr='grid'  ! defaults for latitude strips
    character(4),save::obshdr='obs '
    character(4)::header='none'
    character(NUM_OBS_TYPES)::obslist
    character::tunit
!
!  Announcements
!
    call cpu_time(t0)
    write(LOGFILE,100) "assim_p: Brian's method: version",version
    write(LOGFILE,100) 'background data file:',trim(bfile)
    if(partial_grid) then
       write(LOGFILE,100) 'partial grids'
    else
       write(LOGFILE,100) 'global grids '
    endif
    write(LOGFILE,100) 'observation data file:',trim(ofile)
    call get_observation_types(temp_units=tunit)
    write(LOGFILE,100) 'observation temperature units:',tunit
    if(len_trim(tfile).gt.0) &
       write(LOGFILE,100) 'nature run grid file:',trim(tfile)
100 format(a,1x,a)
!
!  First check that the output files are writable.  No point in doing
!  a lot of work if they aren't!
!
    if(partial_grid) header='true'  ! need slice of true grid with 'true' header
    call check_writable(afile)
    write(LOGFILE,100) 'analysis file:',trim(afile)
    if(len_trim(mfile).gt.0) then
       call check_writable(mfile)
       write(LOGFILE,100) 'analysis mean file:',trim(mfile)
    endif
!
!  Vertical localization input file - this must be read before observation
!  quality control is done so that we apply any requested background
!  covariance inflation.
!
    if(len_trim(zfile).gt.0) then
       write(LOGFILE,101) trim(zfile)
101    format('file defining vertical localization: ',a)
       call vertical_localization_setup(zfile)
    else
       write(LOGFILE,102)
102    format('no vertical localization')
    endif
!
!  Read in the background grid and the observations
!
    if(.not.partial_grid) then    ! whole globe file has no header
       gridhdr='none'
       obshdr='none'
    endif
    nsol=ensemble_size()
    write(LOGFILE,103) 'ensemble size:',nsol
103 format(a,1x,2i8)
    call model_grid_init(bfile,afile,mfile,dfile,gridhdr)
    call obs_init(ofile,obshdr,nobs,nsol)
    call precompute_ensemble_differences
    call model_grid_boundary(south,north)
    write(LOGFILE,104) 'model grid boundaries:',south,RLATS(south),north, &
       RLATS(north)
104 format(a,'  S: ',i3,' (',f7.3,' deg) N: ',i3,' (',f7.3,' deg)')
    write(LOGFILE,105) seqid()
105 format('input grid sequence id: ',i0)
!
!  Parse the diagnostics file and initialize the "true" background
!  grid, if applicable.
!
    if(len_trim(dfile).gt.0) then
       call diag_init(header,dfile,tfile,bfile,seqid(),ofile,afile)
       write(LOGFILE,100) 'diagnostics file:',trim(dfile)
    endif
!
!  Miscellaneous details
!
    write(LOGFILE,106) errthresh
106 format('max tolerable fraction of patches with assimilation errors: ',f6.3)
    call patch_radius(lon=ix,lat=iy)
    write(LOGFILE,101) 'patch width (long) x height (lat):',2*ix+1,2*iy+1
    call assimilation_boundary(south,north)
    write(LOGFILE,104) 'assimilation latitude boundaries:', &
       south,RLATS(south),north,RLATS(north)
    call get_observation_types(obslist)
    write(LOGFILE,107) obslist
107 format('observation types assimilated: ',a)
!
    call cpu_time(t1)
    write(LOGFILE,120) 'initialization time (sec):',t1-t0
120 format(a,f10.2)
    return
!
    contains
       subroutine check_writable(filename)
       character(*),intent(in)::filename
       character(8)::writable
!
       open(unit=OUTNOSWAP,file=filename,form='unformatted')
       inquire(unit=OUTNOSWAP,write=writable)
       if(writable.eq.'NO') call fexit(trim(filename)//': not writable',2)
       close(OUTNOSWAP)
       return
       end subroutine check_writable
    end subroutine readin
!---------------------------------------------------------------------------
    subroutine mainloop(retcode)
! MAINLOOP - do one data assimilation step by looping through all
! longitude,latitude pairs.
!
!  Required packages:
!
!  Argument:
!  RETCODE := return code; 0 for no errors, nonzero otherwise.
!
    integer,intent(out)::retcode
!
!  Local variables.  The arrays are allocatable to avoid stack overflow.
!
    integer::ix,iy,iz,asouth,anorth
    integer::npatch,nfail,nsol,completed
    type(patch_index)::loc
    real(D),allocatable::w(:,:),sv(:)
    real::t0,t1,tcpu
    logical::have_analysis,exceeded_errthresh
    type(error_reason)::errcount
    type(patch_diag)::diag
!
!  Initialize.  ASOUTH and ANORTH give the boundaries of the latitudinal strip
!  over which we center grids for the analysis.
!
    call cpu_time(t0)
    nsol=ensemble_size()
    allocate(w(nsol,nsol))
    call assimilation_boundary(asouth,anorth)
!
!  Main processing loop.  Proceed in latitude-major order.
!  The number of rows in XHAT and AMEAN always equals the number of
!  dynamical variables in the patch.  NFAIL counts the number of assimilation
!  errors for reasons other than a lack of observations in the patch.
!  We quit if assimilation errors occur on more than ERRTHRESH of the
!  patches after we're a quarter of the way through our assigned grid points.
!  IDXGRID lets us run only a single patch if we choose (useful for
!  debugging) without changing any code.
!
    if(any(idxgrid(:,2).eq.0)) idxgrid(:,2)=(/asouth,anorth/)
    npatch=product(idxgrid(2,:)-idxgrid(1,:)+1)
    !!allocate(diag%iceb_sing(nsol))  ! debug: write out singular values
    nfail=0
    completed=0   ! total number of completed patches
    do ix=idxgrid(1,1),idxgrid(2,1)  ! normally 1,LONF
       do iy=idxgrid(1,2),idxgrid(2,2)  ! normally asouth,anorth
          do iz=idxgrid(1,3),idxgrid(2,3)  ! normally 1,LEVS
             loc=patch_index(ix,iy,iz)
             call assimilate(loc,nsol,w,diag)
             call output_patch_diag(loc,diag)
             call update_errcount(diag,errcount,have_analysis)
             completed=completed+1
             if(.not.have_analysis) nfail=nfail+1
             exceeded_errthresh=completed.ge.npatch/4.and.&
                real(nfail)/real(npatch).gt.errthresh
             if(exceeded_errthresh) goto 10
          enddo
       enddo
       if(verbosity().gt.VERB_TERSE) then
          call cpu_time(tcpu)
          write(IERR,130) ix,tcpu,nfail
       endif
    enddo
130 format('mainloop: finished longitude',i4,' time:',f9.2, &
       ' failures so far:',i7)
!
!  All done.
!
 10 continue
!
!  Always report the total CPU time for the assimilation step.
!
    call cpu_time(t1)
    write(LOGFILE,135) 'assimilation processing time',t1-t0
135 format(a,' (sec):',f9.2)
    t0=t1  ! reset the clock
!
!  If we exceeded the error threshold then no analysis grid is generated,
!  as we've terminated the analysis early.
!
    if(exceeded_errthresh) then
       if(verbosity().gt.VERB_TERSE) call finfo('error threshold exceeded')
       write(LOGFILE,136) nfail,npatch,real(nfail)/real(npatch)
136    format('assimilation failed on ',i0,' of ',i0,' patches; ', &
          'fraction: ',f7.4)
       retcode=1
       return
    endif
!
!  Output the analysis grid and optionally the analysis mean, along with
!  observation diagnostics.
!
    call report("Statistics on Brian's method",errcount)
    call model_grid_fini
    call summary_diagnostics(nsol,len_trim(tfile).gt.0)
    if(verbosity().gt.VERB_SILENT) then
       call cpu_time(t1)
       write(LOGFILE,135) 'output time',t1-t0
    endif
!
!  If the method has generated garbage values along the way (as opposed
!  to other non-fatal problems), then flag an error return.
!  For debugging purposes, we still want to have output the analysis grid,
!  which is why this check is postponed till last.
!
    if(errcount%garbage.gt.0) then
       call finfo('Garbage values in analysis grid')
       retcode=1
    else
       retcode=0
    endif
    return
    contains
!-----------
       subroutine report(label,e)
       character(*),intent(in)::label
       type(error_reason),intent(in)::e
!
       write(LOGFILE,140) label
140    format(/a)
       write(LOGFILE,141) e%attempted,'total attempts'
       write(LOGFILE,141) e%no_obs,'patches with no observations'
       write(LOGFILE,141) e%symsqrt_fail,'symmetric sqrt failed'
       write(LOGFILE,141) e%no_pa_ev,'eigenval decomp of P^a failed'
       write(LOGFILE,141) e%neg_cov_ev,'P^a had negative eigenvalues'
       write(LOGFILE,141) e%pa_sing,'P^a was singular'
       write(LOGFILE,141) e%enh_fail,'enhanced inflation failed'
       write(LOGFILE,141) e%garbage,'method generated garbage values'
       write(LOGFILE,141) e%other,'other reason'
       write(LOGFILE,141) e%brian,'Brian''s method failed; original substituted'
141    format(i8,1x,a)
       return
       end subroutine report
!-----------
    end subroutine mainloop
!---------------------------
    subroutine update_errcount(diag,errcount,ok)
!  UPDATE_ERRCOUNT - update the error count for each assimilation method.
!  Arguments:
!   DIAG : structure containing patch diags and return code from assim
!   ERRCOUNT :=: updated error count.
!   OK := set to true if we were still able to update
!        the patch when observations are available.  Also set to true if
!        there are no observations, since we simply
!        return the background as the maximum-likelihood analysis.
!
    type(patch_diag),intent(in)::diag
    type(error_reason),intent(inout)::errcount
    logical,intent(out)::ok
!
    errcount%attempted=errcount%attempted+1
    select case(diag%info)
    case(0)   ! success or not called - no error
       ok=.true.
       return
    case(NO_OBS_UPDATE)
       ok=.true.
       errcount%no_obs=errcount%no_obs+1
    case(INFLATION_FAIL)
       ok=.true.
       errcount%enh_fail=errcount%enh_fail+1
    case(SYM_SQRT_FAIL)
       ok=.false.
       errcount%symsqrt_fail=errcount%symsqrt_fail+1
    case(PA_NO_EV)
       ok=.false.
       errcount%no_pa_ev=errcount%no_pa_ev+1
    case(P_NEG_EV)
       ok=.false.
       errcount%neg_cov_ev=errcount%neg_cov_ev+1
    case(PA_SINGULAR)
       ok=.true.
       errcount%pa_sing=errcount%pa_sing+1
    case(GARBAGE_VALUES)
       ok=.false.
       errcount%garbage=errcount%garbage+1
    case default
       ok=.false.
       errcount%other=errcount%other+1
    endselect
    if(diag%brian_failed) errcount%brian=errcount%brian+1
    return
    end subroutine update_errcount
!---------------------------
    subroutine parseline(partial_grid)
!  PARSELINE - parse execute line.  Getopt works much like the C version
!  (man 3c getopt).
!
!  Arguments:
!  PATCHW, PATCHH    := width and height, respectively, of the grid patches
!  PARTIAL_GRID := set when the -p option is specified; we look for 
!                  a partial grid file with appropriate header
!  Side effects:
!  - may alter the filenames in DRIVER
!  - sets many write-once global variables in grid_manager
!
    logical,intent(out)::partial_grid
!
!  Local variables
!
    character(50),save::options='A:b:c:a:d:e:ghm:n:N:o:q:T:t:v:x:y:z:'
!
    character(MAX_FILE_NAME_LEN)::optarg,nmlfile
    integer::iopt,optind,iseq,nsol,verblev
    integer::lonrad,latrad
    character::opt,obs_temp_units='K'
    character(4)::obstypes='PTUV'
    character(12)::obs_consistency_check='none'
!
!  Optional namelist input can affect most variables.  For convenience in
!  specifying the size of the ensemble, we use a local variable called NSOL
!  and overwrite the value of the corresponding variable in assim_vars
!  as needed.
!
    namelist/nlvars/ bfile,afile,ofile,tfile,dfile,mfile,zfile,idxgrid, &
       verblev,nsol,lonrad,latrad,errthresh,partial_grid,obstypes,iseq, &
       obs_temp_units,obs_consistency_check,dump_loc
!
    call patch_radius(lon=lonrad,lat=latrad)  ! get default value
    partial_grid=.true.  ! assume partial grid by default
    do
       call getopt(iopt,optarg,optind,options)
       if(iopt.lt.0) exit
       opt=char(iopt)
       select case(opt)  
       case('A')  ! decide which observations we'll assimilate (PTUV)
          obstypes=optarg
       case('a')
          afile=optarg
       case('b')
          bfile=optarg
       case('c')
          obs_consistency_check=optarg
       case('d')
          dfile=optarg
       case('e')
          read(optarg,*) errthresh
       case('h')
          call print_usage
          call program_exit(1)
       case('g')
          partial_grid=.false.
       case('m')
          mfile=optarg
       case('n')
          read(optarg,*) nsol
       case('N')
          nmlfile=optarg
          open(IN,file=nmlfile,status='old')
          read(IN,nlvars)
          close(IN)
          call set_observation_types(temp_units=obs_temp_units)
       case('o')
          ofile=optarg
       case('q')
          read(optarg,*) iseq
          call set_sequence_id(iseq)
       case('t')
          tfile=optarg
       case('T')  ! arg must be 'C' or 'K' for Celsius or Kelvin
          call set_observation_types(temp_units=optarg)
       case('v')
          read(optarg,*) verblev
       case('x')
          read(optarg,*) lonrad
       case('y')
          read(optarg,*) latrad
       case('z')
          zfile=optarg
       case default
          write(IERR,110) opt
          call program_exit(1)
       endselect
    enddo
!
    call set_verbosity(verblev)
    call set_observation_types(obstypes)
    call set_ensemble_size(nsol)
    call set_patch_radius(lonrad,latrad)
    call apply_consistency_check(obs_consistency_check)
!
110 format('unknown option ',a)
!
!  Sanity checks on command-line data
!
    if(errthresh.gt.1.0.or.errthresh.lt.0.0) &
       call fexit('errthresh must be in [0,1]',2)
    return
    end subroutine parseline
!-----------------------------------------------------------------------------
    subroutine print_usage
!  PRINT_USAGE - print a brief synopsis of the command line on the standard
!  error stream.
!
    character(MAX_FILE_NAME_LEN)::arg
!
    call getarg(0,arg)
    write(IERR,100) trim(arg),version
100 format(a,': Version ',a,': Options are:')
    write(IERR,110) '-A assim_variable_list (subset of PTUV)'
    write(IERR,110) '-a output_analysis_file'
110 format(a)
    write(IERR,110) '-b background_ensemble_file'
    write(IERR,110) '-d diagnostics_output_file'
    write(IERR,110) '-e error_fraction_threshold'
    write(IERR,110) '-m analysis_mean_file'
    write(IERR,110) '-n number_of_bred_vectors (kprime+1)'
    write(IERR,110) '-N namelist_input_file'
    write(IERR,110) '-o observations_file'
    write(IERR,110) '-q seqid (require specified sequence id for input grid)'
    write(IERR,110) '-T obs_temp_units (C or K); default is K=Kelvin)'
    write(IERR,110) '-t true_grid_file'
    write(IERR,110) '-v verbosity level (0:quiet), (1:default),'
    write(IERR,110) '   (2:show progress), (3:chatty)'
    write(IERR,110) '-W (apply Whitaker''s data consistency check)'
    write(IERR,110) '-x width (of longitude patch)'
    write(IERR,110) '-y width (of latitude patch)'
    write(IERR,110) '   each patch is a (2x+1)*(2y+1) rectangle of grid points'
    write(IERR,110) '-z zrad_dat_file'
    return
    end subroutine print_usage
!--------------------------------------------------------------------------
    end module driver
