real*8 function flux(R,decl,latitude,HA,albedo, & fracir,fracdust,surfaceSlope,azFac) C*********************************************************************** C flux: Insolation on horizontal or sloped surface C with a poor man's atmosphere C C R: distance from sun (AU) C decl: planetocentric solar declination (radians) C latitude: (radians) C HA: hour angle (radians from noon, clockwise) C fracir: fraction of absorption C fracdust: fraction of scattering C surfaceSlope: >0, (radians) C azFac: azimuth of gradient (radians east of north) C C written by Samar Khatiwala and Norbert Schorghofer C*********************************************************************** implicit none real*8 pi, So, d2r parameter (pi=3.1415926535897931, So=1365., d2r=pi/180.) real*8 R,decl,latitude,HA,albedo,fracIR,fracDust real*8 surfaceSlope,azFac real*8 c1,s1,Qo,solarAttenuation,Q real*8 sinbeta,sinbetaNoon,cosbeta,sintheta,azSun,buf ! solarAttenuation = (1.-albedo)*(1.-fracDust) c1=cos(latitude)*cos(decl) s1=sin(latitude)*sin(decl) C beta = 90 minus incidence angle for horizontal surface C beta = elevation of sun above (horizontal) horizon sinbeta = c1*cos(HA) + s1 sinbetaNoon = c1 + s1 sinbetaNoon = max(sinbetaNoon,0.d0) ! horizon cosbeta = sqrt(1-sinbeta**2) C ha -> az (option 1) ! azSun=asin(-cos(decl)*sin(ha)/cosbeta) C ha -> az (option 2) buf = (sin(decl)-sin(latitude)*sinbeta)/(cos(latitude)*cosbeta) if (buf>+1.) buf=1.d0 ! damn roundoff if (buf<-1.) buf=-1.d0 ! damn roundoff azSun = acos(buf) if (sin(HA)>=0) azSun=2*pi-azSun C ha -> az (option 3) without beta ! azSun=sin(latitude)*cos(decl)*cos(ha)-cos(latitude)*sin(decl) ! azSun=atan(sin(ha)*cos(decl)/azSun) ! print *,asin(sinbeta)/d2r,azSun/d2r C theta = 90 minus incidence angle for sloped surface sintheta = cos(surfaceSlope)*sinbeta + & sin(surfaceSlope)*cosbeta*cos(azSun-azFac) sintheta = max(sintheta,0.d0) ! horizon if (sinbeta<0.) sintheta=0. ! horizontal horizon at infinity Qo = So/(R**2) solarAttenuation = (1.-albedo)* & (1.- fracIR - fracDust)**(1./max(sinbeta,0.04)) ! new since May '05 C net flux: solar insolation + IR Q = sintheta*solarAttenuation + & cos(surfaceSlope/2.)**2*sinbetaNoon*fracIR ! new since May '05 if (sinbeta>0.d0) then Q = Q + 0.5*cos(surfaceSlope/2.)**2*fracDust ! new since May '05 endif flux=Qo*Q ! write(6,'(99(1x,f6.2))') decl/d2r,HA/d2r,flux, ! & asin(sintheta)/d2r,asin(sinbeta)/d2r,azSun/d2r,buf end