#/bin/ksh
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# This script computes the specific humidity from the dew point using :
#  
# 1) the approximation e<<P which gives q=0.622(e/P)=0.622(es(Td)/P)
#
# 2) And one of the three formulas :
#
# - the August-Roche-Magnus formula : es=C1 x exp( (A1 x T)/(B1 + T) )
#         A1=17.625, B1=243.04°C, C1=6.1094 hPa
#    whose error <= 0.4% in the range -40°C < T < 50°C. 
#   from Alduchov and Eskridge (1996) Ueber die Berechnung der Expansivkraft 
#        des Wasserdunstes. Ann. Phys. Chem., 13, 122–137.
#
# - the Lawrence formula : es=C2 x exp( -L/(R*(T+273.15)) )
#         C2=2.53e9 hPa, L=2.501e6 J.kg-1, R=461.5 J.K-1.kg-1 
#   from Lawrence (2005), BAMS, 86, 225-233
# 
# - the Tetens formula : es= (D1 + E1*P) x C1 x exp((A1 x T)/(B1 + T) )
#         D1=1.0007 hPa, E1=3.46e-6, C1=6.1121 hPa, B1=240.97°C, A1=17.502      
#   from  Buck, 181, JAM, 20, 1527-1532
#
# Review  : Gibbins (1990) Ann. Geophys., 8, 859–885. 
#
# History : Virginie Guemas - Initial version                         - 2012 
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#                                Arguments
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
year0='2012'
yearf='2014'
equation='Lawrence'
fileD='/cfu/scratch/vguemas/ERAint/d2m_eraint'  # Dew point
fileP='/cfu/scratch/vguemas/ERAint/sp_eraint'   # Surface pressure
fileout='q2_eraint'
varD="d2m"
varP='sp'
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

for ((year=${year0};year<=${yearf};year++)) ; do
  ncap -s ${varD}"=double("${varD}")" ${fileD}_${year}.nc tmp1_${year}.nc
  cdo subc,273.15 tmp1_${year}.nc tmp_${year}.nc
  ncap -s ${varP}"=double("${varP}")" ${fileP}_${year}.nc toto_${year}.nc
  case $equation in 
    'Magnus') cdo expr,'esd=6.1094*exp((17.625*'${varD}')/(243.04+'${varD}'));' tmp_${year}.nc tmp2_${year}.nc ;;
    'Lawrence') cdo expr,'esd=2.53e9*exp(((d2m+273.15)^(-1))*(-2.501e6/461.5));' tmp_${year}.nc  tmp2_${year}.nc ;;
    'Tetens') cdo expr,'esd=(1.0007+3.46e-6*1013)*6.1121*exp((17.502*d2m)/(240.97+d2m));' tmp_${year}.nc  tmp2_${year}.nc ;;
  esac
  ncks -A -v sp toto_${year}.nc tmp2_${year}.nc
  cdo expr,'q2=0.622*esd/(sp/100)' tmp2_${year}.nc ${fileout}_${year}.nc
  rm -f tmp_${year}.nc tmp1_${year}.nc tmp2_${year}.nc toto_${year}.nc
done