source: trunk/src/calc_cloud_vlat.pro

Last change on this file was 204, checked in by pinsard, 10 years ago

fix thanks to coding rules; typo

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