!!! (1) Saturation Vapor Pressure !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for estimation of saturation vapor pressure !!! !!! Inputs: !!! !!! temperature (t, deg.C) !!! !!! Output: !!! !!! saturation vapor pressure (esat, hPa) !!! !!! Version: February 09, 2010 (ANU). !!! !!! Reference: !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function esat(t) esatd = t + 243.5 esatn = 17.67 * t esat = 6.112 * exp(esatn/esatd) End Function esat !!! (2) Vapor Pressure !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for estimation of vapor pressure !!! !!! Inputs: !!! !!! temperature (t, deg.C) !!! !!! relative humidity (rh, %) !!! !!! Output: !!! !!! vapor pressure (e, hPa) !!! !!! Version: February 09, 2010 (ANU). !!! !!! Reference: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function e(t,rh) esatd = t + 243.5 esatn = 17.67 * t esat1 = 6.112 * exp(esatn/esatd) e = esat1 * (rh/100.0) End Function e !!! (3) Mixing Ratio !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for estimation of mixing ratio !!! !!! Inputs: !!! !!! temperature (t, deg.C) !!! !!! relative humidity (rh, %) !!! !!! pressure (p, hPa) !!! !!! Output: !!! !!! mixing ratio (xmixrat, g/kg) !!! !!! Version: February 09, 2010 (ANU). !!! !!! (Put 2 versions) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function r(t,rh,p) eps = 0.622 esatd = t + 243.5 esatn = 17.67 * t esat1 = 6.112 * exp(esatn/esatd) e1 = esat1 * (rh/100.0) r = ((eps*e1)/(p - e1))*1000.0 End Function r !!! (4) Saturation Mixing Ratio !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for estimation of saturation mixing ratio !!! !!! Inputs: !!! !!! temperature (t, deg.C) !!! !!! pressure (p, hPa) !!! !!! Output: !!! !!! saturation mixing ratio (satmrat, g/kg) !!! !!! Version: February 09, 2010 (ANU). !!! !!! (Put 2 versions) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function rsat(t,rh,p) eps = 0.62197 esatd = t + 243.5 esatn = 17.67 * t esat1 = 6.112 * exp(esatn/esatd) e1 = esat1 * (rh/100.0) rsat = ((eps*esat1)/(p - esat1))*1000.0 End Function rsat !!! (5) Specific Humidity !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for estimation of specific humidity !!! !!! Inputs: !!! !!! temperature (t, deg.C) !!! !!! relative humidity (rh, %) !!! !!! pressure (p, hPa) !!! !!! Output: !!! !!! specific humidity (sphum, g/kg) !!! !!! Version: February 09, 2010 (DBS). !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function q(t,rh,p) eps = 0.62197 esatd = t + 243.5 esatn = 17.67 * t esat1 = 6.112 * exp(esatn/esatd) e1 = esat1 * (rh/100.0) q = ((eps*e1)/(p - (0.378 *e1)))*1000.0 End Function q !!! (6) Potential Temperature !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for estimation of potential temperature. !!! !!! Inputs: !!! !!! temperature (t, deg.C) !!! !!! pressure (p, hPa) !!! !!! Output: !!! !!! potential temperature (pt, deg.C) !!! !!! Version: February 09, 2010 (DBS). !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function pt(t,p) p0 = 1000.0 tk = t + 273.15 pt = (tk * ((p0/p)**0.286)) - 273.15 End Function pt !!! (7) Virtual Potential Temperature !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for estimation of virtual potential temperature !!! !!! Inputs: !!! !!! temperature (t, deg.C) !!! !!! relative humidity (rh, %) !!! !!! pressure (p, hPa) !!! !!! Output: !!! !!! virtual potential temperature (vpt, deg.C) !!! !!! Version: February 09, 2010 (DBS). !!! !!! Version: February 09, 2010 (ANU). !!! !!! (Stull Version ?) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function vpt(t,rh,p) tk = t + 273.15 p0 = 1000.0 eps = 0.622 esatd = t + 243.5 esatn = 17.67 * t esat1 = 6.112 * exp(esatn/esatd) e1 = esat1 * (rh/100.0) r1 = (eps*e1)/(p - e1) pt1 = tk * ((p0/p)**0.286) vpt = (pt1 * (1.0 + (0.61*r1))) - 273.15 End Function vpt !!! (8) Virtual Temperature !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for estimation of virtual temperature !!! !!! Inputs: !!! !!! temperature (t, deg.C) !!! !!! relative humidity (rh, %) !!! !!! pressure (p, hPa) !!! !!! Output: !!! !!! virtual temperature (vt, deg.C) !!! !!! Version: February 09, 2010 (DBS). !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function vt(t,rh,p) eps = 0.622 e1 = e(t,rh) tden = 1.0 - ((e1/p)*(1.0 - eps)) vt = t/tden End Function vt !!! (9) Dew Point Temperature !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for estimation of dew point temperature !!! !!! Inputs: !!! !!! temperature (t, deg.C) !!! !!! relative humidity (rh, %) !!! !!! Output: !!! !!! dew point temperature (dpt, deg.C) !!! !!! Version: February 09, 2010 (DBS). !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function dpt(t,rh) If (rh == 0.0) rh = 0.01 esatd = t + 243.5 esatn = 17.67 * t esat = 6.112 * exp(esatn/esatd) e = esat * (rh/100.0) dewn = 243.5 * log(e/6.112) dewd = 17.67 - log(e/6.112) dpt = dewn/dewd End Function dpt !!! (10) Lifting Condensation Level Temperature !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for estimation of temperature at LCL !!! !!! Inputs: !!! !!! temperature (t, deg.C) !!! !!! relative humidity (rh, %) !!! !!! pressure (p, hPa) !!! !!! Output: !!! !!! temperature at LCL (tlcl, deg.C) !!! !!! Version: February 09, 2010 (DBS). !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function tlcl(t,rh) tk = t + 273.15 term1d = 1.0/(tk - 55.0) term2d = log(rh/100.0)/2840.0 tlcl = 55.0 - 273.15 + (1.0/(term1d-term2d)) End Function tlcl !!! (11) Equivalent Temperature !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for estimation of equivalent temperature !!! !!! Inputs: !!! !!! temperature (t, deg.C) !!! !!! relative humidity (rh, %) !!! !!! pressure (p, hPa) !!! !!! Output: !!! !!! equivalent temperature (et, deg.C) !!! !!! Version: February 09, 2010 (DBS). !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function et(t,rh) If (t < -40) xlv = 2600.0 If (t > 25) xlv = 2400.0 cpd = 1004.0 tk = t + 273.15 eps = 0.622 esatd = t + 243.5 esatn = 17.67 * t esat = 6.112 * exp(esatn/esatd) e = esat * (rh/100.0) r = (eps*e)/(p - e) et = t + (xlv*r/cpd) End Function et !!! (11) Lifting Condensation Level Pressure !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for estimation of pressure at LCL !!! !!! Inputs: !!! !!! temperature (t, deg.C) !!! !!! relative humidity (rh, %) !!! !!! pressure (p, hPa) !!! !!! Output: !!! !!! temperature at LCL (plcl, hPa) !!! !!! Version: February 09, 2010 (DBS). !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function plcl(t,rh,p) tk = t + 273.15 term1d = 1.0/(tk - 55.0) term2d = log(rh/100.0)/2840.0 tlcl = 55.0 + (1.0/(term1d-term2d)) plcl = p * ((tlcl/tk)**3.5) End Function plcl !!! (12) Equivalent Potential Temperature !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for estimation of equivalent potential temperature !!! !!! Inputs: !!! !!! temperature (t, deg.C) !!! !!! relative humidity (rh, %) !!! !!! pressure (p, hPa) !!! !!! Output: !!! !!! eq.pot.temp. (ept, deg.C) !!! !!! Version: February 09, 2010 (DBS). !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function ept(t,rh,p) ! tk = t + 273.15 ! p0 = 1000.0 ! eps = 0.622 ! esatd = t + 243.5 ! esatn = 17.67 * t ! esat = 6.112 * exp(esatn/esatd) ! e = esat * (rh/100.0) ! r = (eps*e)/(p - e) ! term1d = 1.0/(tk - 55.0) ! term2d = log(rh/100.0)/2840.0 ! tlcl = 55.0 + (1.0/(term1d-term2d)) ! pt = tk * ((p0/p)**0.286) ! ept = (pt * exp(2.67*1000.0*r/tlcl)) - 273.15 tk = t + 273.15 If (t >= 25.0) xlv = 2400.0 If (t < -40.0) xlv = 2600.0 xlv = 2400.0 p0 = 1000.0 cp = 1004.0 rd = 287.0 eps = 0.622 esatd = t + 243.5 esatn = 17.67 * t esat = 6.112 * exp(esatn/esatd) e = esat * (rh/100.0) r = ((eps*e)/(p - e)) ept1 = tk + (xlv*r*1000.0/cp) ept2 = (p0/p)**(rd/cp) ept = (ept1*ept2) - 273.15 End Function ept !!! (13) Saturation Equivalent Potential Temperature !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for estimation of equivalent potential temperature !!! !!! Inputs: !!! !!! temperature (t, deg.C) !!! !!! relative humidity (rh, %) !!! !!! pressure (p, hPa) !!! !!! Output: !!! !!! sat.eq.pot.temp. (vt, deg.C) !!! !!! Version: February 09, 2010 (DBS). !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function sept(t,rh,p) p0 = 1000.0 tk = t + 273.15 eps = 0.622 esatd = t + 243.5 esatn = 17.67 * t esat = 6.112 * exp(esatn/esatd) e = esat * (rh/100.0) rsat = (eps*esat)/(p - esat) pt = tk * ((p0/p)**0.286) sept = (pt * exp(2.67*1000.0*rsat/tk)) - 273.15 End Function sept !!! (14) Pressure to height conversion !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for conversion of pressure to altitude !!! !!! Inputs: !!! !!! surface pressure (p0,hPa) !!! !!! pressure (p, hPa) !!! !!! Output: !!! !!! Altitude (zp,m) !!! !!! Version: February 16, 2010 (ATJ). !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Real Function zp(p0,p) If (p > p0) Write (*, *) 'Value of P must be smaller than P0' term1 = (p/p0)**0.190263 term2 = 288.15/0.00198122 term3 = (1 - term1) * term2 zp = term3 * 0.3048 end function zp !!! (15) julian day to date conversion !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for conversion of Julian day and year to Date !!! !!! Inputs: !!! !!! Julian day (julday, integer) !!! !!! year (iyyy, integer) !!! !!! Output: !!! !!! Date (jul2dd, integer) !!! !!! Version: February 16, 2010 (ATJ). !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Integer Function jul2dd(iyyy,julday) Dimension incmon(13), leapinc(13) Data incmon /0, 31, 59, 90, 120, 151, 181, 212,243, 273, 304, 334, 365/ Data leapinc /0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366/ icheck = mod (iyyy, 4) ! icheck = leap year checking index ! icheck = 0, if the current year is a leap year ! icheck = 1, if the current year is not a leap year if (icheck.ne.0) then do i = 1, 12 j = i + 1 if ( (julday.gt.incmon(i)).and.(julday.le.incmon(j))) jul2dd = julday-incmon(i) Enddo Else Do i = 1, 12 j = i + 1 if ( (julday.gt.leapinc(i)).and.(julday.le.leapinc(j)) ) jul2dd = julday-leapinc(i) Enddo Endif End function jul2dd !!! (16) julian day to month conversion !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for conversion of Julian Day and Year to Month !!! !!! Inputs: !!! !!! Julian day (julday, integer) !!! !!! year (iyyy, integer) !!! !!! Output: !!! !!! Month (jul2mon, integer) !!! !!! Version: February 16, 2010 (ATJ). !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Integer Function jul2mon(iyyy,julday) Dimension incmon(13), leapinc(13) Data incmon /0, 31, 59, 90, 120, 151, 181, 212,243, 273, 304, 334, 365/ Data leapinc /0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366/ icheck = mod (iyyy, 4) ! icheck = leap year checking index ! icheck = 0, if the current year is a leap year ! icheck = 1, if the current year is not a leap year if (icheck.ne.0) then do i = 1, 12 j = i + 1 if ( (julday.gt.incmon(i)).and.(julday.le.incmon(j))) jul2mon=i Enddo Else Do i = 1, 12 j = i + 1 if ( (julday.gt.leapinc(i)).and.(julday.le.leapinc(j)) ) jul2mon=i Enddo Endif End function jul2mon !!! (17) Calendar to Julian day conversion !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for conversion of Calendar to Julian Day Number !!! !!! Inputs: !!! !!! Date (idd, integer) !!! !!! Month (imo, integer) !!! !!! Year (iyy, integer) !!! !!! Output: !!! !!! Julian Day Number (jday, integer) !!! !!! Version: February 16, 2010 (ATJ). !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Integer Function cal2jul(iddd,immm,iyyy) ! iddd = date (Can be anything between 1 to 31) ! immm = month (Can be anything between 1 to 12) ! iyyy = year (Can be literally anything) ! jday = julian day number (Can be anything between 1 to 366) Dimension incmon(13), leapinc(13) Data incmon /0, 31, 59, 90, 120, 151, 181, 212, & & 243, 273, 304, 334, 365/ Data leapinc /0, 31, 60, 91, 121, 152, 182, 213,& & 244, 274, 305, 335, 366/ icheck = mod (iyyy, 4) ! icheck is leap year checking index If (icheck /= 0) julday = iddd + incmon(immm) If (icheck == 0) julday = iddd + leapinc(immm) cal2jul = julday End Function cal2jul !!! (18) wind speed & direction to u conversion !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for conversion of windspeed & direction to u component !!! !!! Inputs: !!! !!! wind speed (ws, m/s) !!! !!! wind direction (wd, degree) !!! !!! Output: !!! !!! u component (u, m/s) !!! !!! Version: February 22, 2010 (ATJ). !!! Real function u( ws, wd) r_d = 180.0/3.14159265359 ! radian to degree conversion d_r = 3.14159265359/180.0 ! degree to radian conversion u = -1. * ws * Sind(wd) End function u !!! (19) wind speed & direction to v conversion !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Function for conversion of windspeed & direction to v component !!! !!! Inputs: !!! !!! wind speed (ws, m/s) !!! !!! wind direction (wd, degree) !!! !!! Output: !!! !!! v component (v, m/s) !!! !!! Version: February 22, 2010 (ATJ). !!! Real function v( ws, wd) r_d = 180.0/3.14159265359 ! radian to degree conversion d_r = 3.14159265359/180.0 ! degree to radian conversion v = -1. * ws * Cosd(wd) End function v