| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051 |
- #/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
|
