#!/bin/bash
#
# Acquire a set of meteorological forcing data from the North American
#   Regional Reanalysis (NARR) NCEP NOMADS OPeNDAP server for forcing ROMS
#
# John Wilkin
# 13 August 2007
#
# This script is step 2 of 2:
#
# In step 2:
#
#   Concatenate the individual daily files extracted by step 1
#   Convert units if required
#   Change names to ROMS defaults


# BUGS/FEATURES: The very last step in this script is a hack to reformat the
# output so ROMS can read it. For some reason ROMS' get_2dfld complains about
# the file this script creates, but merely dumping the file with ncdump, and 
# reconstituting it with ncgen, solves the problem. Greater netcdf minds than 
# mine (i.e. John Evans) have not been able to solve this problem

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

# Options the user sets -------------------------------------------------------

# The year string is passed in from the command line
echo Starting script $0
if [ -n "$1" ]
    then
    yyyy=$1
    echo for year $yyyy
else
    echo Error in $0   No year to process requested
    echo Usage is, e.g.:
    echo      step1_narr_cat_change_units_names 2008 
    exit
fi

# label to append to varname to name output files
# and baseyear used in shifting time units (see notes below on how 
# to control the time units conversion)  
label=_narr_MAB_and_GoM_$yyyy
baseyear=2006

label=_narr_NENA_$yyyy
baseyear=2001

label=_narr_Delaware_$yyyy
baseyear=1968


# directory in which to find the extracted narr subset data files to
# be combined 
workdir=./frc_narr_$yyyy


# IMPORTANT: You must search for the string "days since" in the code below and 
# change the time units base date. You'll have to use udunits do the math on 
# the days offset. Follow the instructions are given below.

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


# loop through all variables required and concatenate daily files to a 
# single file

echo Concatenating all files in $workdir for variables
for varname in ugrd10m vgrd10m tmp2m rh2m pressfc dlwrfsfc dswrfsfc ulwrfsfc uswrfsfc apcp ; do
  local_file=$varname$label.nc
  # don't create history global attribute for this operation or you get a very long 
  # string (so many files) that can upset some netcdf readers
  echo " " $varname
  ncrcat $workdir/$varname*.nc --history --overwrite $local_file
done


echo
echo It is safe to ignore the same file warnings that follow
echo Changing units or computing net for

# some variables get special treatment ----------------------------------------
# -----------------------------------------------------------------------------  
# (1) convert air temperature in Kelvin to Celsius for ROMS default varinfo.dat
#     TC = TK - 252.15
varname=tmp2m
echo "  "$varname
local_file=$varname$label.nc

# these operations require output to a scratch file which gets renamed back to 
# the original when we're done

scratch=scr.$local_file
ncap -s 'tmp2m=tmp2m-272.15' --overwrite $local_file $scratch
ncatted -a units,tmp2m,c,c,Celsius $scratch
ncatted -a long_name,tmp2m,c,c,"anl 2 m temp. [c]" $scratch
ncatted -a history,tmp2m,c,c,"converted from Kelvin " $scratch
mv $scratch $local_file


# (2) convert air pressure in Pa to millibar for ROMS default varinfo.dat -----
#     millibar = Pa * 0.01
varname=pressfc
echo "  "$varname
local_file=$varname$label.nc

# these operations require output to a scratch file which gets renamed back to 
# the original when we're done

scratch=scr.$local_file
ncap -s 'pressfc=pressfc/100' --overwrite $local_file $scratch
ncatted -a units,pressfc,c,c,millibar $scratch
ncatted -a long_name,pressfc,c,c,"anl surface pressure [millibar] " $scratch
ncatted -a history,pressfc,c,c,"converted from Pa " $scratch
mv $scratch $local_file


# (3) rain is accumulated kg m^-2 over 3 hours so we need to convert 
# to kg m^-2 s^-1 for ROMS default varinfo.dat --------------------------------

varname=apcp
echo "  "$varname
local_file=$varname$label.nc

scratch=scr.$local_file
ncap -s 'apcp=apcp/10800' --overwrite $local_file $scratch
ncatted -a units,apcp,c,c,"kg/m2/s" $scratch
ncatted -a \
    long_name,apcp,c,c,"** 0-3 hr ave precipitation rate [kg/m^2/s^1] " $scratch
ncatted -a \
    history,apcp,c,c,"computed from 3 hr accumulated precip / 10800 " $scratch
mv $scratch $local_file


# (4) subtract upward shortwave from downward to get net ----------------------

varname=nswrfsfc # net shortwave will be named this in the output file
echo "  "$varname
local_file=$varname$label.nc

# ncdiff requires the variables have the same name so we have to rename first
local_fileup=uswrfsfc$label.nc
ncrename -v uswrfsfc,$varname --overwrite $local_fileup $local_fileup
local_filedown=dswrfsfc$label.nc
ncrename -v dswrfsfc,$varname --overwrite $local_filedown $local_filedown
ncdiff $local_filedown $local_fileup $local_file
ncatted -a \
    long_name,nswrfsfc,m,c,"** 0-3 hr ave net shortwave radiation flux [w/m^2]" \
    $local_file 
ncatted -a \
  history,nswrfsfc,c,c,"computed from dswrfsfc minus uswrfsfc" $local_file


# (5) subtract upward longwave from downward to get net ----------------------

varname=nlwrfsfc # net longwave will be named this in the output file
echo "  "$varname
local_file=$varname$label.nc

# will use ncdiff to subtract variables, but ncdiff requires the variables 
# have the same name in the two file, so we have to rename first
local_fileup=ulwrfsfc$label.nc
ncrename -v ulwrfsfc,$varname --overwrite $local_fileup $local_fileup
local_filedown=dlwrfsfc$label.nc
ncrename -v dlwrfsfc,$varname --overwrite $local_filedown $local_filedown
ncdiff $local_filedown $local_fileup $local_file
ncatted -a \
    long_name,nlwrfsfc,m,c,"** 0-3 hr ave net longwave radiation flux [w/m^2]" \
    $local_file 
ncatted -a \
    history,nlwrfsfc,c,c,"computed from dlwrfsfc minus ulwrfsfc" $local_file

# rename downward longwave back to its original name because we need this in
# the event that we use the ROMS options to input downward longwave only 
# (in variable lwrad_down ; activated with #define LONGWAVE_OUT)

ncrename -v nlwrfsfc,dlwrfsfc --overwrite $local_filedown $local_filedown


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

# Final steps are to 
#
# (1) Shift the time base date from NOMADS days since 1-1-1 to
# something more recent.
#
# (2) Add the 'time' and 'coordinates' attributes that inform ROMS getdata
# what the appropriate time variable and spatial coordinates are in order
# to process the data correctly in get_*fld
#
# (3) Rename variables and files to follow the ROMS varinfo.dat defaults
# This is not essential because you can midify varinfo.dat to point ROMS
# to the names used by NCEP, but this way we can run with the default names
# and we already converted to default units anyway
#
# To figure out the time shift use the unix command 'udunits' command to work out 
# the shift
# On the Rutgers UNIX hosts this is /opt/compiled/gnu/opendap/bin/udunits
# DO NOT do the math in Matlab - you will end up off by 2 days.
#
# At the unix shell type:
#
#    udunits
#
# And respond by typing appropriate strings including the days since part
#
#    You have: days since 1-1-1 00:00:0.0
#    You want: days since 2006-01-01 00:00:0.0
#
# You will be given the repsonse:
#
#    <days since 2006-1-1 00:00> = <days since 1-1-1 00:00:0.0>*1 - 732313
# 
# The number at the end is what you need below.


echo Make attribute changes for 

for varname in ugrd10m vgrd10m tmp2m rh2m pressfc dlwrfsfc nlwrfsfc nswrfsfc apcp ; do

  echo "  "$varname

  local_file=$varname$label.nc

  # Shift time and update the time units attribute ----------------------------
  # This loses the accompanying time metadata inherited from the NOMADS files 
  # but it didn't get updated during the extraction and concatenation so it
  # was wrong anyway

  scratch=scr.$local_file

  # (Ideally, the math in the next line could be done by computing the time 
  # offset for a given time units base date, but I don't know how to do 
  # that and don't have the inclination to work it out right now).

  case $baseyear in
      2006 )
	  # ESPreSSO forcing files
	  ncap -s 'time=time-732313' --overwrite $local_file $scratch
	  ncatted -a units,time,c,c,"days since 2006-01-01 00:00:0.0" $scratch
	  ;;
      2001 )
	  # NENA forcing files
	  ncap -s 'time=time-730487' --overwrite $local_file $scratch
	  ncatted -a units,time,c,c,"days since 2001-01-01 00:00:0.0" $scratch
	  ;;
      1968 )
	  # Delaware forcing files
	  ncap -s 'time=time-718576' --overwrite $local_file $scratch
	  ncatted -a units,time,c,c,"days since 1968-05-23 00:00:0.0" $scratch
	  ;;
  esac

  ncatted -a long_name,time,c,c,"time" $scratch
  ncatted -a grads_step,time,c,c,"3hr" $scratch
  ncatted -a minimum,time,c,c,"   " $scratch
  ncatted -a maximum,time,c,c,"   " $scratch
  mv $scratch $local_file

  # Add the time and coordinates attributes
  # The coordaintes will allow ROMS to do the spatial interpolation from the 
  # native NOMADS lon,lat to the ROMS grid, without the user needing to 
  # create a set of forcing files for the particular grid. This means we can 
  # use these files created for MABGOM for all the simulations on inner 
  # nested spatial domains (LaTTE, ESPreSSO, SW06 ...)

  ncatted -a time,$varname,c,c,time $local_file
  ncatted -a coordinates,$varname,c,c,"lon lat" $local_file


  # Who to blame
  ncatted --history -a source,global,c,c,"Generated by script $0 by John Wilkin (jwilkin@rutgers.edu)" $local_file


  # Rename to the default names expected in ROMS varinfo.dat file -------------
  case $varname in 
    ugrd10m   ) newvarname=Uwind ;;
    vgrd10m   ) newvarname=Vwind ;;
    tmp2m     ) newvarname=Tair ;;
    rh2m      ) newvarname=Qair ;;
    pressfc   ) newvarname=Pair ;;
    dlwrfsfc  ) newvarname=lwrad_down ;;
    nlwrfsfc  ) newvarname=lwrad ;;
    nswrfsfc  ) newvarname=swrad ;;
    apcp      ) newvarname=rain ;;
  esac

  ncrename -v $varname,$newvarname $local_file
  final_file=$newvarname$label.nc
  mv $local_file $final_file

  # For some reason the file at this point is not always readable by
  # ROMS get_2dfld.F, yet the files appear to be perfectly acceptable
  # to other netcdf reading tools (Matlab, ncBrowse). 
  #
  # GET_2DFLD complains: unable to find requested variable: time
  #                      in input NetCDF file: ...
  # but time is most certainly there.  
  #
  # Dumping the ENTIRE file with ncdump (yes, all the data) and then
  # reconstituting the file with ncgen -b gives an outcome readable by
  # ROMS.  Go figure.
  #
  # You can try commenting out this final step and see if the files are
  # OK as is. If you figure out what is going wrong with the nco tools
  # operations please let me know.  I have to do this when I build the
  # netcdf files on queequeg.marine.rutgers.edu and run ROMS on
  # moby.marine.rutgers.edu

  echo "   ncdump/ncgen on "$final_file" to solve get_2dfld read problem"
  scratch=scratch_$label.cdl
  ncdump $final_file > $scratch
  ncgen -b $scratch
  /bin/rm $scratch

done