;+ ; ; ==================== ; calc_cloud_vlat.pro ; ==================== ; ; .. function:: calc_cloud_vlat(jday, swr, xlat) ; ; DESCRIPTION ; =========== ; ; Translated directly from Margie McCartys Fortran, ; added XLAT~=0 capability from WEARE et al 1980, ; and lats>20 capability from REED 1977. ; ; CALCULATE DAILY CLOUD FRACTION FROM DAILY SOLAR RADIATION, ; by comparing theoretical clear sky radiation to measured. ; ; INPUT: ; jday: value jday ; SWR: vector of SW rad in Watts/m^2 ; XLAT: corresponding vector of lat ; ; OUTPUT: ; CLD: cloud fraction - value between 0-1. ; Qclsky: clearsky Qsw = Q0 in Reeds notation ; ; SET COEFFICIENTS FOR THE Q0 CALCULATION ; ; EVOLUTIONS ; ========== ; ; $Id$ ; ; $URL$ ; ; - fplod 20120319 ; ; * minimal header ; ;- function calc_cloud_vlat, jday,SWR,XLAT ; ; pi=!pi xlat=double(xlat) NN=n_elements(xlat) ; CLD=SWR*0.+!values.f_nan Qclsky = SWR*0.+!values.f_nan ; A0 = (-32.65 + 674.76*cos(XLAT*pi/180)) A1 = (19.88 + 397.26*cos((XLAT+90.0)*pi/180)) A2 = (-1.32 + 16.10*sin((2*(XLAT-45.0))*pi/180)) B1 = (-6.75 + 224.38*sin(XLAT*pi/180)) B2 = (-1.04 + 29.76*cos((2*(XLAT-5.0))*pi/180)) ; ;A0 = (-32.65 + 674.76*cos(XLAT*pi/180))*(abs(XLAT) le 20.) ;A1 = (19.88 + 397.26*cos((XLAT+90.0)*pi/180))*(abs(XLAT) le 20.) ;A2 = (-1.32 + 16.10*sin((2*(XLAT-45.0))*pi/180))*(abs(XLAT) le 20.) ;B1 = (-6.75 + 224.38*sin(XLAT*pi/180))*(abs(XLAT) le 20.) ;B2 = (-1.04 + 29.76*cos((2*(XLAT-5.0))*pi/180))*(abs(XLAT) le 20.) ; ;A0 = A0+(-15.82 + 326.87*cos(XLAT*pi/180))*((abs(XLAT) gt 20.) and (abs(XLAT) le 40.)) ;A1 = A1+( 9.63 + 192.44*cos((XLAT+90.0)*pi/180))*((abs(XLAT) gt 20.) and (abs(XLAT) le 40.)) ;A2 = A2+( -0.64 + 7.80*sin((2*(XLAT-45.0))*pi/180))*((abs(XLAT) gt 20.) and (abs(XLAT) le 40.)) ;B1 = B1+( -3.27 + 108.70*sin(XLAT*pi/180))*((abs(XLAT) gt 20.) and (abs(XLAT) le 40.)) ;B2 = B2+( -0.50 + 14.42*cos((2*(XLAT-5.0))*pi/180))*((abs(XLAT) gt 20.) and (abs(XLAT) le 40.)) ; ;A0 = A0+(342.61 - 1.97*XLAT - 0.018*XLAT^2)*((abs(XLAT) gt 40.) and (abs(XLAT) le 60.)) ;A1 = A1+( 52.08 - 5.86*XLAT + 0.043*XLAT^2)*((abs(XLAT) gt 40.) and (abs(XLAT) le 60.)) ;A2 = A2+( 1.08 - 0.47*XLAT + 0.011*XLAT^2)*((abs(XLAT) gt 40.) and (abs(XLAT) le 60.)) ;B1 = B1+( -4.80 + 2.46*XLAT - 0.017*XLAT^2)*((abs(XLAT) gt 40.) and (abs(XLAT) le 60.)) ;B2 = B2+(-38.79 + 2.43*XLAT - 0.034*XLAT^2)*((abs(XLAT) gt 40.) and (abs(XLAT) le 60.)) ; DAY=jday ; day is day of that year; XDIV=365.15 ; ; CALCULATE PHI PHI = (360./XDIV * (DAY-21.0))*pi/180. ; ; CALCULATE Q0; MULTIPLY BY .4846 TO CONVERT FROM ; LY/DAY TO W/M**2 Q0 = A0+A1*cos(PHI)+B1*sin(PHI)+A2*cos(2*PHI)+B2*sin(2*PHI); Q0=Q0*.4846; ; SNALPHA=sin(XLAT*pi/180)*sin(23.45*sin((DAY-82)*pi/180)*pi/180) + $ cos(XLAT*pi/180)*cos(23.45*sin((DAY-82)*pi/180)*pi/180) ALPHA = asin(SNALPHA)*180/pi ; ; SEE MARINE CLIMATE ATLAS OF TROPICAL PACIFIC OCEAN, WEARE,STRUB, ; AND SAMUEL -- UC,DAVIS; THIS IS ACTUALLY REEDS 1977 ; FORMULA CLD = (1.0+0.0019*ALPHA-SWR/Q0)/.62 Qclsky = Q0 ; indx=where(CLD ge 1.0) if (indx(0) ne -1) then CLD(indx) = 1.+0.*indx ; ;return, [[CLD],[Qclsky]] return, cld ; end