analogueSearch.f

      ********************************************************************************
      *     This file collects search routines for analogues.
      *     TODO: to be extended
      ********************************************************************************
   5: 
      
            SUBROUTINE translateCoordinates (p_x, p_y, 
           c p_lats, p_longs, p_lat, p_lon, o_x, o_y )
      ************************************************************************
  10: * Returns the location x,y on the grid that is closest to the passed in
      * longitude and lattitude.
      * Input:
      *  - p_x (Integer): grid size of coordinate arrays
      *  - p_y (Integer): grid size of coordinate arrays
  15: *  - p_lats (Real[p_x, p_y]): latitudes
      *  - p_longs (Real[p_x, p_y]): longitudes
      *  - p_lat (Real): Coordinates to be searched for
      *  - p_lon (Real): Coordinates to be searched for
      * Output:
  20: *  - o_x (Integer): grid entry of closest point on grid
      *  - o_y (Integer): grid entry of closest point on grid
      ************************************************************************
            IMPLICIT NONE 
             
  25:       INTEGER p_x, p_y ! grid size x/y
            REAL p_lats(p_x,p_y), p_longs(p_x,p_y) ! longitute/lattitue data
            REAL p_lat, p_lon ! Parameters
            INTEGER o_x, o_y ! Output
            INTEGER x, y ! Vars
  30:       REAL dist, minDist ! Vars
            
            minDist=-1
            do x=1, p_x
              do y=1, p_y
  35:            dist = (p_lon-p_longs(x,y))**2+(p_lat-p_lats(x,y))**2
                 if ( (dist .lt. minDist) .or. (minDist .eq. -1) ) then
                    minDist = dist
                    o_x = x
                    o_y = y
  40:            endif
              enddo
            enddo
            
            print *, "Desired coordinates:",p_lat," /",p_lon
  45:       print *, "Closest grid coordinates:",
           c p_lats(o_x, o_y)," /",p_longs(o_x, o_y)
            print *, "Corresponding grid point:", o_x," ;",o_y
            
            END
  50: 
      
      
            SUBROUTINE avdFindBestAlg (p_x, p_y, 
           c p_longs, p_lats,
  55:      c p_curTemp, p_curPrec, 
           c p_a, p_curYears, 
           c p_cityAridIndex, p_cityHDDs, p_cityCDDs, 
           c p_ddBase, p_decide, 
           c o_bestX, o_bestY, o_d, o_dists, o_probs)
  60: ************************************************************************
      * Analogue search that runs depending on the input either a 3-dimensional
      * K-S test (with Aridity Index, HDD and CDD) to find the best analogue for
      * the passed in city's data, or only a 2-dimensional K-S test (with Aridity 
      * Index and HDD _or_ CDD).
  65: * A preselection of a maximum of 25 best analogues (with a p-value greater
      * than 0.5) are then compared taking their local neighborhood (0.5 degrees
      * latitude/longitude) into consideration and the location of the best
      * analogue is returned along with the K-S distances and p-values of each
      * grid point.
  70: * for more THEORY: see paper
      *
      * Note: temperature data assumed to be in Kelvin, precipitation data assumed
      * to be in mm/day !!
      *
  75: * Input:
      *  - p_x (Integer): grid size of data arrays
      *  - p_y (Integer): grid size of data arrays
      *  - p_longs (Real[p_x, p_y]): longitudes of data grid
      *  - p_lats (Real[p_x, p_y]): latitudes of data grid
  80: *  - p_curTemp (Real[p_x, p_y, p_a, 12]): temperature data that is to be compared with the city
      *  - p_curPrec (Real[p_x, p_y, p_a, 12]): precipitation data that is to be compared with the city
      *  - p_a (Integer): number of years of data
      *  - p_curYears (Integer[p_a]): the years that data is for (real year number likes 1984,1985,etc.)
      *  - p_cityAridIndex (Real[p_a]): annual aridity indices of the city for which to search an alg
  85: *  - p_cityHDDs (Real [p_a]): annual heating degree days of the city for which to search an alg
      *  - p_cityCDDs (Real [p_a]): annual cooling degree days of the city for which to search an alg
      *  - p_ddBase (Real): temperature base (also in Kelvin) for degree days
      *  - p_decide (Integer): which test to run:
      *       1 - 3D K-S test with Aridity Index, HDD and CDD
  90: *       2 - 2D K-S test with Aridity Index and HDD
      *       3 - 2D K-S test with Aridity Index and CDD
      * Output:
      *  - o_bestX (Integer): location on grid of the best analogue found
      *  - o_bestY (Integer): location on grid of the best analogue found
  95: *  - o_d (Real): average K-S distance (taking neighbors into consideration) of winning analogue
      *       [0;1]: analogue found with this distance
      *       -1: no good analogue found (do not use o_bestX and o_bestY either!)
      *  - o_dists (Real[p_x, p_y]): computed K-S distances of all grid points
      *  - o_dists (Real[p_x, p_y]): computed p-values of all grid points
 100: ************************************************************************
            IMPLICIT NONE
      
            ! Constants      
            INTEGER p_x, p_y, p_a ! grid size x/y and number of years p_a
 105:       
            ! Data
            REAL p_longs(p_x, p_y), p_lats(p_x, p_y)
            REAL p_curTemp(p_x,p_y,p_a,12)
            REAL p_curPrec(p_x,p_y,p_a,12)
 110:       INTEGER p_curYears(p_a)
            
            ! Parameters
            INTEGER p_decide !which test to use
            REAL p_cityAridIndex(p_a)
 115:       REAL p_cityHDDs(p_a), p_cityCDDs(p_a)
            REAL p_cityTemps(12), p_cityPrecs(12) ! average city temperatues
            REAL p_ddBase ! temperature base for degree days
            
            !Output
 120:       INTEGER o_bestX, o_bestY
            REAL o_d
            REAL o_dists(p_x, p_y)
            REAL o_probs(p_x, p_y)
            
 125:       ! Vars
            INTEGER i,j,x,y,a !looping integer
            REAL algAridIndex(p_a)
            REAL algHDDs(p_a), algCDDs(p_a)
            REAL d, prob, rr !kol-smir results
 130:       
            INTEGER best, iBest
            PARAMETER(best=25)
            INTEGER bestX(best), bestY(best), worstBest
            REAL bestD(best)
 135:       REAL nD, clusterDistance
            PARAMETER(nD = 0.5)
            
            do i=1,best
                  bestD(i) = 1000.
 140:       enddo
            worstBest=1
             
            o_d = -1
            do x=1, p_x
 145:         do y=1, p_y
                 
                 CALL aridIndex (p_x, p_y, p_a, p_curTemp, p_curPrec, 1., 
           c       x, y, algAridIndex)
            
 150:            CALL yearsDegDays (p_x, p_y, p_a, 
           c       p_curTemp, p_curYears, x, y, 
           c       p_ddBase, algHDDs, algCDDs)
                 
                 if (p_decide == 1) then !3d test   
 155:            CALL ks3d2s(p_cityAridIndex, p_cityHDDs, p_cityCDDs, p_a,
           c            algAridIndex, algHDDs, algCDDs, p_a, d, prob, rr)         
                 elseif (p_decide == 2) then !2d with aridity index + HDD
                 CALL ks2d2s(p_cityAridIndex, p_cityHDDs, p_a,
           c            algAridIndex, algHDDs, p_a, d,prob, rr) 
 160:            elseif (p_decide == 3) then !2d with aridity index + CDD
                 CALL ks2d2s(p_cityAridIndex, p_cityCDDs, p_a,
           c            algAridIndex, algCDDs, p_a, d,prob, rr)
                 endif
                 
 165:            o_dists(x,y) = d
                 o_probs(x,y) = prob
                 
                 if (p_longs(x,y) >= -10 .and. p_longs(x,y) <= 35 .and.
           c        p_lats(x,y) >= 32 .and. p_lats(x,y) <= 70 ) then ! only analogues within within europe are considered
 170:            
                   if (d < bestD(worstBest) .and. prob >= 0.5 ) then ! consider everything with better than 50% probability
                     bestD(worstBest) = d
                     bestX(worstBest) = x
                     bestY(worstBest) = y
 175:              endif
                 
                   do i=1,best
                     if ( bestD(i) > bestD(worstBest) ) worstBest = i
                   enddo
 180:            
                 endif
         
              enddo
            enddo
 185:      
            ! sort the best candidates for neighborhood selection
            do i=1,best
              iBest = i
              do j=i,best
 190:             if (bestD(j) < bestD(iBest)) iBest = j
              enddo
              d = bestD(iBest)
              x = bestX(iBest)
              y = bestY(iBest)
 195:         bestD(iBest) = bestD(i)
              bestX(iBest) = bestX(i)
              bestY(iBest) = bestY(i)
              bestD(i) = d
              bestX(i) = x
 200:         bestY(i) = y
            enddo
           
            do i=1, best
              d = bestD(i)
 205:         if (d == 1000.) goto 2      !no best candidate in this slot
              x = bestX(i)
              y = bestY(i)
              do j=1,best
                if ( .not. (j == i) .and. .not. bestD(j) == 1000.  ) then
 210:             if ( (
           c        ( p_longs(x,y) - p_longs(bestX(j), bestY(j)) ) ** 2 +
           c        ( p_lats(x,y) - p_lats(bestX(j), bestY(j)) ) ** 2 )**0.5
           c        <= nD .and. d > bestD(j) ) then ! do not consider because better neighboring point exists
                        d = -1. 
 215:                   goto 2
                    endif
                endif
              enddo
                 
 220:         ! find cluster distance
              d = clusterDistance(p_x, p_y, p_longs, p_lats, o_dists, 
           c        x, y, 4, nD)
                 
              if ( (o_d == -1) .or. (d < o_d) ) then
 225:           o_d = d
                o_bestX = x
                o_bestY = y
              endif
                     
 230:    2    continue
      *         print *, "point ", i, ": ", x,y, "-",
      *     c    p_longs(x,y),"/",p_lats(x,y), ":", bestD(i), " calc: ", d
           
            enddo
 235:       END
      
      
            REAL FUNCTION clusterDistance (p_x, p_y,
           c      p_longs, p_lats, p_dists, p_curX, p_curY, 
 240:      c      p_weigh, p_neighborDist)
      ************************************************************************
      * Finds the K-S distance average value of a weighed point averaged with its
      * interpolated cardinal neighbors (at a given coordinate distance).
      *
 245: * Input:
      *  - p_x (Integer): grid size of data arrays
      *  - p_y (Integer): grid size of data arrays
      *  - p_longs (Real[p_x, p_y]): longitudes of data grid
      *  - p_lats (Real[p_x, p_y]): latitudes of data grid
 250: *  - p_dists (Real [p_x, p_y]): K-S distances of entire grid
      *  - p_curX (Integer): grid location of point whose weighted+neighbor distance is to be determined
      *  - p_curY (Integer): grid location of point whose weighted+neighbor distance is to be determined
      *  - p_weigh (Integer): how much to weigh the central point
      *  - p_neighborDist (Real): distance (in degree lat/lon) at which to interpolate the neighbors in the cardinal directions
 255: *
      * Return (Real): averaged distance
      ************************************************************************
            implicit none
            
 260:       !params
            INTEGER p_x, p_y ! grid size x/y
            REAL p_longs(p_x, p_y), p_lats(p_x, p_y) !longs and lats data, dah
            REAL p_dists(p_x, p_y) !K-S distances
            INTEGER p_curX, p_curY !point to evluate
 265:       INTEGER p_weigh ! weight of central point
            REAL p_neighborDist ! distance to neighbors
            
            ! process vars
            REAL interpolatePt
 270:       INTEGER cc 
            cc = 0
            
            ! evalute neighbor to the south
            clusterDistance = p_dists(p_curX, p_curY)*p_weigh
 275:       clusterDistance = clusterDistance +
           c      interpolatePt (p_longs(p_curX,p_curY)-p_neighborDist, 
           c        p_lats(p_curX,p_curY), 
           c        p_x, p_y, p_longs, p_lats, p_dists,
           c        1., cc )
 280: 
            ! neighbor to the north
            clusterDistance = clusterDistance +
           c      interpolatePt (p_longs(p_curX, p_curY)+p_neighborDist, 
           c        p_lats(p_curX, p_curY), 
 285:      c        p_x, p_y, p_longs, p_lats, p_dists,
           c        1., cc )
                     
            ! neighbor to the east
            clusterDistance = clusterDistance +
 290:      c      interpolatePt (p_longs(p_curX,p_curY), 
           c        p_lats(p_curX, p_curY)+p_neighborDist, 
           c        p_x, p_y, p_longs, p_lats, p_dists,
           c        1., cc )
      
 295:       ! neighbor to the left
            clusterDistance = clusterDistance +
           c      interpolatePt (p_longs(p_curX, p_curY), 
           c        p_lats(p_curX, p_curY)-p_neighborDist, 
           c        p_x, p_y, p_longs, p_lats, p_dists,
 300:      c        1., cc )
                    
            clusterDistance = clusterDistance / (4+p_weigh)
            END
      
 305: 
            SUBROUTINE sFindBestAlgs (p_x, p_y, p_a, 
           c p_curTemp, p_futTemp, p_curPrec, p_futPrec, 
           c p_curX, p_curY, p_algN, p_precMMax, p_precAMax,
           c o_bestX, o_bestY, o_bestTempDist, o_bestPrecMDist, 
 310:      c o_bestPrecADist)
      ************************************************************************
      * NOTE: this routine is not used in the program at the moment!
      * Selects a number (p_algN) of best analogues on the current grid for a 
      * location on the future grid by minimum temperature distance
 315: * within the allowed threshold of precipitation distance maxima
      * for THEORY: see Hallegatte et al. [2006]: Using climate analogues...
      *
      * Input:
      *  - p_x (Integer): grid size of data arrays
 320: *  - p_y (Integer): grid size of data arrays
      *  - p_a (Integer): number of years of data
      *  - p_curTemp (Real[p_x, p_y, p_a, 12]): temperature data that is to be compared with the city
      *  - p_futTemp (Real[p_x, p_y, p_a, 12]): temperature data from which to retrieve the city climate
      *  - p_curPrec (Real[p_x, p_y, p_a, 12]): precipitation data that is to be compared with the city
 325: *  - p_futPrec (Real[p_x, p_y, p_a, 12]): precipitation data from which retrieve the city climate
      *  - p_curX (Integer): location of the city on the future data grid
      *  - p_curY (Integer): location of the city on the future data grid
      *  - p_algN (Integer): number of best analogues to select
      *  - p_precMMax (Real): relative monthly precipitation distance threshold
 330: *  - p_precAMax (Real): relative annual precipitation distance threshold
      *
      * Output:
      *  - o_bestX (Integer[p_algN]): locations on grid of the best analogues found
      *  - o_bestY (Integer[p_algN]): locations on grid of the best analogues found
 335: *  - o_bestTempDist (Real[p_algN]): temperature distances of the best analogues found
      *  - o_bestPrecMDist (Real[p_algN]): relative monthly prec distances of the best analogues found
      *  - o_bestPrecADist (Real[p_algN]): relative annual prec distances of the best analogues found
      ************************************************************************
            IMPLICIT NONE
 340:             
            ! Constants      
            INTEGER p_x, p_y, p_a ! grid size x/y and number of years p_a
            INTEGER p_algN ! number of analogues to select
            REAL p_precMMax, p_precAMax ! max ratios allowed for precipitation distances
 345:       
            ! Data
            REAL p_curTemp(p_x,p_y,p_a,12), p_futTemp(p_x,p_y,p_a,12)
            REAL p_curPrec(p_x,p_y,p_a,12), p_futPrec(p_x,p_y,p_a,12)
            
 350:       ! Parameters
            REAL p_curX, p_curY ! city location on the grid
            
            ! Output
            INTEGER o_bestX(p_algN),o_bestY(p_algN)
 355:       REAL o_bestTempDist(p_algN)
            REAL o_bestPrecMDist(p_algN), o_bestPrecADist(p_algN)
            
            ! Vars
            INTEGER i,x,y !looping integer
 360:       INTEGER c !found analogues
            REAL tempD, precMD, precAD
            INTEGER worstI ! among the best scenarios for the output, this id's the worst      
            REAL worstTempD
            
 365:       print *, "Seach", p_algN, 
           c    " analogues with relative monthly precipitation max:",
           c    int(100.*p_precMMax),"%",
           c    "and relative annual precipitation max:",
           c    int(100.*p_precAMax),"% ..."
 370:       
            do i=1, p_algN
                  o_bestTempDist(i) = -1 ! comparison array initialization
            enddo
            
 375:       ! Search analogues ...
            worstI = 1
            worstTempD = o_bestTempDist(worstI)
            c=0
            do x=1, p_x
 380:         do y=1, p_y
              
                 CALL tempDist ( p_x, p_y, p_a, p_curTemp, x, y,
           c             p_futTemp, p_curX, p_curY, tempD) ! find temp dist
                 CALL precDist ( p_x, p_y, p_a, p_curPrec, x, y,
 385:      c             p_futPrec, p_curX, p_curY, precMD, precAD) ! find prec dists
      
                 if ( (precMD < p_precMMax) .and. (precAD < p_precAMax) ) then
                     ! analogue fits the criteria
                     c = c + 1
 390:                if ( (tempD < worstTempD) .or.
           c                  (worstTempD .eq. -1) ) then ! better analogue found
                        o_bestX(worstI) = x
                        o_bestY(worstI) = y
                        o_bestTempDist(worstI) = tempD
 395:                   o_bestPrecMDist(worstI) = precMD
                        o_bestPrecADist(worstI) = precAD
                        worstTempD = tempD
                     end if
                     
 400:                do i=p_algN, 1, -1 ! look which one is the worst analogue again
                        if ( (o_bestTempDist(i) .eq. -1) .or.
           c              ( (worstTempD < o_bestTempDist(i)) .and.
           c               (worstTempD .ne. -1) ) ) then
                            worstI = i
 405:                       worstTempD = o_bestTempDist(i)
                        endif
                     enddo
                 endif
           
 410:         enddo
            enddo
            
            if (c < p_algN) then
                  print *, 
 415:      c        "Only",c," analogues found that fit the criteria. ",
           c        "Please select a smaller number of analogues!" 
            else
                  print *, 
           c        p_algN, " out of", c, " found suitable analogues selected"
 420:       endif
                  
            END


Info Section
Warning: externals (function calls) may not be acurate calls: aridindex, ks2d2s, ks3d2s, precdist, tempdist yearsdegdays
back to top
f2html v0.3 (C) 1997,98 Beroud Jean-Marc. Fri Aug 11 17:54:58 CEST 2006