! POSTANALYSIS - combine one or a family of output files from the
! assim or assim_p programs into a collection of sigma files, one
! per ensemble solution.
!  Modification history:
!  08/25/05 - Added compatibility to both 2004 and 2001 GFS versions.
!  06/17/05 - Added capability to convert physical perturbation files
!             directly to sigma file format.
!  06/16/05 - Changes to accommodate overhauled spectral code.
!  05/25/05 - Changes to handle sigma file format for 2004 GFS; added option
!             for namelist input (for debugging); simplified interface:
!             -f file now specifies a sigma file from which all unassimilated
!             variables are taken.
!  03/07/05 - Uses new spectral transform routines to fix N/S hemisphere
!             indexing.
!  12/14/04 - altered array indexing of input model grid format for more
!             cache-friendly processing.  Also put burden of rescaling
!             onto the data assimilation code.
!  10/12/04 - changes to read the new output model grid format.
!  08/25/04 - made compatible with revised grid manager code and
!             eliminated support for longitude-major formats (except
!             for the debug output model grid).
!  07/27/04 - added handling of surface analysis files from GFS code.
!  05/17/04 - added -q and -z options for selective copying of humidity and
!             ozone.
!  02/20/04 - Fixed the handling of humidity and ozone
!             when a sigma file for the analysis mean is requested.
!             When a mean sigma file is requested, these
!             values are averaged from those in each of the reference sigma
!             files.  (Note: divergence and vorticity are computed from
!             the analyses.)
!  10/17/03 - version by Eric that can input grid formats in either
!             longitude- or latitude-major format
!  08/28/03 - version by Eric to handle new grid format that slices
!             latitudinally rather than longitudinally.
!  05/30/03 - version by Eric to handle families of grid files.
!  10/04/02 - initial version by Istvan Szunyogh.
!-----------------------------------------
!  Implementation note on the Karman/Keckwulf clusters:
!  The input and output sigma files need to be in big-endian format,
!  as this is the format that the GFS code uses.  For correct operation,
!  the environment variable FORT90L must be set as follows:
!     export FORT90L="-Wl,-T1,-T2"
!  This is a Lahey-specific environment variable.
!-----------------------------------------
! COMMAND LINE OPTIONS:
!  postanalysis -f name [options] [input_grid_file]
!  -f name : (no default) This is a required option.  NAME specifies
!            a sigma file containing scaling coefficients and all unassimilated
!            data.
!            Two cases:
!            (1) If name is a template name, then name refers to a family
!            of sigma files, one for each ensemble member, from which
!            unassimilated values are copied.
!            (2) If name is not a template name, then unassimilated variables
!            are copied from this file into EACH of the output sigma files,
!            one per ensemble member (as well as to the mean sigma file
!            if -m is specified).  In this case, name should refer
!            to the appropriate sigma file from the nature run or NCEP analysis.
!  -d date : The date, in yyyymmddhh format, used to document the output
!            sigma files.  (Otherwise y=m=d=h=0 is used.)
!  -v lev  : (default: lev=0) Set the verbosity level to LEV.
!            0 gives silent operation (unless an error occurs).
!            1 gives brief status, 2 shows many details.
!  -m mean : (default: none) The name of the sigma file containing the
!            spectral transform of the analysis mean.
!  -n nsol : (default: nsol=1) number of ensemble members.
!  -p np   : (default: np=1) number of input grid files to expect.
!            The value for NP _must_ agree with that used in the
!            corresponding invocation of preanalysis; otherwise, some grid
!            points may have garbage values.
!  -s name : (default: name="sigma.out") Base name for the output sigma
!            files.  If NSOL=1, then name is the actual name of the
!            output file.  If NSOL > 1, then name is a template name.
!            NSOL output files are created, named name1, ..., nameNSOL.
!            If name contains the substring %n, then %n is replaced with
!            the numbers 1 ... NSOL to create the sequence of output
!            sigma files.  Example: -s pnew%n.00 would yield files named
!            pnew1.00, pnew2.00, etc.
!  -S name : (no default)  If present, name must be a template name that
!            specifies a family of surface analysis files, one for each
!            ensemble solution, that are output from the GFS model.
!            These are averaged together to create a mean file that is
!            indicated with %n replaced by 0.  If omitted, the program does
!            not process any surface analysis files.
!  -h      : print a summary of the command-line options and halt.
!  input_grid_file :  (default "egrid.out")  The name of the file from
!            the data assimilation step.  If -p is specified, then there
!            are NP slices to the data, and input_grid_file should specify
!            a template file.  Substitutions in the name are performed as
!            described in the -s option above.
!
    module postanalysis
    use spectral
    use fsplit
    use surface
    use command_line_support
    use gfs_ids
    implicit none
!
    character(*),parameter::VERSION='08/25/05 (2001/04 models, energy units)'
!
    character(MAX_FILE_NAME_LEN),save::phys_grid_fn='egrid.out'
!
!  Output file, which is a template name if nsol > 1
!
    character(MAX_FILE_NAME_LEN),save::sigma_out_fn='sigma.out'
!
!  Optional output file containing the analysis mean in spectral coords
!
    character(MAX_FILE_NAME_LEN),save::mean_sigma_fn=' '
!
!  Required sigma reference file for spectral data and any unassimilated
!  variables (-f flag).
!
    character(MAX_FILE_NAME_LEN),save::ref_sig_fn=' '  ! no default
!
!  Optional surface analysis file family (-S flag).
!
    character(MAX_FILE_NAME_LEN),save::sfc_fn=' '
!
    integer,parameter::HEADER_TAG_LEN=4
    character(HEADER_TAG_LEN),save::expected_tag='aout'
!
!  Number of partitions of output grid files from data assimilation
!
    integer,save::nsplit=1
!
!  By default, the input consists of latitude strips processed by the
!  data assimilation code.  But if -P is specified, then the input
!  consists of physical grid perturbation files that will be transformed
!  directly to sigma.
!
    logical,save::input_is_pert=.false.
    contains
!------------------------------------------------------------------------------
    subroutine do_everything
!  DO_EVERYTHING - so we don't have to go to the end of the file to
!  see the main program
!
    character(10)::cdate=' '   ! in yyyymmddhh format
    type(sfcanl)::avg
    character(MAX_FILE_NAME_LEN)::sfc_out
    real(R)::fhour
    integer::nsol
!
    if(verbosity().gt.VERB_TERSE) write(IOUT,100) version
100 format('postanalysis: version ',a)
    call parseline(cdate,nsol)
    call convert_to_sigma(cdate,nsol,ref_sig_fn)
!
!  Surface average is member 0
!
    if(len_trim(sfc_fn).gt.0) then
       call read_surface_files(sfc_fn,avg,fhour,nsol)
       call numeric_filename(sfc_fn,'%n',0,sfc_out)
       call write_sfcanl(sfc_out,avg,0.0,cdate) ! new initial condition
    endif
    return
    end subroutine do_everything
!-----------------------------------------------------------------------------
    subroutine convert_to_sigma(cdate,nsol,ufn)
!  CONVERT_TO_SIGMA - convert from physical to spectral coordinates and
!  write out a family of sigma files, one for each ensemble member, suitable
!  for input to the GFS model to run a forecast.
!  CDATE is a time stamp of the form yyyymmddhh that is written into
!  each of the sigma files.
!  UFN : the name of a family of sigma files from which to copy
!     the coefficients of any unassimilated variables.  This name is a 
!     template name if it contains "%n".  Otherwise, ALL output files
!     will contain the SAME values of any unassimilated variables.
!
    character(*),intent(in)::cdate,ufn
    integer,intent(in)::nsol
!
!  Local variables.
!
    type(sigma_file_header)::header
    type(spectral_coeffs)::spect    ! coeffs for assimilated variables
    type(spectral_coeffs)::unspect  ! coeffs for unassimilated variables
    type(spectral_coeffs)::avgspect  ! coeffs for averaged variables
    real(GRIDKIND),allocatable::egrid(:,:,:,:,:),mean(:,:,:,:)
    integer::nc,ns,verblev
    character(MAX_FILE_NAME_LEN)::sfn,refname,fn
    logical::template_ref_name
!
!  Initialize.  The reference sigma file has the same endian convention
!  as the output sigma files from this program: it may require byte-swapping.
!  We regard the reference sigma file as a template name if it contains
!  the substring "%n".  Copy all unassimilated variables from it.
!
    template_ref_name=index(ufn,'%n').gt.0
    if(template_ref_name) then  ! use the first one for initialization
       call numeric_filename(ufn,'%n',1,refname)
    else
       refname=ufn
    endif
    call spectral_code_init(refname,header,unspect)
!
    verblev=verbosity()
    if(verblev.gt.VERB_SILENT) &
       write(LOGFILE,'(2a)') 'postanalysis: output family: ',trim(sigma_out_fn)
!
!  Read the model grids in physical coordinates.
!
    allocate(egrid(LONF,NAVARS,LATG,LEVS,nsol))
    ns=nsplit
    if(len_trim(mean_sigma_fn).gt.0) then
       allocate(mean(LONF,NAVARS,LATG,LEVS))
       avgspect%gz=0.0
       avgspect%rq=0.0
       avgspect%oz=0.0
    endif
!
    if(input_is_pert) then
       if(allocated(mean)) mean=0.0
       do nc=1,nsol
          if(nsol.eq.1) then
             fn=phys_grid_fn
          else
             call numeric_filename(phys_grid_fn,'%n',nc,fn)
          endif
          call read_perturbation_file(fn,egrid(1,1,1,1,nc))
          if(allocated(mean)) mean=mean+egrid(:,:,:,:,nc)
       enddo
       if(allocated(mean)) mean=mean/nsol
    else if(allocated(mean)) then
       call read_grid_files(phys_grid_fn,egrid,nsol,ns,mean)
    else
       call read_grid_files(phys_grid_fn,egrid,nsol,ns)
    endif
! 
!  Convert assimilated variables to spectral coordinates,
!  merge unassimilated variables from reference sigma files,
!  and write out the merged coefficients.  Compute the mean of the
!  spectral coefficients as we go along if the mean field is requested.
!
    do nc=1,nsol 
       call from_physical_space(egrid(1,1,1,1,nc),spect)
       call numeric_filename(sigma_out_fn,'%n',nc,sfn)
       if(verblev.gt.VERB_TERSE) &
          write(LOGFILE,'(2a)') 'postanalysis: writing ',trim(sfn)
       if(template_ref_name.and.nc.gt.1) then  ! 1st file was read previously
          call numeric_filename(ref_sig_fn,'%n',nc,refname)
          call read_sigma_file(refname,header,unspect)
       endif
       call writeout(header,unspect,spect,cdate,sfn)
       if(allocated(mean)) then  ! average unassimilated humidity and ozone
          avgspect%gz=avgspect%gz+spect%gz
          avgspect%rq=avgspect%rq+spect%rq
          avgspect%oz=avgspect%oz+spect%oz
       endif
    enddo
!
!  Transform and output the analysis mean as needed.  
!
    if(allocated(mean)) then
       if(verblev.gt.0) &
          write(LOGFILE,'(2a)') 'postanalysis: writing ',trim(mean_sigma_fn)
       call from_physical_space(mean,avgspect)
       unspect%gz=avgspect%gz/nsol
       unspect%rq=avgspect%rq/nsol
       unspect%oz=avgspect%oz/nsol
       call writeout(header,unspect,avgspect,cdate,mean_sigma_fn)
    endif
    return
    end subroutine convert_to_sigma
! ----------------------------------------------------------------------------
    subroutine writeout(header,unspect,spect,cdate,fn)
!  WRITEOUT - write a sigma file, taking the spectral coefficients of
!  unassimilated variables from UNSPECT and merging them with SPECT,
!  which contains the spectral coefficients of the assimilated variables.
!  Arguments:
!  HEADER :=: header to place in the output sigma file.  The date field
!  is updated with CDATE.
!  UNSPECT : spectral coefficients of unassimilated variables, currently
!     ozone and relative humidity.
!  SPECT :=: on input, spectral coefficients of assimilated variables;
!     unassimilated components are copied from UNSPECT on return.
!  CDATE  : character date in yyyymmddhh format; placed in HEADER>
!  FN  : name of the sigma file to which we write the data.
!  VERBLEV  : if > 1, then we output the rms mean of each dynamical variable
!     at level 5 to the standard output.
!
!  SIDE EFFECTS
!    Opens and closes OUTSWAP.
!
    type(sigma_file_header),intent(inout)::header
    type(spectral_coeffs),intent(in)::unspect
    type(spectral_coeffs),intent(inout)::spect
    character(*),intent(in)::fn,cdate
!
!  Local variables
!
    integer::idate(4)
!
!  Copy coeffs of unassimilated humidity and ozone variables; also retain
!  the geopotential height!
!
    spect%gz=unspect%gz
    spect%rq=unspect%rq
    spect%oz=unspect%oz
!
!  The input date string has the format yyyymmddhh, but the date
!  format in the GFS sigma files is hh mm dd yyyy.
!
    if(len_trim(cdate).eq.10) then
       read(cdate,'(i4,3i2)') idate(4),idate(2),idate(3),idate(1)
    else
       idate=0
    endif
    header%fhour4=0.0
    header%idate=idate
    call write_sigma_file(fn,header,spect)
    return
    end subroutine writeout
!-----------------------------------------------------------------------------
    subroutine read_grid_files(abase,analysis,nsol_expect,npart,mean)
!  READ_GRID_FILES - read one or more files containing the model grid
!  from the data assimilation code.
!  ABASE  : this is the base name for the analysis files.  This is a
!           template name if NS > 1.
!  ANALYSIS  := holds the global analysis grid.
!  NSOL_EXPECT  : number of ensemble solutions to expect
!  NPART  : The number of partitions of the assimilation output files;
!     may be adjusted on return according to the same algorithm used
!     to compute the number of splits in preanalysis.
!  MEAN  := If present, the ensemble mean is computed and stored here
!
    implicit none
    character(*),intent(in)::abase
    integer,intent(in)::nsol_expect
    integer,intent(inout)::npart
    real(GRIDKIND),intent(out)::analysis(LONF,NAVARS,LATG,LEVS,nsol_expect)
    real(GRIDKIND),intent(out),optional::mean(LONF,NAVARS,LATG,LEVS)
!
!  Local variables
!
    real(IOKIND),allocatable::buf(:,:,:,:,:)
    integer::k,n,nsol,nstrips,verblev,seqno,nassim,north,south
!
!  COVERED(i,j) is the number of times model location (i,j) is covered
!  by some assimilation grid; should be exactly 1 for every (i,j).  
!
    integer::covered(LONF,LATG)
    integer,allocatable::list(:),startlat(:),endlat(:),psouth(:),pnorth(:)
    character(len(abase)+4)::nthfile
    real(GRIDKIND),allocatable::grid(:,:,:,:,:)
!
!  Initialize.  Use the same algorithm for splitting as in preanalysis,
!  so that we compute the correct number of analysis vector files to expect.
!
    nstrips=LATG
    call imod(npart,list,nstrips,0,startlat,endlat,psouth,pnorth)
    verblev=verbosity()
    if(verblev.gt.VERB_SILENT) write(IERR,100) trim(abase)
100 format('postanalysis: input: ',a)
!
!  Loop through all partial analysis files.
!
    if(npart.eq.1) then
       call read_grid_header(abase,'none',south,north,psouth(1),pnorth(1), &
          nsol,nassim,seqno)
       if(nsol.ne.nsol_expect) then
          write(IERR,103) trim(abase),nsol,nsol_expect
103       format(a,': contains ',i0,' solutions; expected ',i0,' members')
          call fexit(trim(abase)//': unexpected number of solutions',2)
       else if(seqno.ne.seqid()) then
          write(IERR,104) trim(abase),seqno,seqid()
104       format(a,' has sequence number: ',i0,' expected: ',i0)
       endif
       allocate(buf(LONF,NAVARS,LATG,LEVS,nsol_expect))
       read(INNOSWAP) buf
       close(INNOSWAP)
       covered=1
       analysis=buf
    else
       analysis=0.0
       covered=0
       do n=1,npart
          call numeric_filename(abase,'%n',n,nthfile)
          if(verblev.gt.VERB_TERSE) write(IERR,105) trim(nthfile)
105       format('postanalysis: reading ',a)
          call read_grid_header(nthfile,expected_tag,south,north, &
             psouth(n),pnorth(n),nsol,nassim,seqno)
          if(nsol.ne.nsol_expect) then
             write(IERR,103) trim(nthfile),nsol,nsol_expect
             call fexit(trim(nthfile)//': unexpected number of solutions',2)
          else if(seqno.ne.seqid()) then
             write(IERR,104) trim(abase),seqno,seqid()
          endif
          if(nassim.ne.NAVARS) &
             call fexit(trim(nthfile)//': assimilated var types must be 4',2)
          if(verblev.gt.VERB_TERSE) write(IERR,107) psouth(n),pnorth(n)
!
          allocate(grid(LONF,NAVARS,south:north,LEVS,nsol))
          call read_model_grid(grid,south,north,nsol)
          analysis(:,:,south:north,:,:)=grid
          covered(:,south:north)=covered(:,south:north)+1
          deallocate(grid)
       enddo
    endif
107 format('assimilation region: S:',i4,' N:',i4)
!
!  Ensure that every grid point is covered exactly once.
!
    if(any(covered.ne.1)) &
       call fexit('missing or overlapping grid centers in data',errcode=2)
!
!  Compute the mean of the analysis fields if requested.
!
    if(present(mean)) then
       mean=analysis(:,:,:,:,1)
       do k=2,nsol
          mean=mean+analysis(:,:,:,:,k)
       enddo
       mean=mean/nsol
    endif
    return
    end subroutine read_grid_files
!----------------------------------------------------------------------------
    subroutine read_perturbation_file(fn,grid)
!  READ_PERTURBATION_FILE - handle the special startup case where the
!  initial condition is a physical grid in GFS format.
!  These files have wind in meters per second and the surface pressure
!  in ln(dPa).
!
    use mrf_ids   ! from support directory
    character(*),intent(in)::fn
    real(GRIDKIND),intent(out)::grid(LONF,NAVARS,LATG,LEVS)
!
    real(D)::obsgrid(LONFD,2,NVAR,LATG2,LEVS)
    integer::ierr,j,lev
!
    open(INSWAP,file=fn,form='unformatted',status='old',iostat=ierr)
    if(ierr.ne.0) call fexit('read_perturbation_file: cannot open '//trim(fn))
    if(verbosity().gt.VERB_TERSE) call finfo('preanalysis: reading '//trim(fn))
!
    read(INSWAP) obsgrid
    close(INSWAP)
!
    do lev=1,LEVS
       do j=1,LATG2
          !
          !  South Pole to Equator
          !
          grid(:,PRES_TYPE,j,lev)= &
             real(10.0*exp(obsgrid(1:LONF,SH,Q,j,lev)),kind(grid))
          grid(:,TEMP_TYPE,j,lev)=real(obsgrid(1:LONF,SH,TE,j,lev),kind(grid))
          grid(:,UWIND_TYPE,j,lev)=real(obsgrid(1:LONF,SH,UG,j,lev),kind(grid))
          grid(:,VWIND_TYPE,j,lev)=real(obsgrid(1:LONF,SH,VG,j,lev),kind(grid))
          !
          !  Equator to North Pole; otherwise identical conversions
          !
          grid(:,PRES_TYPE,j+LATG2,lev)= &
             real(10.0*exp(obsgrid(1:LONF,NH,Q,LATG2-j+1,lev)),kind(grid))
          grid(:,TEMP_TYPE,j+LATG2,lev)= &
             real(obsgrid(1:LONF,NH,TE,LATG2-j+1,lev),kind(grid))
          grid(:,UWIND_TYPE,j+LATG2,lev)= &
             real(obsgrid(1:LONF,NH,UG,LATG2-j+1,lev),kind(grid))
          grid(:,VWIND_TYPE,j+LATG2,lev)= &
             real(obsgrid(1:LONF,NH,VG,LATG2-j+1,lev),kind(grid))
       enddo
    enddo
    return
    end subroutine read_perturbation_file
!--------------------------------------------------------------------------
    subroutine read_surface_files(template,avg,fhour,nsol)
!  READ_SURFACE_FILES - read a family of surface analysis files from the MRF
!  code and compute a mean surface analysis file.
!  Arguments:
!  TEMPLATE : template name for surface analysis files.
!  AVG := averaged surface analysis, computed from all the inputs.
!  FHOUR := forecast hour represented in the input surface analysis files.
!  NB : number of ensemble solutions (and therefore surface analysis files)
!     to expect.
!  In the case of the MRF, all surface analysis files are in big-endian format
!  and so may require byte-swapping on little-endian architectures.
!
    use surface
    implicit none
    character(*),intent(in)::template
    type(sfcanl),intent(out)::avg
    integer,intent(in)::nsol
    real(R),intent(out)::fhour
!
!  Local variables
!
    type(sfcanl)::x
    real(R)::hour
    integer::j,idate(4),refdate(4)
    character(MAX_FILE_NAME_LEN)::fn
!
!  Use the first file for initialization, then sum over the remaining ones.
!
    call numeric_filename(template,'%n',1,fn)
    call read_sfcanl(fn,avg,fhour,refdate)
    do j=2,nsol
       call numeric_filename(template,'%n',j,fn)
       call read_sfcanl(fn,x,hour,idate)
       if(hour.ne.fhour.or.any(idate.ne.refdate)) &
          call fexit('read_surface_files: '//trim(fn)// &
             ' forecast hour and/or dates disagree',2)
       avg=avg+x
    enddo
    avg=avg/nsol
    return
    end subroutine read_surface_files
!------------------------------------------------------------------------------
    subroutine parseline(cdate,nsol)
!  PARSELINE - parse an execute line each of whose arguments is 
!  of the form option=value.
!
    implicit none
    character(*),intent(inout)::cdate
    integer,intent(out)::nsol
!
!  Local variables
!
    character(*),parameter::OPTLIST='d:f:hH:m:n:N:Pp:S:s:v:'
    character(MAX_FILE_NAME_LEN)::optarg,nmlfile
    integer::nargs,iopt,optind,verblev
    integer,external::iargc  ! g95 intrinsic
    character(10)::chardate
    namelist/input/ chardate,ref_sig_fn,mean_sigma_fn,nsol,nsplit, &
       sigma_out_fn,sfc_fn,verblev,phys_grid_fn
!
    nargs=iargc()  ! number of command-line arguments EXCLUDING program name
    nmlfile=' '
    chardate=cdate
    if(nargs.lt.2) call usage
    do 
       call getopt(iopt,optarg,optind,OPTLIST)
       if(iopt.lt.0) exit
       select case(char(iopt))
       case('d')         ! date to use in yyyymmddhh format
          cdate=optarg
       case('f')
          ref_sig_fn=optarg
       case('h')
          call usage
       case('H')
          expected_tag=optarg
       case('m')
          mean_sigma_fn=optarg
       case('N')    ! namelist input comes from the named file
          nmlfile=optarg
          open(IN,file=nmlfile,status='old')
          read(IN,input)
          close(IN)
          cdate=chardate
       case('n') 
          read(optarg,*) nsol
       case('P')   ! input is physical perturbation files
          input_is_pert=.true.
       case('p')
          read(optarg,*) nsplit
       case('s')
          sigma_out_fn=optarg
       case('S')
          sfc_fn=optarg
       case('v')
          read(optarg,*) verblev
          call set_verbosity(verblev)
       endselect       ! otherwise, getopt prints message
    enddo
!
!  Get the required argument if namelist input isn't specified
!
    if(len_trim(nmlfile).eq.0) then
       if(nargs.lt.optind) call usage
       call getarg(nargs,phys_grid_fn)
    endif
!
!  Brief sanity checks
!
    if(len_trim(ref_sig_fn).eq.0) &
       call fexit('-f must be used to specify a sigma file')
    if(nsol.lt.1) call fexit('nsol must be > 0')
!
!  Brief announcement if requested
!
    if(verblev.ge.VERB_TERSE) then
       write(LOGFILE,100) 'postanalysis: version ',VERSION
100    format(2a)
       write(LOGFILE,100) 'input grid: ',trim(phys_grid_fn)
       write(LOGFILE,100) 'ref sigma file (-f): ',trim(ref_sig_fn)
       write(LOGFILE,100) 'output sigma files: ',trim(sigma_out_fn)
       if(len_trim(mean_sigma_fn).gt.0) &
          write(LOGFILE,100) 'mean sigma file: ',trim(mean_sigma_fn)
       if(len_trim(sfc_fn).gt.0) &
          write(LOGFILE,100) 'surface analysis files: ',trim(sfc_fn)
    endif
    return
!---------
    contains
       subroutine usage
       write(IERR,100) VERSION
100    format('postanalysis: version ',a/ &
          'usage: [options] grid_file, where options are'/ &
          ' [-f ref_sigma_file] (required option; gives spectral coeffs)'/ &
          ' [-N input_file] - all options taken from namelist in input_file'/ &
          ' [-n number_of_ensembles]'/ &
          ' [-d date in yyyymmddhh format]'/ &
          ' [-s output_sigma_file_template]'/ &
          ' [-P (input is physical perturbation file)]'/ &
          ' [-S surface_analysis_file_template; mean placed in 0]'/ &
          ' [-H expected_tag; default "aout"]'/ &
          ' [-v verblev]'/ &
          ' [-p nsplit (number of pieces for input physical grid)]'/ &
          ' [-m output_mean_sigma_file]'/ &
          '-s specifies a family of sigma files for output, one per'/ &
          '   ensemble member; default: sigma.outN'/ &
          'grid_file specifies the name of a physical model grid file'/ &
          '   in latitude-major format')
       call program_exit(1)
       end subroutine usage
    end subroutine parseline
! ----------------------------------------------------------------------------
    end module postanalysis
! -----------------------------------------------------------------------------
!  MAIN.
!
    use postanalysis
!
    call do_everything
    call program_exit(0)
    end
