!  NAT2OBS - convert nature-run, synthetic observations to the raw_obs
!  file format used in newer versions of the data assimilation code.
!  This program can handle the obsgrid format in /raid/szunyogh/observations.
!  Note: the input files are big-endian.
!
    use grid_types
    use grid_obs,only:raw_obs
    use mrf_ids
    use ferror
    use io
    implicit none
    character(MAX_FILE_NAME_LEN)::infn,outfn
    integer::j,k,nrandom,verblev,idate,ioerr,nobs
    type(patch_index)::loc(LONF,LATG),loc1(LONF*LATG)
    equivalence(loc,loc1)
    real(D)::obsgrid(LONFD,2,8,LATG2,LEVS)
    integer::idx(LONF*LATG)=(/(j,j=1,LONF*LATG)/)
!
    call parseline(infn,outfn,idate,nrandom,verblev)
    open(unit=INSWAP,file=infn,form='unformatted',status='old',iostat=ioerr)
    if(ioerr.ne.0) call fexit(trim(infn)//': cannot open')
    open(unit=OUTNOSWAP,file=outfn,form='unformatted',status='replace', &
       iostat=ioerr)
    if(ioerr.ne.0) call fexit(trim(outfn)//': cannot create')
!
    read(INSWAP,err=911) obsgrid
    write(OUTNOSWAP) idate
    do k=1,LATG
       do j=1,LONF
          loc(j,k)=patch_index(j,k,1)
       enddo
    enddo
    nobs=nrandom
    if(nrandom.lt.LONF*LATG) call randperm(idx,LONF*LATG)
    if(verblev.gt.0) then
       write(IERR,100) infn,outfn,nrandom
       write(IERR,101) minval(loc1(idx(1:nrandom))%y), &
          maxval(loc1(idx(1:nrandom))%y)
    endif
100 format('nat2obs: input: ',a/'output: ',a/'nobs: ',i0)
101 format('nat2obs: check latitudes: min: ',i0,' max: ',i0)
!
    do j=1,nobs
       call process(verblev,obsgrid,loc1(idx(j)),Q)  ! surface pressure
       call process(verblev,obsgrid,loc1(idx(j)),TE)  ! temperature
       call process(verblev,obsgrid,loc1(idx(j)),UG)  ! uwind
       call process(verblev,obsgrid,loc1(idx(j)),VG)  ! vwind
    enddo
    stop
!
911 continue
    call finfo('i/o error; on Lahey, export FORT90L="-Wl,-T1,-T2"')
    end
!------------------------------------------------------------------------------
    subroutine parseline(infn,outfn,idate,nrandom,verblev)
!  PARSELINE - parse the command line.
!
    use resolution
    use io
    use command_line_support
    implicit none
    character(*),intent(out)::infn,outfn
    integer,intent(out)::nrandom,idate,verblev
!
!  Local variables
!
    character(*),parameter::OPTLIST='d:hi:o:n:v:'
    character(MAX_FILE_NAME_LEN)::optarg
    integer::nargs,iopt,optind
    integer,external::iargc
!
!  NARGS is the number of command-line arguments EXCLUDING the program name
!
    nargs=iargc()
    infn=' '  ! no defaults; -i and -o must be specified
    outfn=' ' 
    idate=0  ! -d should be specified
    verblev=0  ! 0 => silent, increasing chattiness with level
    nrandom=LONF*LATG  ! full observing network by default
    do 
       call getopt(iopt,optarg,optind,OPTLIST)
       if(iopt.lt.0) exit
       select case(char(iopt))
       case('d')
          read(optarg,*) idate   !  in yyyymmddhh format
       case('i')  ! list of variables to be assimilated
          infn=optarg
       case('o')
          outfn=optarg
       case('h')
          call usage
       case('n')    ! N random locations
          read(optarg,*) nrandom
       case('v')
          read(optarg,*) verblev
       case default
          call program_exit(1)     ! getopt prints message
       endselect
    enddo
!
!  Get the required argument if no namelist input is given
!
    if(len_trim(infn).eq.0.or.len_trim(outfn).eq.0) then
       call usage
       call program_exit(1)
    endif
!
    contains
       subroutine usage
       write(IERR,103)
103    format('nat2obs -i infile -o outfile [-n #random] [-v verblev]'/ &
          '[-d date (yyyymmddhh format)]'/ &
          '-i and -o must be specified.'/ &
          'If -n is omitted or value is 0 then all observations are copied;'/ &
          'otherwise N randomly-chosen (x,y) points are chosen'/ &
          'and data for PTUV are output at every level')
       call program_exit(1)
       end subroutine usage
    end subroutine parseline
!------------------------------------------------------------------------------
    subroutine process(verblev,obsgrid,loc,var)
!  PROCESS - construct and output an observation in raw_obs format from
!  the synthetic observations on the model grid for variables of type VAR.
!  We output values at all vertical levels for all dynamical variables
!  for the (x,y) coordinate given in LOC.
!
    use grid_types
    use grid_obs,only:raw_obs
    use resolution
    use mrf_ids
    use ferror
    use io
    implicit none
    integer,intent(in)::verblev,var
    real(D),intent(in)::obsgrid(LONFD,2,8,LATG2,LEVS)
    type(patch_index),intent(in)::loc
!
    integer::ix,ih,iy,iz,nobs,ioerr,lat
    integer,parameter::NH=1,SH=2  ! Northern & Southern hemispheres
    type(raw_obs)::robs
!
!  GFS/obsgrid indexing:  IY=1 <=> pole, IY=LATG2 <=> equator
!  Our physical grid indexing:  Y=1 <=> South Pole, Y=LATG <=> North Pole
!
    ix=loc%x
    if(loc%y.le.LATG2) then  ! our Southern Hemisphere
       ih=SH
       iy=loc%y
    else                     ! our Northern Hemisphere
       ih=NH
       iy=LATG-loc%y+1
    endif
!
!  Dummy values for irrelevant fields; fake values for others
!
    robs%typecode=0
    robs%lev=0
    robs%time_diff=0.0
    robs%p_lev=0.0  ! special flag value for synthetic observations
    robs%qcm=1.0  ! high quality obs
    robs%cat=0
    robs%prss=1000.0
    robs%grp='NATURE'
    robs%sid='NATU'
!
    select case(var)
    case(Q)
       robs%field='P'
    case(TE)
       robs%field='T'
    case(UG)
       robs%field='U'
    case(VG)
       robs%field='V'
    case default
       call fabort('unknown variable type')
    endselect
    robs%x_grd=loc%x  ! same internal (x,y) indexing
    robs%y_grd=loc%y
    !
    !  LONSPERLAT is indexed from 1 to LATG2
    !  RLATS is indexed from 1 to LATG
    !  See ~/support/resolution.f90 for details.
    !
    robs%lambda=(360.0/LONSPERLAT(iy))*(ix-1) ! 0 degrees <=> ix=1
    robs%phi=RLATS(loc%y)
    if(var.eq.Q) then  ! surface pressure on level 1 only
       robs%sigma=1.0
       robs%value=obsgrid(ix,ih,Q,iy,1) ! log(dPa)
       robs%sdev=1.0e-04  ! dummy value so grid_obs will substitute 4e-04
       write(OUTNOSWAP) robs
    else
       robs%sdev=merge(1.0,1.1,robs%field.eq.'T') ! T=>1.0 K, wind=> 1.1 m/s
       do iz=1,LEVS
          robs%sigma=iz
          robs%value=obsgrid(ix,ih,var,iy,iz)
          write(OUTNOSWAP) robs
       enddo
    endif
    return
    end subroutine process
!------------------------------------------------------------------------------
    subroutine randperm(ix,n)
!  RANDPERM - create a random permutation of the integers 1 to N and place
!  into IX.
!
    implicit none
    integer,intent(in)::n
    integer,intent(out)::ix(n)
!
    integer::j,idx,itemp
    real::r
!
    !!call random_seed  - always use the same observing network
    do j=1,n
       call random_number(r)
       idx=min(int(r*n)+1,n)
       itemp=ix(n)
       ix(n)=ix(idx)
       ix(idx)=itemp
    enddo
    return
    end subroutine randperm
