climateClassification.f
********************************************************************************
* This file collects several functions and subroutines for the
* calculation of parameters for climate classification
* TODO: to be extended
5: ********************************************************************************
SUBROUTINE aridIndex (p_x, p_y, p_a, p_temp, p_prec, p_precConv,
c p_curX, p_curY, o_aridIndex)
10: ************************************************************************
* Calculates the arridity indices of a location on the grid for a number
* of years from the monthly average temperatures and precipitations according
* to the definition:
* aridity index = 100 d/n where d is the sum of the monthly differences
15: * between precipitation and pet for water deficient months (i.e. month
* where the pet is higher than the precipitation) and n is the sum of
* the pet for these months
* NOTE: passed in temperatures assumed to be in kelvin!
* NOTE: precipitation should be converted into mm/day (or kg/m^2/day which is
20: * the same) by the convertion factor!
* for THEORY: see paper
* Input:
* - p_x (Integer): grid size of data arrays
* - p_y (Integer): grid size of data arrays
25: * - p_a (Integer): number of years of data
* - p_temp (Real[p_x, p_y, p_a, 12]): monthly temperatures
* - p_prec (Real[p_x, p_y, p_a, 12]): monthly precipitations
* - p_precConv (Real): precipiation convertion factor (to get unit mm/day)
* - p_curX (Integer): location of position for which Aridity incices to be calculated
30: * - p_curY (Integer): location of position for which Aridity incices to be calculated
* Output:
* - o_aridIndex (Real[p_a]): Aridity Indices of every year for the location
************************************************************************
IMPLICIT none
35: INTEGER p_x, p_y, p_a !constants
REAL p_temp(p_x, p_y, p_a, 12), p_prec(p_x, p_y, p_a, 12) !data
REAL p_precConv !precipitation converter
INTEGER p_curX, p_curY !working parameters
REAL o_aridIndex(p_a) !output
40: REAL constI, constA, pet, calcPet, d, n !process params
INTEGER a, i
EXTERNAL pet
do a=1, p_a
45: CALL petConsts(p_x, p_y, p_a,
c p_temp, p_curX, p_curY, a, constI, constA)
d=0
n=0
50: do i=1,12
calcPet = (10./30.) * !adjustment factory, cm to mm, standard month to day
c pet(p_temp(p_curX, p_curY, a, i), constI, constA)
if (calcPet -
55: c p_precConv*p_prec(p_curX, p_curY, a, i) > 0) then
d = d + calcPet -
c p_precConv*p_prec(p_curX, p_curY, a, i)
n = n + calcPet
endif
60: enddo
if (n .EQ. 0) then
o_aridIndex(a) = 0
else
65: o_aridIndex(a) = 100 * d/n
endif
enddo
END
70:
REAL FUNCTION pet (p_temp, p_I, p_A)
************************************************************************
* Calculates the unadjusted potential evapotranspiration of a monthly
75: * mean temperature (NOTE: temperate assumed to be in kelvin)
* PET is returned in cm monthly potential evapotranspiration
* of a standrad month (30 days, each day 12 possible hours of sunshine)
* for THEORY: see paper
* Input:
80: * - p_temp (Real): monthly mean temperature in kelvin
* - p_I (Real): 'I' constant for the calculation (see routine petConsts)
* - p_a (Real): 'a' constant for the calculation (see routine petConsts)
* Return (Real): monthly PET (in cm for 30 day standard month)
************************************************************************
85: IMPLICIT none
REAL p_temp, p_I, p_a, tempK2C
if (tempK2C(p_temp) <= 0 .or. p_I == 0 .or. p_A == 0 ) then
pet = 0.
else
90: pet = 1.6 * ( (10 * tempK2C(p_temp)) / p_I ) ** p_A
endif
END
95: SUBROUTINE petConsts (p_x, p_y, p_a, p_temp,
c p_curX, p_curY, p_curA,
c o_I, o_A)
************************************************************************
* Calculates the pet constants p_I and p_a from a years monthly
100: * average temperatures (NOTE: temperatures assumed to be in Kelvin)
* for THEORY: see paper
* Input:
* - p_x (Integer): grid size of data arrays
* - p_y (Integer): grid size of data arrays
105: * - p_a (Integer): number of years of data
* - p_temp (Real[p_x, p_y, p_a, 12]): monthly temperatures
* - p_curX (Integer): location of position for which PET constants to be calculated
* - p_curY (Integer): location of position for which PET constants to be calculated
* - p_curA (Integer): year for which PET constants to be calculated
110: * Output:
* - o_I (Real): Constant I for the PET calculation
* - o_a (Real): Constant a for the PET calculation
************************************************************************
IMPLICIT none
115: INTEGER p_x, p_y, p_a !constants
REAL p_temp(p_x, p_y, p_a, 12) !data
INTEGER p_curX, p_curY, p_curA !working parameters
REAL o_I, o_A !output
INTEGER i !vars
120: REAL tempK2C !functions
o_I=0.
do i=1,12
if (tempK2C(p_temp(p_curX, p_curY, p_curA, i)) > 0)
125: c o_I = o_I +
c (tempK2C(p_temp(p_curX, p_curY, p_curA, i))/5.)**1.514
enddo
if (o_I == 0 ) then
o_A=0.
130: else
o_A=0.000000675*o_I**3 - 0.0000771*o_I**2 + 0.01792*o_I + 0.49239
endif
RETURN
END
135:
SUBROUTINE yearsDegDays (p_x, p_y, p_a,
c p_temp, p_years, p_curX, p_curY, p_base, o_hdds, o_cdds)
************************************************************************
140: * Calculates the annual heating and cooling degree days for a data set of
* monthly mean temperatures over several years.
* NOTE: all temperatures assumed to be in the same unit (also the base!!)
* I.e. all in Kelvin, all in Celsius or all in Fahrenheit (Kelvin is recommended)
* for THEORY: see paper
145: * Input:
* - p_x (Integer): grid size of data arrays
* - p_y (Integer): grid size of data arrays
* - p_a (Integer): number of years of data
* - p_temp (Real[p_x, p_y, p_a, 12]): monthly temperatures for all years (Kelvin)
150: * - p_years (Integer[p_a]): the years (i.e. e.g. 1975,1976...) that are processed
* - p_curX (Integer): location of position for which Degree Days are to be calculated
* - p_curY (Integer): location of position for which Degree Days are to be calculated
* - p_base (Real): the temperature base (traditionally 18 °C -> 291.15 K)
* Output:
155: * - o_hdds (Real(p_a)): Heating degree days for each year
* - o_cdds (Real(p_a)): Cooling degree days for each year
* NOTE: if the temperatures input was in °C (or K for that matter), the
* resulting degree days will have the unit degrees centigrade/kelvin,
* but if input is in Fahrenheit, the resulting degree days will have
160: * the unit degress fahrenheit (convertion factor from the prior to the
* latter is 9/5 (not *9/5+32!!) because this is a measure of temperature
* DIFFERENZ
************************************************************************
IMPLICIT none
165: INTEGER p_x, p_y, p_a !constants
REAL p_temp(p_x, p_y, p_a, 12) !data
INTEGER p_years(p_a) !years
INTEGER p_curX, p_curY !working parameters
REAL p_base !param
170: REAL o_hdds(p_a), o_cdds(p_a) !output
INTEGER monthDN, daysInMonth !working vars and functoins
REAL stdDev(12), truncPt, truncCoef, monthDegDays
INTEGER i, m
REAL test
175: CALL tempStdDev(p_x, p_y, p_a, p_temp, p_curX, p_curY, stdDev)
do i=1, p_a
o_hdds(i) = 0
o_cdds(i) = 0
180: do m = 1, 12
monthDN = daysInMonth(p_years(i),m)
o_hdds(i) = o_hdds(i) +
c monthDegDays(monthDN,
c p_temp(p_curX, p_curY, i,m), p_base, -1, stdDev(m))
185: o_cdds(i) = o_cdds(i) +
c monthDegDays(monthDN, p_temp(p_curX, p_curY, i,m),
c p_base, 1, stdDev(m))
enddo
enddo
190: END
REAL FUNCTION monthDegDays (p_days, p_temp, p_base, p_direction,
c p_stdDev)
195: ************************************************************************
* Calculates a months normal degree days above or below any base from the
* monthly mean temperature.
* Note: all temperatures need to be in the same unit!! (C, K, F)
* for THEORY: see paper
200: * Input:
* - p_days (Integer): number of days in the particular month
* - p_temp (Real): mean temperature of the month
* - p_base (Real): base temperature for the degree days (18C traditionally)
* - p_direction (Integer) {-1;1}: indicates whether the degree days ABOVE the base
205: * (such as Cooling Degree Days (1)) or BELOW the base (such as Heating Degree
* Days (-1)) are to be calculated
* - p_stdDev (Real): standard deviation of the monthly average temperature
* Return (Real): normal degree days for the month
************************************************************************
210: IMPLICIT none
INTEGER p_days, p_direction
REAL p_temp, p_base, p_stdDev
REAL truncPt, truncCoef
215: if ( p_direction*(p_temp - p_base) > 0 ) then
monthDegDays = p_days * (
c truncCoef(-1*p_direction*
c truncPt(p_temp,p_base,p_days,p_stdDev)
c )*real(p_days)**0.5*p_stdDev +
220: c p_direction*(p_temp - p_base) )
else
monthDegDays = 0
endif
END
225:
SUBROUTINE tempStdDev (p_x, p_y, p_a, p_temp, p_curX, p_curY,
c o_stdDev)
************************************************************************
230: * Calculates the standard deviation for each month from the set of
* temperature data over several years.
* SOURCE (very modified): Numerical Recipes in Fortran 77 Chapter 14.1
* Input:
* - p_x (Integer): grid size of temperature data array
235: * - p_y (Integer): grid size of temperature data array
* - p_a (Integer): number of years of data
* - p_temp (Real[p_x, p_y, p_a, 12]): monthly temperatures for all years
* - p_curX (Integer): location of position for which Standard Devs are to be calculated
* - p_curY (Integer): location of position for which Standard Devs are to be calculated
240: * Output:
* - o_stdDev (Real[12]): Standard deviation of monthly average temperature for each month
************************************************************************
INTEGER p_x, p_y, p_a !constants
REAL p_temp(p_x, p_y, p_a, 12), avg(12), o_stdDev(12) !data and output
245: INTEGER p_curX, p_curY, p_base !working parameters
do m=1, 12
avg(m) = 0
do i=1, p_a
avg(m) = avg(m) + p_temp(p_curX, p_curY, i, m)
250: enddo
avg(m) = avg(m) / REAL(p_a)
o_stdDev(m) = 0
do i=1, p_a
o_stdDev(m) = o_stdDev(m) +
255: c (p_temp(p_curX, p_curY, i,m) - avg(m))**2
enddo
o_stdDev(m) = ( o_stdDev(m) / REAL(p_a - 1) )**0.5
enddo
END
260:
REAL FUNCTION truncPt (p_temp, p_base, p_daysN, p_stdDev)
************************************************************************
* Calculates the truncation point (h or x_0) for a monthly avg.
265: * temperature and standard deviation.
* Note: all temperatures needed in same unit!
* for THEORY: see paper
* Input:
* - p_temp (Real): Average temperature of the month
270: * - p_base (Real): Degree Days base (usually 18C)
* - p_daysN (Integer): number of days in the month
* - p_stdDev (Real): Standard deviation of the month's average temperature
* Return (Real): truncation point
************************************************************************
275: INTEGER p_daysN
REAL p_temp, p_base, p_stdDev
truncPt = (p_base - p_temp)/( real(p_daysN)**0.5 * p_stdDev)
RETURN
END
280:
REAL FUNCTION truncCoef (p_truncPt)
************************************************************************
* Calculates the truncation coefficient l from the truncation point (h or x_0)
* Input:
285: * - p_truncPt (Real): truncation point to be converted
* Return (Real): truncation coefficient
************************************************************************
REAL p_truncPt, e
PARAMETER(e=2.71828)
290: truncCoef = 0.34*e**(-4.7*p_truncPt) - 0.15*e**(-7.8*p_truncPt)
if (p_truncPt .LT. 0) truncCoef = truncCoef + p_truncPt
END
Info Section
Warning: externals (function calls) may not be acurate
calls: petconsts, tempstddev
back to top
f2html v0.3 (C) 1997,98 Beroud Jean-Marc. Fri Aug 11 17:54:58 CEST 2006