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