Index: wflow-py/wflow/pcrut.py
===================================================================
diff -u -r3d8b5fe212103dd76367b768bb968449e5022d39 -rc7476868c506252d099ab57c9909e317b43aeeee
--- wflow-py/wflow/pcrut.py	(.../pcrut.py)	(revision 3d8b5fe212103dd76367b768bb968449e5022d39)
+++ wflow-py/wflow/pcrut.py	(.../pcrut.py)	(revision c7476868c506252d099ab57c9909e317b43aeeee)
@@ -253,107 +253,3 @@
 
     return rmat
 
-
-
-# def correctrad(Day,Hour,Lat,Lon,Slope,Aspect,Altitude):
-#     """ Determines radiation over a DEM assuming clear sky"""
-#
-#     Sc  = 1367.0          # Solar constant (Gates, 1980) [W/m2]
-#     Trans   = 0.6             # Transmissivity tau (Gates, 1980)
-#     pi = 3.1416
-#     AtmPcor = pow(((288.0-0.0065*Altitude)/288.0),5.256)
-#     Lat = Lat * pi/180
-#     ##########################################################################
-#     # Calculate Solar Angle and correct radiation ############################
-#     ##########################################################################
-#     # Solar geometry
-#     # ----------------------------
-#     # SolDec  :declination sun per day  between +23 & -23 [deg]
-#     # HourAng :hour angle [-] of sun during day
-#     # SolAlt  :solar altitude [deg], height of sun above horizon
-#     # SolDec  = -23.4*cos(360*(Day+10)/365);
-#     # Now added a new function that should work on all latitudes!
-#     theta    =(Day-1)*2 * pi/365  # day expressed in radians
-#
-#     SolDec =0.006918-0.399912 * cos(theta)+0.070257 * sin(theta) -   0.006758 * cos(2*theta)+0.000907 * sin(2*theta) -  0.002697 * cos(3*theta)+0.001480 * sin(3*theta)
-#
-#     HourAng = 180/pi * 15*(Hour-12.01)
-#     SolAlt  = scalar(asin(scalar(sin(Lat)*sin(SolDec)+cos(Lat)*cos(SolDec)*cos(HourAng))))
-#
-#     # Solar azimuth
-#     # ----------------------------
-#     # SolAzi  :angle solar beams to N-S axes earth [deg]
-#     SolAzi = scalar(acos((sin(SolDec)*cos(Lat)-cos(SolDec)* sin(Lat)*cos(HourAng))/cos(SolAlt)))
-#     SolAzi = ifthenelse(Hour <= 12, SolAzi, 360 - SolAzi)
-#
-#     # Surface azimuth
-#     # ----------------------------
-#     # cosIncident :cosine of angle of incident; angle solar beams to angle surface
-#     cosIncident = sin(SolAlt)*cos(Slope)+cos(SolAlt)*sin(Slope)*cos(SolAzi-Aspect)
-#
-#
-#     # Critical angle sun
-#     # ----------------------------
-#     # HoriAng  :tan maximum angle over DEM in direction sun, 0 if neg
-#     # CritSun  :tan of maximum angle in direction solar beams
-#     # Shade    :cell in sun 1, in shade 0
-#     # NOTE: for a changing DEM in time use following 3 statements and put a #
-#     #       for the 4th CritSun statement
-#     HoriAng   = cover(horizontan(Altitude,directional(SolAzi)),0)
-#     #HoriAng   = horizontan(Altitude,directional(SolAzi))
-#     HoriAng   = ifthenelse(HoriAng < 0, scalar(0), HoriAng)
-#     CritSun   = ifthenelse(SolAlt > 90, scalar(0), scalar(atan(HoriAng)))
-#     Shade   = SolAlt > CritSun
-#
-#     # Radiation outer atmosphere
-#     # ----------------------------
-#     OpCorr = Trans**((sqrt(1229+(614*sin(SolAlt))**2) -614*sin(SolAlt))*AtmPcor)    # correction for air masses [-]
-#     Sout   = Sc*(1+0.034*cos(360*Day/365)) # radiation outer atmosphere [W/m2]
-#     Snor   = Sout*OpCorr                   # rad on surface normal to the beam [W/m2]
-#
-#     # Radiation at DEM
-#     # ----------------------------
-#     # Sdir   :direct sunlight on a horizontal surface [W/m2] if no shade
-#     # Sdiff  :diffuse light [W/m2] for shade and no shade
-#     # Stot   :total incomming light Sdir+Sdiff [W/m2] at Hour
-#     # Radiation :avg of Stot(Hour) and Stot(Hour-HourStep)
-#     # NOTE: PradM only valid for HourStep & DayStep = 1
-#     SdirCor   = ifthenelse(Snor*cosIncident*scalar(Shade)<0,0.0,Snor*cosIncident*scalar(Shade))
-#     Sdir   = ifthenelse(Snor*cosIncident<0,0.0,Snor*cosIncident)
-#     Sdiff  = ifthenelse(Sout*(0.271-0.294*OpCorr)*sin(SolAlt)<0, 0.0, Sout*(0.271-0.294*OpCorr)*sin(SolAlt))
-#     AtmosDiffFrac = ifthenelse(Sdir > 0, Sdiff/Sdir, 1)
-#
-#     # Stot   = cover(Sdir+Sdiff,windowaverage(Sdir+Sdiff,3));     # Rad [W/m2]
-#     Stot   = Sdir + Sdiff                                     	    # Rad [W/m2]
-#     StotCor   = SdirCor + Sdiff                                   # Rad [W/m2]
-#
-#
-#     return scalar(SolAlt)
-#
-#
-# def GenRadMaps(SaveDir,SaveName,Lat,Lon,Slope,Aspect,Altitude,debug):
-#     """ Generates daily radiation maps for a whole year.
-#     It does so by running correctrad for a whole year with hourly
-#     steps and averaging this per day."""
-#
-#
-#     for Day in range(1,366):
-#     	avgrad = 0.0 * Altitude
-#     	for Hour in range(6,19):
-#     		print "day: " + str(Day) + " Hour: " + str(Hour)
-#     		crad = correctrad(Day,Hour,Lat,Lon,Slope,Aspect,Altitude)
-#     		if debug:
-#     			nr = "%0.3d" % Hour
-#     			report(crad,SaveDir + "/" + str(Day) + "/" + SaveName + "00000." + nr)
-#     		avgrad=avgrad + crad
-#
-#     	nr = "%0.3d" % Day
-#     	report(avgrad/13.0,SaveDir + "/" + SaveName + "00000." + nr)
-#
-
-
-
-    
-     
-
-    
\ No newline at end of file