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