modelProcessing.f
SUBROUTINE processModel(
c overview, table, visualisation, wantedMethod,
c resultsP, modelName,
c x, y, futX, futY, tX, tY, a,
5: c curYearOffset, futYearOffset,
c citynum, cityName, cityLat, cityLon)
IMPLICIT none
************************************************************************
* Process one particular climate model and find the best analogues for
10: * all processed cities. I/O output is described thoroughly in INSTALL.TXT
*
* TODO: could netCDF files be potentially figured out automatically by 'cdo sinfo'?
*
* Input:
15: * - overview (Boolean): whether the overview maps should be created or not
* - table (Boolean): whether the p-value summary table should be created or not
* - visualisation (Boolean): whether K-S test visualization maps should be created or not
* - wantedMethod (Integer): which method should to be processed
* 0 - all methods used (prequisite for the p-value summary table being created)
20: * 1 - only 3D method w/ AI, HDD, CDD used (0 or 1 prerequisite for overview maps)
* 2 - only 2D method w/ AI, HDD used
* 3 - only 2D method w/ AI, CDD used
* - resultsP (String): folder in which to store the results
* - modelName (String): name of the model, currently supported:
25: * HC - HC climate simulation model
* CNRM - CNRM climate simulation model (ARPEGE)
* DKRZ+PRE - preindustrial CO2 conditions from DKRZ REMO model
* DKRZ+THC - preindustrial CO2 conditions w diminished THC from DKRZ REMO model
* - x, y (Integer, Integer): todays climate data grid
30: * - futX, futY (Integer, Integer): simulated climate data grid
* - tX, tY (Integer, Integer): grid for visualisation (to be interpolated into)
* - a (Integer): number of years data is from
* - curYearOffset (Integer): with which year the years of todays data coverage start (e.g. 1984)
* - futYearOffset (Integer): with which year the years of simulated data coverage start (e.g 2085)
35: * - citynum (Integer): number of cities to be processed
* - cityName (String[citynum]): Names of the cities to be processed
* - cityLat (Real[citynum]): Locations of the cities to be processed
* - cityLon (Real[citynum]): Locations of the cities to be processed
*
40: ************************************************************************
! methods available and used
INTEGER methodsN
PARAMETER(methodsN=3)
45: CHARACTER*20 methods(methodsN)
INTEGER method, method2, methodsStart, methodsEnd ! method working vars
! control vars
LOGICAL overview
50: LOGICAL table ! whether to calc and print p values in entire table
LOGICAL visualisation ! whether to create nc files of distance gradient
INTEGER wantedMethod ! which method is desired (0 is all, 1 is gonna be processed in any case if table is
REAL clusterDistance ! how far away is neighborhood
55: REAL precScaling ! how to convert to mm/day
REAL multi! multiplication factor in pvalues table
* *******************************************
* Sample cities (passed in from main routine)
60: * Data arrays, Calculation funcs and consts
* *******************************************
INTEGER citynum, i, j, k
CHARACTER*50 cityName(citynum)
REAL cityLat(citynum), cityLon(citynum)
65:
* Constants for the data
INTEGER x, y, a, m, curYearOffset, futYearOffset
INTEGER tX, tY ! vars for grids transposition
INTEGER futX, futY ! future data grid
70: CHARACTER resultsP*(*), resultsPath*100 ! results folder
CHARACTER modelName*(*), model*50 !model name
CHARACTER*(100)
c curCoordFile, futCoordFile,
c curTempFile, futTempFile, curPrecFile, futPrecFile,
75: c tempBgFile
PARAMETER (tempBgFile="data/europe_cur_temp.nc") !TODO: maybe put into configfile also
* Constants and functions for the calculation
REAL tempC2K
80: REAL ddBase
* Data arrays
REAL lats(x,y), longs(x,y), tLats(tY), tLongs(tX)
REAL futLats(futX, futY), futLongs(futX, futY)
85: REAL curTemp(x,y,a,12), futTemp(futX,futY,a,12) !, tTemp(tX,tY,a,m)
REAL curPrec(x,y,a,12), futPrec(futX,futY,a,12) !, tPrec(tX,tY,a,m)
INTEGER curYears(a), futYears(a)
* Vars and functions for the calculations
90: INTEGER strlen
INTEGER curX, curY ! processed city's current location on grid
REAL cityAridIndex(a) ! city's aridity index
REAL cityHDDs(a), cityCDDs(a) ! citiy's degree days
INTEGER algX, algY ! best analogue location on grid
95: REAL algWDist ! best analogue neighbor weighed distance
REAL dists(x,y), probs(x,y) ! results
REAL tDists(tX, tY) !interpolated transformed distances
INTEGER ic, iwmo, imod ! analogue parallel city search
100: CHARACTER station*20, country*10, minStation*20
REAL rlat, rlong, minLat, minLon, dist, minDist
CHARACTER citiesDb*100
PARAMETER(citiesDb='data/cities/cities.inv') !TODO: maybe put into configfile also
105: CHARACTER*100 overviewFName, fName ! working variables for composing file names
INTEGER b,c,d,e !looping var integers
REAL tempVar !working variable
LOGICAL fExist !working var for file existance checks
110: * Collected results from the calculastions
REAL ksdists(citynum, methodsN, methodsN)
REAL methodsDists(methodsN,x,y), methodsProbs(methodsN,x,y)
REAL algDists(methodsN), algProbs(methodsN)
INTEGER algXs(methodsN), algYs(methodsN)
115: CHARACTER*50 algCities(citynum, methodsN)
!*********************
!** Set config vars **
120: !*********************
model=modelName
resultsPath = resultsP
methods(1) = "_AI+HDD+CDD_"
methods(2) = "_AI+HDD_"
125: methods(3) = "_AI+CDD_"
ddBase=tempC2K(18.) ! degree days base from °C to Kelvin
if (model == "DKRZ+THC" .or. model=="DKRZ+PRE" ) then
precScaling = 1/30. !TODO should this be more exact? from config file? 1/30 for dkrz data where prec in mm/m
else
130: precScaling = 1. !HC and CNRM where prec is already in mm/d
endif
multi=9. ! multiplication factor for table summary graph size
if (wantedMethod == 0) then ! process all
135: methodsStart = 1
methodsEnd = methodsN
else !process only 1 method
methodsStart = wantedMethod
methodsEnd = wantedMethod
140: endif
if (wantedMethod > 0) table = .false. ! table can only be created if all methods are processed
if (wantedMethod > 1) overview = .false. ! make sure to create no overview if method 1 is not processed
145: !*********************
!*** Begin routine ***
!*********************
print *, ""
150: print *, "****************************"
print *, "Log for Analogue search"
print *, "Model: ", model
if (wantedMethod == 0) print *, "Method: all"
if (wantedMethod > 0) print *, "Method: ", methods(wantedMethod)
155: if (overview) print *, "Generate overview map: yes"
if (.not.overview) print *, "Generate overview map: no"
if (table) print *, "Generate summary table: yes"
if (.not.table) print *, "Generate summary table: no"
if (visualisation) print *,
160: c "Generate missing visualisations: yes"
if (.not.visualisation) print *,
c "Generate missing visualisations: no"
print *, "Process", citynum, " cities: ",
c (cityName(i)(1:strlen(cityName(i))+1), i=1, citynum)
165: print *, ""
!*********************
!*** Read the data ***
!*********************
170: ! TODO: make more dynamic, everything from config files!!
! maybe even directly from make written into configs?
if (model == "HC" .or.
c model == "DKRZ+PRE" .or. model == "DKRZ+THC" ) then
curCoordFile = "data/hc/t2m.HC.achgi.1960-1990.nc"
175: curTempFile = "data/hc/t2m.HC.achgi.1960-1990.nc"
curPrecFile = "data/hc/precip.HC.achgi.1960-1990.nc"
if (model == "HC") then !HC
futCoordFile = curCoordFile
futTempFile = "data/hc/t2m.HC.ackda.2070-2100.nc"
180: futPrecFile = "data/hc/precip.HC.ackda.2070-2100.nc"
elseif ( model == "DKRZ+PRE") then !DKRZ preindustrial
futCoordFile = "data/dkrz/temp_control.nc"
futTempFile = "data/dkrz/temp_control.nc"
futPrecFile = "data/dkrz/prec_control.nc"
185: else ! DKRZ THC
futCoordFile = "data/dkrz/temp_fut.nc"
futTempFile = "data/dkrz/temp_fut.nc"
futPrecFile = "data/dkrz/prec_fut.nc"
endif
190: elseif (model == "CNRM") then
curCoordFile = "data/cnrm/t2m.CNRM.DA9.nc"
curTempFile = "data/cnrm/t2m.CNRM.DA9.nc"
curPrecFile = "data/cnrm/precip.CNRM.DA9.nc"
futCoordFile = curCoordFile
195: futTempFile = "data/cnrm/t2m.CNRM.DE6.nc"
futPrecFile = "data/cnrm/precip.CNRM.DE6.nc"
endif
! current data
200: CALL getNetCdfCoordinates (
c curCoordFile(1:strlen(curCoordFile)), "Normal", x, y,
c "rlon", "rlat","lon","lat", longs, lats )
CALL getNetCdfVarVals (
205: c curTempFile(1:strlen(curTempFile)), "Normal", x, y, a,
c "rlon", "rlat", "time", "t2m", curTemp )
CALL getNetCdfVarVals (
c curPrecFile(1:strlen(curPrecFile)), "Normal", x, y, a,
210: c "rlon", "rlat", "time", "precip", curPrec )
if ( model == "HC" .or. model == "CNRM" ) then ! reading a normal netcdf file
CALL getNetCdfCoordinates (
215: c futCoordFile(1:strlen(futCoordFile)), "Normal", futX, futY,
c "rlon", "rlat","lon","lat", futLongs, futLats )
CALL getNetCdfVarVals (
c futTempFile(1:strlen(futTempFile)), "Normal", futX, futY, a,
220: c "rlon", "rlat", "time", "t2m", futTemp )
CALL getNetCdfVarVals (
c futPrecFile(1:strlen(futPrecFile)), "Normal", futX, futY, a,
c "rlon", "rlat", "time", "precip", futPrec )
225:
else ! reading standard self-defined netcdf file, i.e. DKRZ
CALL getNetCdfCoordinates (
c futCoordFile(1:strlen(futCoordFile)), "SDF", futX, futY,
230: c "lon", "lat","lon","lat", futLongs, futLats )
CALL getNetCdfVarVals (
c futTempFile(1:strlen(futTempFile)), "SDF", futX, futY, a,
c "lon", "lat", "time", "var167", futTemp )
235:
CALL getNetCdfVarVals (
c futPrecFile(1:strlen(futPrecFile)), "SDF", futX, futY, a,
c "lon", "lat", "time", "var260", futPrec )
240: endif
* years array creation
CALL createYearsArray (a, curYearOffset, curYears)
CALL createYearsArray (a, futYearOffset, futYears)
245:
!*********************
!*** Results files ***
!*********************
250: * if overview is to be created, open the overview grads files
if (overview) then
overviewFName = resultsPath(1:strlen(resultsPath)) //
c 'map_' // model(1:strlen(model))
CALL openSimpleGS (73, ! black and white map
255: c overviewFName(1:strlen(overviewFName)) // "_bw",
c tempBgFile, "var")
WRITE (73, *) "'c'"
WRITE (73, *) "'set black -100 100'"
WRITE (73, *) "'set grads off'"
260: WRITE (73, *) "'d var'"
CALL openSimpleGS (74, ! color map
c overviewFName(1:strlen(overviewFName)) // "_col",
c tempBgFile, "var")
265: WRITE (74, *) "'run grads/hcbar'"
endif
* if table is to be created, open the table files
if (overview) then
270: OPEN(21, ! tex table
c file=resultsPath(1:strlen(resultsPath)) // 'table_' //
c model(1:strlen(model)) // ".tex")
WRITE(21,*) "\\begin{table}"
WRITE(21,*) "\\begin{tabular}{llll}"
275: WRITE(21,*) "City & Analogue 3D & Analogue 2DH",
c " & Analogue 2DC \\\\ \\hline"
OPEN(22, ! plain text log table
c file=resultsPath(1:strlen(resultsPath)) // 'table_' //
280: c model(1:strlen(model)) // ".log")
endif
!*********************
285: !***** City loop *****
!*********************
do i=1, citynum
print *, ""
290: print *, "**************************************"
print *, "Find analogues for city",i," : ", cityName(i)
! coordinate translation
CALL translateCoordinates (futX, futY, futLats, futLongs,
295: c cityLat(i), cityLon(i), curX, curY)
! citys aridity index and HDD /CDDs
CALL aridIndex (futX, futY, a, futTemp, futPrec, precScaling,
c curX, curY, cityAridIndex)
300:
CALL yearsDegDays (futX, futY, a,
c futTemp, futYears, curX, curY,
c ddBase, cityHDDs, cityCDDs)
305: !*********************
!** Process Methods **
!*********************
do method=methodsStart, methodsEnd
310: print *, "Search with method ",
c methods(method)(1:strlen(methods(method))), " ..."
!find analogue
CALL avdFindBestAlg (x, y, longs, lats,
315: c curTemp, curPrec,
c a, curYears,
c cityAridIndex,cityHDDs,cityCDDs,
c ddBase, method,
c algX, algY, algWDist, dists, probs)
320:
algDists(method) = algWDist
algProbs(method) = probs(algX, algY)
if (algWDist > -1. ) then
325: algXs(method) = algX
algYs(method) = algY
if (table) then
! store all results for later usage for the summary table
do k=1,y
330: do j=1,x
methodsDists(method,j,k) = dists(j,k)
methodsProbs(method,j,k) = probs(j,k)
enddo
enddo
335: endif
endif
! log file entry
if ( algDists(method) > -1. ) then
340: print *,
c " Best analogue: ", algX, ";", algY,
c " (", lats(algX, algY),
c " /", longs(algX, algY), ")"
print *,
345: c " Probability: ", int(probs(algX, algY)*100),"%"
print *,
c " Test distance: ", dists(algX, algY)
print *,
c " Weighted distance: ", algWDist
350: else
print *, " No good analogue found"
endif
! search for city corresponding to the grid point
355: if ( algWDist > -1. ) then
minDist = -1.
open(unit=1,file=citiesDb(1:strlen(citiesDb)) )
100 continue
read(1,102,end=200) ic,iwmo,imod,station,country,rlat,rlong
360:
if ( rlong >= -12 .and. rlong <= 33.5
c .and. rlat >= 32.25 .and. rlat <= 70 ) then
! location approximately within europe, check for distance to best analogue
dist = (rlong-longs(algX, algY))**2+
365: c (rlat-lats(algX,algY))**2
if ( (dist .lt. minDist) .or. (minDist .eq. -1.) ) then
minDist = dist
minLat = rlat
minLon = rlong
370: minStation = station
endif
endif
102 format(i3,i5,i3,1x,a20,a10,1x,f6.2,1x,f7.2)
go to 100
375: 200 continue
close(1)
CALL wordCaseConvert (minStation) ! case convert
print *, " Closest station to analogue: ",
380: c minStation(1:strlen(minStation)),
c " at location", minLat, minLon
algCities(i,method) = minStation
else
algCities(i,method) = "no analogue"
385: endif
! visualize method result if visualisation is turned on the file does not yet exist
if (visualisation) then
390: print *, "Visualize result... "
fName = resultsPath(1:strlen(resultsPath))
c // model(1:strlen(model))
c // methods(method)(1:strlen(methods(method)))
395: c // cityName(i)
INQUIRE (file=fName(1:strlen(fName))//".nc",exist=fExist)
if (fExist) then
print *, "Visualization files already exist, skip."
400: goto 13 ! file exists skip the fisualisation
endif
! interpolation and plotting of distances in netCdf file
!FIXME not ncessary here, tLats, tLongs and tDists already known
405: print *, "Interpolate results to visualizable grid..."
CALL transform2RegGrid(x, y, lats, longs, dists,
c tX, tY, tLats, tLongs, tDists )
! create net cdf data file
410: print *, "Store results in netCdf data file..."
CALL netCDFCreate(fName, tX, tY,
c tLats, tLongs, tDists)
print *, "Create GrADS script file for visualizsation..."
415: ! grads script for net cdf with city location
CALL openSimpleGS(17, fName,
c fName(1:strlen(fName))//".nc", "var")
WRITE (17, *) "'c'"
420: WRITE (17, *) "'set rbcols 6 2 12 7 3 5 11 0 0'"
WRITE (17, *) "'set rbrange 0.3 1'"
WRITE (17, *) "'set grads off'"
WRITE (17, *) "'d var'"
425: ! original city
CALL drawGSPt (17, cityLat(i), cityLon(i),
c 0.1, 3, cityName(i))
! cross best analogue if there is one
430: if (algWDist > -1.) then
WRITE (17, *) "'set line 0'"
CALL drawGSPt (17,lats(algX,algY),longs(algX,algY),
c 0.2, 1, "")
endif
435:
CALL closeSimpleGS(17, fName)
13 continue
endif
440: enddo
! visualization of the 3D analogue on the overview map
if (overview) then
method = 1
445: if ( algDists(method) > -1.) then
print *, "Adding found 3D analogue to the overview map..."
! the general map (bw and color)
CALL drawGSPt (73, lats(algXs(method),algYs(method)),
c longs(algXs(method),algYs(method)),
450: c 0.1, 3, cityName(i))
CALL drawGSPt (74, lats(algXs(method),algYs(method)),
c longs(algXs(method),algYs(method)),
c 0.1, 3, cityName(i))
endif
455: endif
! create the comparison table of the different methods
! if this is desired
if (table) then
460: print *, "Comparing results of the processed methods..."
do method2=1,methodsN
! reverse translation, method2 distance results used to check how good other methods best distance results are on this one
do k=1,y
do j=1,x
465: probs(j,k) = methodsProbs(method2,j,k)
enddo
enddo
do method=1, methodsN
if (algDists(method) > -1. ) then ! analogue found for the method
470: ksdists(i,method, method2) =
c probs(algXs(method), algYs(method))
endif
enddo
enddo
475:
! generate tech p value table
print *, "Generating .tex table and log file entry..."
WRITE (21, *) cityname(i)(1:strlen(cityname(i)))
do method=1, methodsN
480: if (ksdists(i,method,method) == 0. ) then
WRITE (21, *) " & no good analogue"
else
WRITE (21, *) " & \\barres"
do method2=1, methodsN
485: tempVar = ksdists(i,method,method2)
if (tempVar == 0.) tempVar = 0.05
WRITE (21, *) "{", multi*tempVar, "}"
enddo
WRITE (21, *) "~~"
490: WRITE (21, 69)
c algCities(i,method),int(algProbs(method)*100)
69 format(a15,1x,i3,"\\%")
endif
enddo
495: WRITE (21, *) "\\\\"
! normal log file for easy overview and later reference if necessary
WRITE (22, *), cityname(i)(1:strlen(cityname(i)))
do method=1, methodsN
500: if (ksdists(i,method,method) == 0. ) then
WRITE (22, *) "no good analogue"
else
WRITE(22,55), algCities(i,method),
c (ksdists(i,method,method2), method2=1,methodsN)
505: endif
enddo
WRITE(22, *) "************************************"
55 format(a15,2x,f5.3,1x,f5.3,1x,f5.3)
endif
510:
enddo
!*********************
!**** End of loop ****
!****** Cleanup ******
515: !*********************
if (table) then ! close table files
close(22)
WRITE(21,*) "\\hline"
520: WRITE(21,*) "\\end{tabular}"
WRITE(21,*) "\\caption{Relative Performances}"
WRITE(21,*) "\\end{table}"
close(21)
endif
525:
if (overview) then ! close overview files
CALL closeSimpleGS (74,
c overviewFName(1:strlen(overviewFName)) // "_col")
CALL closeSimpleGS (73,
530: c overviewFName(1:strlen(overviewFName)) // "_bw")
endif
print *, "*******"
print *, "analogue search completed!"
535:
END
Info Section
Warning: externals (function calls) may not be acurate
calls: aridindex, avdfindbestalg, closesimplegs, createyearsarray, drawgspt
getnetcdfcoordinates, getnetcdfvarvals, netcdfcreate, opensimplegs, transform2reggrid
translatecoordinates, wordcaseconvert, yearsdegdays
back to top
f2html v0.3 (C) 1997,98 Beroud Jean-Marc. Fri Aug 11 17:54:58 CEST 2006