!  GRIDDUMP - dump the contents of a global grid file.
!  Optionally compute and display relative or absolute differences
!  between two global grid files.
!  Modification history:
!  08/08/05 - Now relies on GRIDKIND definition in support modules
!  07/26/05 - added -d option
!  07/24/05 - added old-style formatting for easier comparison with 2004 version
!  07/16/05 - version to handle energy units.
!  06/19/05 - added options to display ensemble mean and std error.
!             Deleted option for scaled displays since all output is in
!             physical units nowadays.
!  06/17/05 - changes to make the program a little smarter.
!  03/07/05 - change to accommodate new model grid i/o format.
!  11/29/04 - change to accommodate new model grid format.
!             Changed -N to -s to dump scaled values (default is unscaled).
!  05/18/04 - change the default to output values in physical units.
!             Added -N option to dump unscaled data.
!  08/31/03 - version to handle latitude-oriented indexing of grid
!             arrays.  However, for compatibility with older versions,
!             we output data with longitude indexes most rapidly varying.
!--------------------------------------------------------------------
    module dumpdefs
    use grid_io
    use spectral
    implicit none
    character(*),parameter::VERSION='08/08/05'
    integer,save::verblev=1
    logical,save::output_stats=.false.  ! by default, output each solution
    character(MAX_FILE_NAME_LEN),save::sigmafn=' '  ! no ref file by default
    character(MAX_FILE_NAME_LEN),save::diff_fn=' '  ! no diff file by default
    logical,save::model_obs_format=.false.  ! true => output like 2004 version
    logical,save::physical_units=.false.  ! default: energy units
    contains
!-------------------------------------------------------------------
    subroutine read_grid(fn,expect,nsol,south,north,grid)
!  If NSOL > 0, then it is an error if the input file does not contain
!  exactly NSOL solutions.
!
    integer,intent(inout)::nsol
    integer,intent(out)::south,north
    character(*),intent(in)::fn,expect
    real(GRIDKIND),allocatable,intent(out)::grid(:,:,:,:,:)
!
    integer::nassim,psouth,pnorth,seqno,nb
    character(80)::msg
!
!  Read grid.  This requires care: if the number of ensemble solutions
!  represented to the program differs from the number in the files,
!  then an unexpected EOF may occur.
!
    call read_grid_header(fn,expect,south,north,psouth,pnorth,nb,nassim,seqno)
!
    if(nassim.ne.NAVARS) then
       write(msg,101) nassim
101    format('griddump: header alleges ',i0,' variables; expected 4')
       call fexit(trim(msg),2)
    else if(nsol.le.0) then
       nsol=nb
    else if(nb.ne.nsol) then
       write(msg,102) nsol,nb
102    format('griddump: expected ',i3,' ensemble solutions; got ',i3)
       call fexit(trim(msg),2)
    endif
!
    allocate(grid(LONF,NAVARS,south:north,LEVS,nsol))
    call read_model_grid(grid,south,north,nsol)
    if(verblev.gt.0) then
       write(LOGFILE,110) trim(fn)
       write(LOGFILE,111) south,north,psouth,pnorth,nsol,seqno
    endif
110 format('file: ',a)
111 format('south: ',i0,' north: ',i0,' psouth: ',i0,' pnorth: ',i0/ &
       ' nsol: ',i0,' seqno: ',i0)
    return
    end subroutine read_grid
!-------------------------------------------------------------------
    subroutine dump_grid_stats(grid,south,north,nsol,lon,lat)
    integer,intent(in)::south,north,nsol,lon,lat
    real(GRIDKIND),intent(in)::grid(LONF,NAVARS,south:north,LEVS,nsol)
!
    integer::k,n,lev
    real(GRIDKIND)::mean(LONF,NAVARS,south:north,LEVS)
    real(GRIDKIND)::sdev(LONF,NAVARS,south:north,LEVS)
    real(GRIDKIND)::big(LONF,NAVARS,south:north,LEVS)
    real(GRIDKIND)::small(LONF,NAVARS,south:north,LEVS)
!
    write(LOGFILE,200)
200 format('#what lon  lat  lev  P  T  U  V')
!
    mean=0.0
    sdev=0.0
    big=-huge(big)
    small=huge(small)
!
    if(lon.gt.0) then
       if(lat.gt.0) then  ! dump stats at (lon,lat)
          if(lat.lt.south.or.lat.gt.north) &
             call fexit('specified latitude outside available strip',1)
          do n=1,nsol
             mean(lon,:,lat,:)=mean(lon,:,lat,:)+grid(lon,:,lat,:,n)
             sdev(lon,:,lat,:)=sdev(lon,:,lat,:)+grid(lon,:,lat,:,n)**2
             big(lon,:,lat,:)=max(big(lon,:,lat,:),grid(lon,:,lat,:,n))
             small(lon,:,lat,:)=min(small(lon,:,lat,:),grid(lon,:,lat,:,n))
          enddo
          mean(lon,:,lat,:)=mean(lon,:,lat,:)/nsol
          sdev(lon,:,lat,:)=sqrt(sdev(lon,:,lat,:)/nsol-mean(lon,:,lat,:)**2)
          write(LOGFILE,100)
100       format('#what lon  lat   lev   P  T  U  V')
          do lev=1,LEVS
             write(LOGFILE,101) 'mean',lon,lat,lev,mean(lon,:,lat,lev)
             write(LOGFILE,101) 'sdev',lon,lat,lev,sdev(lon,:,lat,lev)
             write(LOGFILE,101) 'max ',lon,lat,lev,big(lon,:,lat,lev)
             write(LOGFILE,101) 'min ',lon,lat,lev,small(lon,:,lat,lev)
          enddo
101       format(a,3i4,4es11.3)
       else                   ! dump stats at (lon,:)
          do n=1,nsol
             mean(lon,:,:,:)=mean(lon,:,:,:)+grid(lon,:,:,:,n)
             sdev(lon,:,:,:)=sdev(lon,:,:,:)+grid(lon,:,:,:,n)**2
             big(lon,:,:,:)=max(big(lon,:,:,:),grid(lon,:,:,:,n))
             small(lon,:,:,:)=min(small(lon,:,:,:),grid(lon,:,:,:,n))
          enddo
          mean(lon,:,:,:)=mean(lon,:,:,:)/nsol
          sdev(lon,:,:,:)=sqrt(sdev(lon,:,:,:)/nsol-mean(lon,:,:,:)**2)
          write(LOGFILE,100)
          do lev=1,LEVS
             do k=south,north
                write(LOGFILE,101) 'mean',lon,k,lev,mean(lon,:,k,lev)
                write(LOGFILE,101) 'sdev',lon,k,lev,sdev(lon,:,k,lev)
                write(LOGFILE,101) 'max ',lon,k,lev,big(lon,:,k,lev)
                write(LOGFILE,101) 'min ',lon,k,lev,small(lon,:,k,lev)
             enddo
          enddo
       endif
    else if(lat.gt.0) then  ! dump stats at (:,lat)
       do n=1,nsol
          mean(:,:,lat,:)=mean(:,:,lat,:)+grid(:,:,lat,:,n)
          sdev(:,:,lat,:)=sdev(:,:,lat,:)+grid(:,:,lat,:,n)**2
          big(:,:,lat,:)=max(big(:,:,lat,:),grid(:,:,lat,:,n))
          small(:,:,lat,:)=min(small(:,:,lat,:),grid(:,:,lat,:,n))
       enddo
       mean(:,:,lat,:)=mean(:,:,lat,:)/nsol
       sdev(:,:,lat,:)=sqrt(sdev(:,:,lat,:)/nsol-mean(:,:,lat,:)**2)
       write(LOGFILE,100)
       do lev=1,LEVS
          do k=south,north
             write(LOGFILE,101) 'mean',lon,k,lev,mean(lon,:,k,lev)
             write(LOGFILE,101) 'sdev',lon,k,lev,sdev(lon,:,k,lev)
             write(LOGFILE,101) 'max ',lon,k,lev,big(lon,:,k,lev)
             write(LOGFILE,101) 'min ',lon,k,lev,small(lon,:,k,lev)
          enddo
       enddo
    else              ! dump stats at all locations
       do n=1,nsol
          mean=mean+grid(:,:,:,:,n)
          sdev=sdev+grid(:,:,:,:,n)**2
          big(:,:,:,:)=max(big(:,:,:,:),grid(:,:,:,:,n))
          small(:,:,:,:)=min(small(:,:,:,:),grid(:,:,:,:,n))
       enddo
       mean=mean/nsol
       sdev=sqrt(sdev/nsol-mean**2)
       do n=1,LONF
          do k=south,north
             do lev=1,LEVS
                write(LOGFILE,101) 'mean',n,k,lev,mean(n,:,k,lev)
                write(LOGFILE,101) 'sdev',n,k,lev,sdev(n,:,k,lev)
                write(LOGFILE,101) 'max ',n,k,lev,big(n,:,k,lev)
                write(LOGFILE,101) 'min ',n,k,lev,small(n,:,k,lev)
             enddo
          enddo
       enddo
    endif
    return
    end subroutine dump_grid_stats
!-------------------------------------------------------------------
    subroutine dump_grid(grid,nsol,south,north,lon,lat)
    integer,intent(in)::nsol,south,north,lon,lat
    real(GRIDKIND),intent(in)::grid(LONF,NAVARS,south:north,LEVS,nsol)
!
    integer::j,k,n
!
    if(lon.gt.0) then
       if(lat.gt.0) then  ! dump solutions at (lon,lat)
          if(lat.lt.south.or.lat.gt.north) &
             call fexit('specified latitude outside available strip',1)
          do n=1,nsol
             write(LOGFILE,110) n
110          format('Ensemble solution ',i0)
             call dump_grid_point(grid,nsol,south,north,lon,lat,n)
          enddo
       else                   ! dump solutions at (lon,:)
          do n=1,nsol
             write(LOGFILE,110) n
             do k=south,north
                call dump_grid_point(grid,nsol,south,north,lon,k,n)
             enddo
          enddo
       endif
    else if(lat.gt.0) then  ! dump solutions at (:,lat)
       do n=1,nsol
          write(LOGFILE,110) n
          do k=1,LONF
             call dump_grid_point(grid,nsol,south,north,k,lat,n)
          enddo
       enddo
    else              ! dump all solutions
       do n=1,nsol
          write(LOGFILE,110) n
          do j=south,north
             do k=1,LONF
                call dump_grid_point(grid,nsol,south,north,k,j,n)
             enddo
          enddo
       enddo
    endif
    return
    end subroutine dump_grid
!-----------------------------------------------------------------------
    subroutine dump_grid_point(g,nsol,south,north,ix,iy,n)
    integer,intent(in)::ix,iy,south,north,nsol,n
    real(GRIDKIND),intent(in)::g(LONF,NAVARS,south:north,LEVS,nsol)
!
    integer::lev
    logical,save::first=.true.
!
    if(model_obs_format) then  ! output all levels together
       write(LOGFILE,100) ix,iy
100    format(' Longitude, latitude',2i4)
       write(LOGFILE,101) 'surface pressure',g(ix,PRES_TYPE,iy,1,n)
       write(LOGFILE,101) 'temp',g(ix,TEMP_TYPE,iy,:,n)
       write(LOGFILE,101) 'uwind',g(ix,UWIND_TYPE,iy,:,n)
       write(LOGFILE,101) 'vwind',g(ix,VWIND_TYPE,iy,:,n)
101    format(3x,a/(3x,7es11.3))
    else
       if(first) write(LOGFILE,110)
110    format('#lon  lat  lev  var values (for all ensemble members)')
       do lev=1,LEVS
          if(lev.eq.1) write(LOGFILE,111) ix,iy,lev,'P',g(ix,PRES_TYPE,iy,lev,:)
          write(LOGFILE,111) ix,iy,lev,'T',g(ix,TEMP_TYPE,iy,lev,:)
          write(LOGFILE,111) ix,iy,lev,'U',g(ix,UWIND_TYPE,iy,lev,:)
          write(LOGFILE,111) ix,iy,lev,'V',g(ix,VWIND_TYPE,iy,lev,:)
       enddo
111    format(3i4,1x,a1,1x,5es11.3/('> ',7es11.3))
    endif
    first=.false.
    return
    end subroutine dump_grid_point
!-----------------------------------------------------------------------------
    subroutine parseline(fn,header,nsol,locations)
!  PARSELINE - parse execute line.  Getopt works much like the C version
!  (man 3c getopt).
!  The arguments are the background, analysis, observations, and analysis mean
!  filenames, respectively.
!  Side effects:
!  - calls patch_init to initialize patch-related variables
!
    use command_line_support
    integer,intent(out)::nsol
    character(*),intent(out)::fn,locations
    character(*),intent(inout)::header
!
!  Local variables
!
    character(*),parameter::OPTIONS='d:f:mn:phl:v:H:s'
    character(80)::optarg
    integer::iopt,optind
    character::opt
!
!  Copy defaults
!
    fn=' '
    locations=' '
    nsol=0  ! get the number of solutions from input file unless -n is specified
    do
       call getopt(iopt,optarg,optind,options)
       if(iopt.lt.0) exit
       opt=char(iopt)
       select case(opt)
       case('d')  ! compute differences from selected grid
          diff_fn=optarg
       case('f')  ! implies physical units; take sigma file header from here
          sigmafn=optarg
          physical_units=.true.
       case('m')
          model_obs_format=.true.
       case('n')
          read(optarg,100) nsol
       case('h')
          write(IERR,120) VERSION
          call exit(1)
       case('p')
          physical_units=.true.  ! and use default sigma file header
       case('l')
          locations=optarg
       case('s')
          output_stats=.true.
       case('v')
          read(optarg,100) verblev
       case('H')
          header=optarg
       case default
          write(IERR,110) opt
          call exit(1)
       end select
    enddo
100 format(i12)
110 format('unknown option ',a)
!
    call getarg(optind,fn)
120 format('griddump: version ',a/ &
       ' usage: options gridfile'/ &
       ' options include: [-l long,lat] (dump at selected location)'/ &
       ' [-f ref_sigma_file (used to convert to physical units)]'/ &
       ' [-H header] (header to expect; normally "grid" or "aout")'/ &
       ' [-m (model_obs format: all levels together)]'/ &
       ' [-s (output stats)] (output mean and std error only)')
    return
    end subroutine parseline
!--------------------------------------------------------------------
    subroutine get_locations(str,ix,iy)
!  GET_LOCATIONS - str must have the form ix,iy or :,iy or ix,:
!
    character(*),intent(in)::str
    integer,intent(out)::ix,iy
!
    integer::icom 
!
    icom=index(str,',')
    if(icom.eq.0) then   ! bad string; forget it
       ix=0
       iy=0
       return
    endif
    if(str(1:icom-1).eq.':') then
       ix=0
    else
       read(str(1:icom-1),*) ix
    endif
    if(str(icom+1:).eq.':') then
       iy=0
    else
       read(str(icom+1:),*) iy
    endif
    return
    end subroutine get_locations
!--------------------------------------------------------------------
    end module dumpdefs
!--------------------------------------------------------------------
! MAIN.
    use dumpdefs
    implicit none
    character(128)::fn,locations
    character(4)::header='none'  ! by default, assume global grid
    integer::nsol,ix,iy,south,north,j,ndiff
    real(GRIDKIND),allocatable::grid(:,:,:,:,:),diffgrid(:,:,:,:,:)
       ! (LONF,NAVARS,south:north,LEVS,nsol)
    type(sigma_file_header)::sfhdr
!
    call parseline(fn,header,nsol,locations)
    if(len_trim(locations).gt.0) call get_locations(locations,ix,iy)
!
!  The grid is assumed to be in energy units; so convert to physical units
!  (pressure in hPa, temperature in K, and wind in m/s).  Both sets of units
!  maintain pressure in ln(dPa); for output purposes, convert to hPa.
!
    call read_grid(fn,header,nsol,south,north,grid)
    if(len_trim(diff_fn).gt.0) then  ! do analysis on difference grid
       call read_grid(diff_fn,header,ndiff,south,north,diffgrid)
       if(ndiff.eq.1) then  ! DIFFGRID is nature or mean grid, presumably
          do j=1,nsol
             grid(:,:,:,:,j)=grid(:,:,:,:,j)-diffgrid(:,:,:,:,1)
          enddo
       else if(ndiff.ne.nsol) then
          call fexit(trim(diff_fn)//': incompatible ensemble size')
       else
          grid=grid-diffgrid
       endif
    endif
    if(physical_units) then
       if(len_trim(sigmafn).gt.0) then
          call spectral_code_init(sigmafn,sfhdr)
       else
          call spectral_code_init  ! use defaults
       endif
       do j=1,nsol
          call from_energy_units(grid(:,:,:,:,j),south,north)
       enddo
       if(allocated(diffgrid)) then  ! show zero as 0.0, not 10 = 10*exp(0)
          where(grid(:,PRES_TYPE,:,:,:).gt.1.0e-03)
             grid(:,PRES_TYPE,:,:,:)=10.0*exp(grid(:,PRES_TYPE,:,:,:)) ! to hPa
          elsewhere
             grid(:,PRES_TYPE,:,:,:)=10.0*grid(:,PRES_TYPE,:,:,:)
          endwhere
       else
          grid(:,PRES_TYPE,:,:,:)=10.0*exp(grid(:,PRES_TYPE,:,:,:))  ! to hPa
       endif
    endif
    if(output_stats) then
       call dump_grid_stats(grid,south,north,nsol,ix,iy)
    else
       call dump_grid(grid,nsol,south,north,ix,iy)
    endif
    call exit(0)
    end
!--------------------------------------------------------------------
