#!/bin/sh
#
#  NCLetasubber.sh lonmin lonmax latmin latmax NETCDFINFILE NETCDFOUTFILE
#
# Assumes that the full domain GRIB file has been converted
# into the file NETCDFINFILE  for use by just
# this sub region extraction.
#
# Script for reading the full domain ETA file 
# Use lonmin,lonmax,latmin,latmax to subset the 
# netcdf into a much smaller file
# Output file is NETCDFOUTFILE
# 
# example to run from directory work
# etatime=`MostRecentEta.sh meso-eta_218`
# source NCLwebeta.sh "$etatime" meso-eta_218_
# 
# echo "source NCLwebeta.sh    exports name of ETA netcdf >>"$ETAFILE"<<"
# NCLetasubber.sh -78 -74 36 40 $ETAFILE windsout.nc
#

#cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
#c    NCLetasubber.sh: 
#c    Copyright (C) 2003,  Thomas F. Gross
#c
#c    This program is free software; you can redistribute it and/or modify
#c    it under the terms of the GNU General Public License as published by
#c    the Free Software Foundation; either version 2 of the License, or
#c    (at your option) any later version.
#c
#c    This program is distributed in the hope that it will be useful,
#c    but WITHOUT ANY WARRANTY; without even the implied warranty of
#c    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#c    GNU General Public License for more details.
#c
#c    You should have received a copy of the GNU General Public License
#c    along with this program; if not, write to the Free Software
#c    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
#ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc




lonmin=$1
lonmax=$2
latmin=$3
latmax=$4
NETCDFINPUT=$5
NETCDFOUTFILE=$6

echo "Start NCLetasubber.sh "
echo $lonmin $lonmax $latmin $latmax
echo $NETCDFINPUT
echo $NETCDFOUTFILE


##############################################################################

echo " Start extraction to "$NETCDFOUTFILE
/usr/local/bin/ncl << EOD
;  extractncd.ncl
; extract subgrid from original netcdf file
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
  in = addfile("$NETCDFINPUT","r")
  system("rm $NETCDFOUTFILE")
  out = addfile("$NETCDFOUTFILE","c")
  
lonmin = $lonmin
lonmax = $lonmax
latmin = $latmin
latmax = $latmax
           if(latmin.gt.latmax)then
             print("Error latmin is greater than latmax")
             exit()
           end if
           if(lonmin.gt.lonmax)then
             print("Error lonmin is greater than lonmax")
             exit()
           end if

numdims=dimsizes(dimsizes(in->lon))
print(numdims)
if (numdims.eq.3) then
lon=in->lon(0,:,:)
lat=in->lat(0,:,:)
end if
if (numdims.eq.2) then
lon=in->lon(:,:)
lat=in->lat(:,:)
end if

   r=(lon-lonmin)^2 + (lat-latmin)^2
ijmin=local_min(r,False,0.0)
   r=(lon-lonmax)^2 + (lat-latmax)^2
ijmax=local_min(r,False,0.0)

print(ijmin)
print(ijmax)
; reversed the order of the indices ? 
  iymin=ijmin@xi
  ixmin=ijmin@yi
  iymax=ijmax@xi
  ixmax=ijmax@yi


; Grib names used in the odaas eta.ncf file are translated here:
;sugrd=u-component of wind  m/s 
;svgrd=v-component of wind  m/s
;stmp=Temperature  K
;temp=water temperature   K
;prmsl=pressure reduced to MSL,   Pa
;  Translate those to
; uwind, vwind, temp=stmp C, press  mb

  out->lon=lon(ixmin:ixmax,iymin:iymax)
  out->lat=lat(ixmin:ixmax,iymin:iymax)
  
  out->uwind=in->uwind(:,ixmin:ixmax,iymin:iymax)
  out->vwind=in->vwind(:,ixmin:ixmax,iymin:iymax)
  out->temp=in->temp(:,ixmin:ixmax,iymin:iymax)
  temp=out->temp
  temp=temp-273.15
  out->temp=temp
  out->temp@units="C"
  out->press=in->press(:,ixmin:ixmax,iymin:iymax)
  press=out->press
  press=press/100.
  out->press=press
  out->press@units="mb"
  out->time=in->time
  
  
end

EOD