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