    module preanalysis
!  PREANALYSIS - convert from sigma files to physical grid, and split
!  along latitude strips for efficient parallel operation on a Beowulf cluster.
!  This version is intended for the 2001 or 2004 versions of the GFS model.
!
! SYNPOSIS:
!    [options] input_file
!
!  There are two modes of operation.
!  1)  We need to create a set of observations, where we specify -o.
!      The input_file argument is the name of an observation file.
!      Only very simple preprocessing is done: we simply copy those observations
!      that belong to the appropriate latitudinal strip.
!  2)  We need to create a grid of model variables from a sigma file.
!      In this case, the required input_file argument
!      is a sigma file, which is assumed to contain the required spectral
!      coefficients and model variables in spectral coordinates.
!      input_file is a template name if NSOL > 1 (see -n below).
!      The output is a grid file containing the data structures
!      used in the data assimilation code.
!
!  A template name describes how the program should examine a numerical
!  sequence of files.  Given a string like "input", the program looks
!  for files name input1, input2, etc., if more than one input is expected.
!  Alternatively, the template name can be something like input%n.dat,
!  where the string '%n' is substituted with 1, 2, ..., as required.
!  In this case, we look for files named input1.dat, input2.dat, etc.
!
!  OPTIONS:
!  -b name  : NAME is the name of the output file containing the model grid 
!      data.  It is a template name if NP > 1 (see -p below).
!  -P : specifies that the input is not a GFS sigma file but a physical grid 
!      perturbation file.
!  -o name  : NAME is the name of the output file containing the observations in
!      the format needed by the data assimilation code.
!      It is a template name if NP > 1 (see -p below).
!  -s name  : NAME is the name of an input sigma file to which the various 
!      coefficients are written for conversion between spectral and 
!      physical units.  The output is in text format.
!  -v lev  : (default: lev=0)  set the verbosity to level LEV.
!      0 gives silent operation, 1 gives short details, 2 shows the status
!      of every operation.
!  -y yrad :  Make output files suitable for use with patches whose
!      vertical (latitudinal) radius is <= YRAD.
!  -n nsol : (default: nsol=40)  The number of ensemble solutions.
!      Incompatible with -o.  If NSOL > 1, then input_file is a template name.
!  -H header : (default: none)  Use the indicated 4-char header when splitting
!      files for mulitple processors.  The defaults are 'obs ' for observations
!      and 'grid' for the model grid.
!  -p np    : (default: np=2)  The number of partitions that we should make
!      of the output files.  This value is the number of processors on which
!      we expect to run the data assimilation step.
!  -q seqid : sequence id (usually the date in yyyymmddhh format) to expect
!      from the input files and the sequence id to place in the output.
!
!   input_file is a sigma file if -b (and/or -n) are specified.
!   input_file is an observation file of -o is specified.
!
! OTHER INFO:
!  -This program uses INNOSWAP and INSWAP for unformatted input and OUTNOSWAP
!  for unformatted output, as defined in module io.  Use the FORT90L
!  environment variable to specify whether byte-swapping should occur on
!  logical unit INNOSWAP (Lahey only).
!
!  -Although the GFS model uses arrays with leading dimension 1:LONFD,
!  we create, write, and use only arrays with longitude dimension 1:LONF.
!  Zero pads are written as needed on output, and the information in
!  rows LONF+1:LONFD is stripped on input.
!
! ATTRIBUTES:
!    Language: Fortran 95
!    Machine : Generic (except for optional FORT90L usage, specific to Lahey)
!
!  Modification history
!    09/10/05  Altered documentation and version to reflect compatibility
!              with both the 2001 and 2004 versions of the GFS (I hope).
!    06/29/05  Output grids in energy units instead of physical units.
!    06/17/05  New version that calls overhauls spectral code and re-introduces
!              support for reading the physical grid perturbation files.
!    05/29/05  Added option for namelist input (for debugging) and made changes
!              to reflect alterations of spectral code to handle the 2004 GFS.
!    03/07/05  Used revised spectral code to fix N/S hemisphere indexing
!    01/13/05  Added -A flag for simple selection of observation types to
!              distribute to the nodes
!    12/15/04  Modifications for more efficient GFS -> model grid processing:
!              Changed array subscript ordering for model grid.
!              Also removed special initialization step, since I can't figure
!              out how the original perturbation grid was scaled.
!              Instead, we ALWAYS start with a set of sigma files.
!              Runs with this version will require the sigma files output
!              by the initialization step of the old (~/nws) code version.
!    11/24/04  Initial version by Eric and Istvan for real observations.
!              The program processes raw_obs formats and simply copies
!              each observation to the appropriate latitude strip for parallel
!              operation of the data assimilation program.
!    07/22/04  Eric Kostelich, changed handling of unit numbers for
!              byte-swapping in the Lahey Linux implementation.
!    10/10/03  Eric Kostelich, optionally writes out coefficients for
!              conversion to physical units; for diagnostics package
!    08/23/03  Eric Kostelich, removed support for averaging and changed
!              to latitude-oriented strips
!    04/29/03  Eric Kostelich, bug fixes
!    04/17/03  Eric Kostelich, modified version for Beowulf cluster
!    10/04/02  Istvan Szunyogh, original version
!
    use spectral
    use readobs
    use resolution
    use io
    use constant
    use command_line_support
    implicit none
!
    character(*),parameter::VERSION='09/10/05 (2001/04 models, energy units)'
!
!  By default, the input is a sigma file; but if -P is specified, then
!  the input is a physical grid perturbation file.
!
    logical,save::input_is_sigma=.true.
!
!  By default, the output is an array of model variables for use in the
!  data assimilation code (-b option).
!
    logical,save::output_is_grid=.true.
!
!  Required input file specification - either a sigma file or a file
!  of observations.  It is a template name if processing more than one 
!  sigma file in an ensemble of solutions.
!
    character(MAX_FILE_NAME_LEN),save::infn=' '    ! no default
!
!  The name of the output file (-b or -o)
!
    character(MAX_FILE_NAME_LEN),save::outfn='egrid'
!
!  Vertical (latitudinal) radius for which outputs can be used
!
    integer,save::ywin=3
!
!  For -H option
!
    character(4),save::optional_header=' '
!
!  Number of partitions of the output data for the assimilation step.
!
    integer,save::nsplit=2
    contains
!----------------------------------------------------------------------------
    subroutine do_everything
!  DO_EVERYTHING - so we don't have to look for the main program at the
!  end of the file.
! 
    integer::nsol=40  ! default number of ensemble solutions
    character(4)::list='PTUV'  ! types of observations to assimilate
! 
    call parseline(nsol,list)
    if(verbosity().gt.VERB_SILENT) write(LOGFILE,100) version
100 format('preanalysis: version ',a)
!
!  Process the model grid or the observations as required.
!
    if(output_is_grid) then
       if(verbosity().gt.VERB_TERSE) write(LOGFILE,101) nsol
101    format('preanalysis: nsol: ',i0)
       call create_grid(nsol)
    else
       call read_observations(infn,outfn,nsplit,ywin,list)
    endif
    return
    end subroutine do_everything
!---------------------------------------------------------------------------
    subroutine create_grid(nsol)
!  CREATE_GRID - create the physical model grid from a collection of
!  sigma files.
!  Side effects:  Opens and closes INSWAP and OUTNOSWAP.
!
    use fsplit
    integer,intent(in)::nsol
!
!  Local variables
!
    integer::k
    real(GRIDKIND),allocatable::egrid(:,:,:,:,:)
    character(MAX_FILE_NAME_LEN)::fn
    type(sigma_file_header)::header
    type(spectral_coeffs)::spect
!
!  Initialize the spectral transform package using the first (or only)
!  sigma file in the series.
!
    if(input_is_sigma) then
       if(nsol.eq.1) then    ! the given sigma file will do
          fn=infn
       else
          call numeric_filename(infn,'%n',1,fn)
       endif
       call spectral_code_init(fn,header)
    endif
!
    allocate(egrid(LONF,NAVARS,LATG,LEVS,nsol))
    if(verbosity().gt.VERB_SILENT) write(LOGFILE,111) trim(infn)
111 format('preanalysis: input sigma file(s): ',a)
!
!  Process all input files to get the physical grid for all ensemble members.
!
    do k=1,nsol
       if(nsol.eq.1) then  ! not a template name
          fn=infn
       else
          call numeric_filename(infn,'%n',k,fn)
       endif
       if(input_is_sigma) then
          call read_sigma_file(fn,header,spect)
          call to_physical_space(spect,egrid(1,1,1,1,k))
       else
          call read_perturbation_file(fn,egrid(1,1,1,1,k))
       endif
    enddo
!
    if(len_trim(optional_header).eq.0) optional_header='grid'
    if(verbosity().gt.VERB_SILENT) then
       write(LOGFILE,112) trim(outfn)
112    format('preanalysis: output: ',a)
       if(len_trim(optional_header).gt.0) &
          write(LOGFILE,113) trim(optional_header)
113    format('preanalysis: header label: ',a)
    endif
!
!  Split the physical grid into pieces for the other processors
!
    call grid_split(egrid,outfn,nsol,nsplit,ywin,optional_header)
    return
    end subroutine create_grid
!----------------------------------------------------------------------------
    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('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 parseline(ensemble_size,list)
!  PARSELINE - parse the command line.
!
    implicit none
    integer,intent(inout)::ensemble_size
    character(*),intent(inout)::list
!
!  Local variables
!
    character(*),parameter::OPTLIST='A:b:H:hN:n:o:Pp:q:v:y:'
    character(MAX_FILE_NAME_LEN)::optarg,nmlfile
    character(8)::assimlist
    integer::nargs,iopt,optind,verblev,nsol,iseq
    logical::have_b,have_o
    namelist/input/ assimlist,infn,outfn,nsol,nsplit,verblev,ywin, &
       optional_header,output_is_grid
    integer,external::iargc  ! g95 intrinsic
!
!  NARGS is the number of command-line arguments EXCLUDING the program name
!
    nargs=iargc()
    have_b=.false.
    have_o=.false.
    nmlfile=' '
    assimlist=list
    nsol=ensemble_size
    if(nargs.lt.2) call usage
    do 
       call getopt(iopt,optarg,optind,optlist)
       if(iopt.lt.0) exit
       select case(char(iopt))
       case('A')  ! list of variables to be assimilated
          list=optarg
       case('b')
          outfn=optarg
          have_b=.true.
          output_is_grid=.true.
       case('H')
          optional_header=optarg
       case('h')
          call usage
       case('N')    ! namelist input comes from the named file
          nmlfile=optarg
          open(unit=IN,file=nmlfile,status='old')
          read(IN,input)
          close(IN)
       case('n') 
          read(optarg,*) nsol
       case('o')
          outfn=optarg
          have_o=.true.
          output_is_grid=.false.
       case('P')
          input_is_sigma=.false.
       case('p')
          read(optarg,*) nsplit
       case('q')
          read(optarg,*) iseq
          call set_sequence_id(iseq)
       case('v')
          read(optarg,*) verblev
          call set_verbosity(verblev)
       case('y')
          read(optarg,*) ywin
       case default
          call program_exit(1)     ! getopt prints message
       endselect
    enddo
!
!  Get the required argument if no namelist input is given
!
    if(len_trim(nmlfile).eq.0) then
       if(nargs.lt.optind) call usage
       call getarg(nargs,infn)
    endif
    list=assimlist
    ensemble_size=nsol
!
!  Brief sanity checks
!
    if(have_b.and.have_o) call fexit('-b and -o are mutually exclusive options')
    if(nsol.lt.1) call fexit('nsol must be > 0')
    if(nsplit.lt.1) call fexit('nsplit must be > 0')
    if(len_trim(infn).eq.0) &
       call fexit('input observation or sigma file must be specified')
    return
!
    contains
       subroutine usage
       write(IERR,101) VERSION
101    format('preanalysis: version ',a/ &
          'usage: three basic command lines:'/ &
          '1. preanalysis -b egrid [-p nproc] [-n nsol] sigma,'/ &
          '   where sigma holds the spectral coords of the '/ &
          '   GFS model grid (one file per ensemble member)'/ &
          '   and the output file family egrid#n holds the physical grid'/ &
          '   (NRPOC slices, each with NSOL members)')
       write(IERR,102)
102    format('2. preanalysis -P -b egrid [-p nproc] [-n nsol] grid,'/ &
          '   where sigref is a reference GFS sigma file'/ &
          '   and grid is a template name that specifies a family of'/ &
          '   physical grid perturbation files'/ &
          'Specify -P with -b only for the first step in a cycle.')
       write(IERR,103)
103    format('3. preanalysis -o obsgrid [-p nproc] obsdata'/ &
          '   where obsdata holds a file of "real" observations'/ &
          '   and the output file family obsgrid#n holds the NPROC'/ &
          '   slices of the obsdata, one per processor in latitude strips')
       write(IERR,104)
104    format('Other options:'/ &
          '[-A obslist] (P:pressure, T:temperature, U,V:wind components)'/ &
          '[-v verblev] (0:silent, 2:verbose)'/ &
          '[-N input_file] (read namelist data from input_file)'/ &
          '[-y padding] (default: 3)'/ &
          '[-p nsplit] (number of splits of output files: default 2)'/ &
          '[-q seqid] (sequence id, e.g., yyyymmddhh)')
       call program_exit(1)
       end subroutine usage
    end subroutine parseline
    end module preanalysis
!-------------------------------------------------------------------------
! MAIN.
!
    use preanalysis
! 
    call do_everything
    call program_exit(0)
    end
