#!/bin/bash
#
# Acquire a set of meteorological forcing data from the North American
# Mesoscale (NAM) WRF model from NCEP NOMADS OPeNDAP server for forcing ROMS
#
# John Wilkin (jwilkin@rutgers.edu) December 2007
# with help from Rich Signell (rsignell@usgs.gov)
#
# This script is step 1 of 3:
#
# In step 1: 
#   Acquire the forcing data in individual files for each day for a user 
#   defined spatial subset from an NCEP NOMADS OPeNDAP server and prepare 
#   them for concatenation
#
#
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# There are multiple NOMADS archives for NAM and NARR models spanning 
# different time ranges at different time resolutions
#
# For NARR access see script:
# step1_narr_get_data_from_nomads
# server=http://nomads.ncdc.noaa.gov:9091/dods/NCEP_NARR_DAILY
#
# For NAM Analysis archive a 6-hour intervals see script:
# step1_get_nam_analysis_from_nomads
# The onlone analysis archive starts 18Z02MAR2004 to recent 
# (approximately 3 month lag) 
# http://nomads.ncdc.noaa.gov:9091/dods/NCEP_NAM_ANALYSIS/Anl_Complete
#
# For NAM 3-hourly 3.5 day forecasts (9 month archive) see script:
# step1_get_nam_9mon3hour_from_nomads
# server=http://nomads.ncdc.noaa.gov:9091/dods/NCEP_NAM
#
# For NAM 1-hourly 1.5 day forecasts see script:
# step1_get_nam_1hour_from_nomads
# server=http://nomad3.ncep.noaa.gov:9090/dods/nam
#
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

echo --------------------------------------------------------------------------
echo
echo Starting script $0
if [ -n "$1" ] 
    then
    archive=$1
    echo for archive option $archive
else
    echo Error in $0   No archive option 1hour, 3hour ... given
    echo Usage is, e.g.:
    echo      step1_get_nam  3hour
    exit
fi

echo 
echo --------------------------------------------------------------------------
echo --------------------------------------------------------------------------
    
# --------------------------------------------------------------------------
# Options the user needs to set to control the spatial/temporal subsetting
# These might be better passed in at the command line, but those
# niceties can come later

# Use to decide whether to overwrite existing files or skip
# =1 causes overwrite
CLOBBER=0

case $archive in 

    1hour )

        # label to append to varname to name output files and directory 
        # in which to Find the extracted subset (from step1) of data 
        # files to be combined 
	label=_nam_1hourly_MAB_and_GoM 
	workdir=./frc_nam_1hourly 

        # Specify the server and elements of the data url:

	# this server was down Dec 20 to April?
	# server=http://nomad3.ncep.noaa.gov:9090/dods/nam

	# this server has only a 2 day archive
	server=http://nomad5.ncep.noaa.gov:9090/dods/nam

	file=nam
	cycle=nam1hr_00z

        # real values extract by coordinate value
	lon_constraint=-82.05,-63.10  # 82 W to 60 W
	lat_constraint=30.90,47.10    # 31 N to 47 N

        # These indices are 1-based (because we use ncks -F below)
	time_constraint=2,25      # hours [1:1:24] for each day
                          # We don't want initial time since precip=0 always

        # days ago counter starts at:
        n=2

	# remove any previous forecast files lest they get merged into 
        # the new archive
	echo rm $workdir/*forecast.nc
	rm $workdir/*forecast.nc
	;;

    3hour )

	label=_nam_3hourly_MAB_and_GoM
	workdir=./frc_nam_3hourly

        # Specify the server and elements of the data url:
	server=http://nomads.ncdc.noaa.gov:9091/dods/NCEP_NAM
	file=nam_218_
	cycle=_0000_fff

        # real values extract by coordinate value
	lon_constraint=-82.00,-60.00  # 82 W to 60 W
	lat_constraint=31.00,47.00    # 31 N to 47 N

        # These indices are 1-based (because we use ncks -F below)
	time_constraint=2,9      # hours [1:1:24] for each day
                          # We don't want initial time since precip=0 always

        # days ago counter starts at:
        n=11 

	# remove any previous forecast files lest they get merged into 
        # the new archive
	echo rm $workdir/*forecast.nc
	rm $workdir/*forecast.nc
	;;

    analysis )

        echo WARNING this option not debugged 
	label=_narr_MAB_and_GoM_2007
	workdir=./frc_nam_analysis 
	;;

    narr )

        echo WARNING this option not debugged 
	label=_nam_analysis_MAB_and_GoM 
	workdir=./frc_narr_2007
	;;

esac

# ----------------------------------------------------------------------
# ----------------------------------------------------------------------

# Record last successful url processed in case we want to go back and 
# get full forecast cycle
last_url=none
last_date=none

# Make use of 'date' commands like
#   date --date="2 days ago"
#   date --date="2 days later"
# with output formatting options to count files
# Here, step through available forecasts for a list of "days ago" 
# relative to the time when the script runs (so this is designed for
# appending to an existing archive) 

# days ago count stops at tomorrow (-1 days ago)
tomorrow=-1

#for n in 4 3 2 1 0 -1 ; do
while [ "$n" -ge "$tomorrow" ] ; do

    n=`expr $n - 1`

    # this format date string is used to build the input url
    yyyymmdd=$(exec date --date="$n days ago" +%Y%m%d)
    yyyymm=$(exec date --date="$n days ago" +%Y%m)
  
    case $archive in
        # url for the relevant cycle on NOMADS
	1hour )
	    url=$server/$file$yyyymmdd/$cycle ;;

	3hour )
	    url=$server/$yyyymm/$yyyymmdd/$file$yyyymmdd$cycle ;;
    esac
    
    echo
    echo ------ $yyyymmdd -------------------------------------
    echo ------ trying data url $url 
    echo
    
    # send global metadata to stdout
    ncks -M $url
    # trap event that the url does not exist (because this day of 
    # month does not exist or archive runs out of data)
    if [ $? -gt 0 ]; then
	echo
	echo _________ URL not found. Probably a non existent date - keep going
	echo 
	continue # move on to next dd
    fi
    
    last_url=$url
    last_date=$yyyymmdd
    echo
    echo Last URL presently $last_url
    echo Last file suffix presently $last_date
    echo
    
    # loop Through variables required
    for varname in ugrd10m vgrd10m tmp2m rh2m pressfc dlwrfsfc dswrfsfc ulwrfsfc uswrfsfc apcpsfc ; do

      # build output file name
      local_file=$workdir/$varname$yyyymmdd.nc

      # Check if file exists before downloading. This way the script
      # can simply be run again if for some reason it dies or the server
      # times out, and it will just skip over the files that were already 
      # downloaded
      
      if [ ! -f $local_file ] || [ $CLOBBER -eq 1 ] ; then
	  echo
	  echo __________  $varname
	  echo
	  
        # extract regional subset for each variable and keep in separate 
        # file (to be concatenated later)
	  echo ncks -F \
	      --dimension lon,$lon_constraint \
	      --dimension lat,$lat_constraint \
	      --dimension time,$time_constraint \
	      --variable $varname \
	      --overwrite $url $local_file
	  ncks -F \
	      --dimension lon,$lon_constraint \
	      --dimension lat,$lat_constraint \
	      --dimension time,$time_constraint \
	      --variable $varname \
	      --overwrite $url $local_file
	  
        # NOTE: the metadata for min,max values of lon,lat are not
        # updated in this ncks extraction step 
	  
        # document the original NCEP variable name and the source
	  ncatted --history -a NCEP_variable,$varname,c,c,$varname $local_file
	  
        # The next three operations convert time to an unlimited dimension 
        # so that the files can be combined later with ncrcat
        # create an unlimited record dimension
	  ncecat --history --overwrite $local_file $local_file 
	  
        # switch time and the record dimension
	  ncpdq -a time,record --history --overwrite $local_file $local_file
	  
        # average to remove the spurious singleton dimension record
	  ncwa -a record --history --overwrite $local_file $local_file
	  
        # special treatment for some variables 
	  case $varname in
	      "apcpsfc" )
            # note analysis interval is 1 hour so we can normalize by this later
		  ncatted --history -a analysis_interval_hours,$varname,c,f,1 $local_file
		  ;;
	  esac
	  
      else	
	  echo
	  echo ------ $local_file already exists ... Skipping

      fi # close file exist check
      
    done # close varname loop
done # close date loop


echo Done processing date list loop 
echo 
echo ______________________________
echo ______________________________


# Now get the full forecast for the lastest available url
# I couldn't think of a clever way to do this, so repeat the entire code 
# segment here
# This should probably be reworked to grab the sequence of cycles up to 
# whatever the latest is, but for now just look for the 12Z cycle if it exists 

echo
echo Repeating extraction at latest url to get forecast $last_url
echo _____________________________________________
echo

# loop through variables required
for varname in ugrd10m vgrd10m tmp2m rh2m pressfc dlwrfsfc dswrfsfc ulwrfsfc uswrfsfc apcpsfc ; do
    
    case $archive in

	1hour )
	    forecast_url=${last_url/00z/12z/}
            # build output file names
	    local_file=$workdir/$varname$last_date.nc
	    nowcast=$workdir/$varname$last_date"_00Znowcast.nc"
	    forecast=$workdir/$varname$last_date"_12Zforecast.nc"
	    ftime_constraint=2,37
	    ;;

	3hour )
	    forecast_url=$last_url
            # build output file names
	    local_file=$workdir/$varname$last_date.nc
	    forecast=$workdir/$varname$last_date"_forecast.nc"
	    ftime_constraint=2,29
	    ;;
    esac

    ncks -M $forecast_url
    # trap event that the 12z cycle url does not exist
    if [ $? -gt 0 ]; then
	echo
	echo ---------- URL $forecast_url not found. 
        echo ---------- 12Z cycle probably not available yet
	echo 
	break
    fi

    echo
    echo ----------- Building $forecast 
    echo
    
    case $archive in

	1hour )
   
            # Trim last file keeping only first 12 hours
	    echo ncks -F --dimension time,1,12 $local_file $nowcast
	    echo
	    ncks -F --dimension time,1,12 $local_file $nowcast
	    ;;

	3hour )

            # remove the current file, to be replaced with complete forecast cycle
	    rm $local_file
	    ;;

    esac

    # Extract forecast with appropriate time constraint
    echo ncks -F --dimension lon,$lon_constraint --dimension lat,$lat_constraint \
	--dimension time,$ftime_constraint \
	--variable $varname \
	--overwrite $forecast_url $forecast

    ncks -F --dimension lon,$lon_constraint --dimension lat,$lat_constraint \
	--dimension time,$ftime_constraint \
	--variable $varname \
	--overwrite $forecast_url $forecast

    ncatted --history -a NCEP_variable,$varname,c,c,$varname $forecast
    ncecat --history --overwrite $forecast $forecast 
    ncpdq -a time,record --history --overwrite $forecast $forecast
    ncwa -a record --history --overwrite $forecast $forecast
    
    # special treatment for some variables 
    case $varname in
	"apcpsfc" )
            # note analysis interval is 1 hour so we can normalize by this later
	    ncatted --history -a analysis_interval_hours,$varname,c,f,1 $forecast
	    ;;
    esac

    case $archive in

	1hour )

            # concatenate the nowcast and forecast 
	    echo ncrcat --overwrite $nowcast $forecast $workdir/$varname$last_date"_latest.nc" 
	    echo
	    ncrcat --overwrite $nowcast $forecast $workdir/$varname$last_date"_latest.nc" 
	    
            # remove temporary files so that ncrcat on directory has desired outcome
	    rm $local_file
	    rm $nowcast
	    rm $forecast
	    ;;

	3hour )

            # do nothing because forecast file is what we want
	    ;;

    esac

done

echo DONE with step1 latest url $last_url
echo job step1 completed at $(exec date) 
echo