source: trunk/src/calc_cloud_vlat.pro @ 152

Last change on this file since 152 was 99, checked in by pinsard, 13 years ago

suppress successive semi-column lines

  • Property svn:keywords set to URL
File size: 2.8 KB
RevLine 
[88]1;+
2;-
[87]3 function calc_cloud_vlat, jday,SWR,XLAT
[97]4;
[87]5; Translated directly from Margie McCartys Fortran,
6; added XLAT~=0 capability from WEARE et al 1980,
7; and lats>20 capability from REED 1977.
8;
9;     CALCULATE DAILY CLOUD FRACTION FROM DAILY SOLAR RADIATION,
10;     by comparing theoretical clear sky radiation to measured.
11;
[97]12;  INPUT:
13;       jday: value jday
14;       SWR: vector of  SW rad in Watts/m^2
[87]15;       XLAT: corresponding vector of lat
16;
17;  OUTPUT:
18;       CLD:  cloud fraction - value between 0-1.
19;       Qclsky: clearsky Qsw = Q0 in Reeds notation
20;
21; SET COEFFICIENTS FOR THE Q0 CALCULATION
[97]22;
[87]23pi=!pi
[97]24xlat=double(xlat)
[94]25NN=n_elements(xlat)
[97]26;
[87]27CLD=SWR*0.+!values.f_nan
28Qclsky = SWR*0.+!values.f_nan
[97]29;
[87]30A0 = (-32.65 + 674.76*cos(XLAT*pi/180))
31A1 = (19.88 + 397.26*cos((XLAT+90.0)*pi/180))
32A2 = (-1.32 +  16.10*sin((2*(XLAT-45.0))*pi/180))
33B1 = (-6.75 + 224.38*sin(XLAT*pi/180))
34B2 = (-1.04 +  29.76*cos((2*(XLAT-5.0))*pi/180))
[97]35;
[87]36;A0 = (-32.65 + 674.76*cos(XLAT*pi/180))*(abs(XLAT) le 20.)
37;A1 = (19.88 + 397.26*cos((XLAT+90.0)*pi/180))*(abs(XLAT) le 20.)
38;A2 = (-1.32 +  16.10*sin((2*(XLAT-45.0))*pi/180))*(abs(XLAT) le 20.)
39;B1 = (-6.75 + 224.38*sin(XLAT*pi/180))*(abs(XLAT) le 20.)
40;B2 = (-1.04 +  29.76*cos((2*(XLAT-5.0))*pi/180))*(abs(XLAT) le 20.)
41;
42;A0 = A0+(-15.82 + 326.87*cos(XLAT*pi/180))*((abs(XLAT) gt 20.) and (abs(XLAT) le 40.))
43;A1 = A1+(  9.63 + 192.44*cos((XLAT+90.0)*pi/180))*((abs(XLAT) gt 20.) and (abs(XLAT) le 40.))
44;A2 = A2+( -0.64 +   7.80*sin((2*(XLAT-45.0))*pi/180))*((abs(XLAT) gt 20.) and (abs(XLAT) le 40.))
45;B1 = B1+( -3.27 + 108.70*sin(XLAT*pi/180))*((abs(XLAT) gt 20.) and (abs(XLAT) le 40.))
46;B2 = B2+( -0.50 +  14.42*cos((2*(XLAT-5.0))*pi/180))*((abs(XLAT) gt 20.) and (abs(XLAT) le 40.))
47;
48;A0 = A0+(342.61 - 1.97*XLAT - 0.018*XLAT^2)*((abs(XLAT) gt 40.) and (abs(XLAT) le 60.))
49;A1 = A1+( 52.08 - 5.86*XLAT + 0.043*XLAT^2)*((abs(XLAT) gt 40.) and (abs(XLAT) le 60.))
50;A2 = A2+(  1.08 - 0.47*XLAT + 0.011*XLAT^2)*((abs(XLAT) gt 40.) and (abs(XLAT) le 60.))
51;B1 = B1+( -4.80 + 2.46*XLAT - 0.017*XLAT^2)*((abs(XLAT) gt 40.) and (abs(XLAT) le 60.))
52;B2 = B2+(-38.79 + 2.43*XLAT - 0.034*XLAT^2)*((abs(XLAT) gt 40.) and (abs(XLAT) le 60.))
[97]53;
[87]54  DAY=jday ; day is day of that year;
55  XDIV=365.15
[97]56;
[87]57; CALCULATE PHI
58  PHI = (360./XDIV * (DAY-21.0))*pi/180.
[97]59;
[87]60; CALCULATE Q0; MULTIPY BY .4846 TO CONVERT FROM
61; LY/DAY TO W/M**2
62  Q0 = A0+A1*cos(PHI)+B1*sin(PHI)+A2*cos(2*PHI)+B2*sin(2*PHI);
63  Q0=Q0*.4846;
[97]64;
[87]65  SNALPHA=sin(XLAT*pi/180)*sin(23.45*sin((DAY-82)*pi/180)*pi/180) + $
66           cos(XLAT*pi/180)*cos(23.45*sin((DAY-82)*pi/180)*pi/180)
67  ALPHA = asin(SNALPHA)*180/pi
[97]68;
69; SEE MARINE CLIMATE ATLAS OF TROPICAL PACIFC OCEAN, WEARE,STRUB,
[87]70; AND SAMUEL -- UC,DAVIS;  THIS IS ACTUALLY REEDS 1977
71; FORMULA
72  CLD = (1.0+0.0019*ALPHA-SWR/Q0)/.62
73  Qclsky = Q0
[97]74;
75indx=where(CLD ge 1.0)
[87]76if (indx(0) ne -1) then CLD(indx) = 1.+0.*indx
[97]77;
78;return, [[CLD],[Qclsky]]
[87]79return, cld
[97]80;
[87]81end
Note: See TracBrowser for help on using the repository browser.