!  SDUMP - dump the contents of a sigma file.
!  Usage:  sdump [options] file
!  Options:  -l lev  : dump only for level LEV
!            -v var  : one of gz, q, te, di, sz, rq, oz

    module sdump
    use spectral
    use io
    use command_line_support
    implicit none
    integer,parameter::HEADONLY=1,RMSONLY=2,PRINTVALS=3
    contains
!-------------------------------------------------------------------------
    subroutine parseline(fn,lev,var,action)
    character(*),intent(in)::fn
    character(*),intent(out)::var
    integer,intent(out)::lev,action
!
    character(*),parameter::OPTLIST='hHl:v:'
    integer::iopt,optind,nargs
    character(MAX_FILE_NAME_LEN)::optarg
    character::opt
    integer,external::iargc
!
    var='all'
    lev=0
    action=0
    nargs=iargc()
    if(nargs.lt.1) then
       write(IERR,100)
       call program_exit(1)
    endif
100 format('sdump [-H] [-v var] [-l level] file'/ &
       'Default is RMS output by level for all variables'/ &
       '-H dumps the sigma file header only'/ &
       '-v dumps individual coefficients for the indicated level (if given)'/ &
       'var is one of gz, q, te, di, ze, rq, or oz'/ &
       'If -l is unspecified then all levels are shown (lots of output!)'/ &
       'On Lahey, export FORT90L="-Wl,-T1"')
!
    do
       call getopt(iopt,optarg,optind,OPTLIST)
       if(iopt.lt.0) exit
       opt=char(iopt)
       select case(opt)  
       case('h')
          write(IERR,100)  ! print usage
          call program_exit(1)
       case('H')
          if(action.ne.0) call fexit('only one of -H or -v may be specified')
          action=HEADONLY
       case('l')
          read(optarg,*) lev
       case('v')
          if(action.ne.0) call fexit('only one of -H or -v may be specified')
          var=optarg
          action=PRINTVALS
       case default
          write(IERR,110) opt
          call program_exit(1)
       endselect
    enddo
    if(action.le.0) action=RMSONLY
110 format('unknown option ',a)
!
    call getarg(nargs,fn)
    return
    end subroutine parseline
!-------------------------------------------------------------------------
    subroutine rmsoutput(a,nlev,label)
    character(*),intent(in)::label
    integer,intent(in)::nlev
    real(D),intent(in)::a(LNTD,2,nlev)
!
    integer::k
    real(D)::rms(nlev)
!
    do k=1,nlev
       rms(k)=rmsof(a(:,:,k))
    enddo
    write(LOGFILE,100) label
100 format(a,' rms by level')
    write(LOGFILE,101) rms
101 format(7es11.3)
    return
    end subroutine rmsoutput
!-------------------------------------------------------------------------
    subroutine output(a,lev,label,action)
    character(*),intent(in)::label
    integer,intent(in)::lev,action
    real(D),intent(in)::a(LNTD,2,LEVS)
!
    integer::h,j,k,idx
!
    write(LOGFILE,'(a)') label
    if(action.eq.RMSONLY) then
       do k=merge(lev,1,lev.gt.0),merge(lev,LEVS,lev.gt.0)
          write(LOGFILE,102) 'level',k,rmsof(a(:,:,k))
       enddo
       return
    endif
!
    do h=1,2
       write(LOGFILE,101) 'hemisphere',h
       do k=merge(lev,1,lev.gt.0),merge(lev,LEVS,lev.gt.0)
          write(LOGFILE,102) 'level',k,rmsof(a(:,:,k))
          idx=1
          do j=1,LNTD,7
	  write(LOGFILE,100) idx,a(idx:idx+6,h,k)
	  idx=idx+7
          enddo
       enddo
    enddo
100 format(i4,': ',7es11.3)
101 format(a,1x,i3)
102 format(a,1x,i3,' rms ',es11.3)
    return
    end subroutine output
!-------------------------------------------------------------------------
    subroutine process_file(fn,lev,var,action)
    character(*),intent(in)::fn,var
    integer,intent(in)::lev,action
!
    type(spectral_coeffs)::coeff
    type(sigma_file_header)::header
    integer::k
!
    call read_sigma_file(fn,header,coeff)
    if(action.eq.HEADONLY) then
       call dump_header(fn,LOGFILE,header,VERB_CHATTY)
       return
    else if(action.eq.RMSONLY) then
       call rmsoutput(coeff%gz,1,'gz: geopotential height')
       call rmsoutput(coeff%q,1,'q: surface pressure')
       call rmsoutput(coeff%te,LEVS,'te: temperature')
       call rmsoutput(coeff%di,LEVS,'di: divergence')
       call rmsoutput(coeff%ze,LEVS,'ze: vorticity')
       call rmsoutput(coeff%rq,LEVS,'rq: relative humidity')
       call rmsoutput(coeff%oz,LEVS,'oz: ozone')
       return
    endif
!
    select case(var)
    case('gz')
       call output(coeff%gz,1,'gz: geopotential height',action)
    case('q')
       call output(coeff%q,1,'q: surface pressure',action)
    case('te')
       call output(coeff%te,lev,'te: temperature',action)
    case('di')
       call output(coeff%di,lev,'di: divergence',action)
    case('ze')
       call output(coeff%ze,lev,'ze: vorticity',action)
    case('rq')
       call output(coeff%rq,lev,'rq: relative humidity',action)
    case('oz')
       call output(coeff%oz,lev,'oz: ozone',action)
    case default
       call fexit('unknown variable '//var,1)
    endselect
    return
    end subroutine process_file
!-------------------------------------------------------------------------
    end module sdump
!-------------------------------------------------------------------------
    use sdump
    character(MAX_FILE_NAME_LEN)::fn
    integer::lev,action
    character(3)::var
!
    call parseline(fn,lev,var,action)
    call process_file(fn,lev,var,action)
    call program_exit(0)
    end
