Ocean Color Science Software

 ocssw V2022
getglint.f
Go to the documentation of this file.
1 C Incident ray
2 C i = ( sin(PHI) , 0 , cos(PHI) )
3 C
4 C Reflected ray
5 C o = ( sin(MU)cos(NU) , sin(MU)sin(NU) , cos(MU) )
6 C
7 C Reflection angle
8 C i . o = cos(2 OMEGA) = sin(PHI)sin(MU)cos(NU)+cos(PHI)cos(MU)
9 C
10 C Normal to plane of reflection
11 C n = ( -sin(BETA)cos(ALPHA) , -sin(BETA)sin(ALPHA) , cos(BETA) )
12 C
13 C Reflection plane orientation
14 C n = i+o
15 C
16 C n = ( sin(PHI)+sin(MU)cos(NU) , sin(MU)sin(NU) , cos(PHI)+cos(MU) )/
17 C (2 cos(OMEGA))
18 C
19 C Tilt of plane
20 C cos(BETA) = (cos(PHI)+cos(MU))/(2 cos(OMEGA))
21 C
22 C Rotation of plane
23 C (i+n)^2 = 4 cos(OMEGA/2)
24 C
25 C cos(OMEGA) = cos(BETA)cos(PHI)-sin(BETA)sin(PHI)cos(ALPHA)
26 C
27 C See papers by Cox and Munk (1954a,b; 1956) for derivation of equations.
28
29
30  subroutine getglint(X1,X2,X3,X4,X5,X6)
31 C
32 C Calculate sun gliter coefficient
33 C
34 C X1 Angle MU (sensor zenith angle)
35 C X2 Angle PHI (solar zenith angle)
36 C X3 Angle NU (sensor-sun azimuth)
37 C X4 Wind speed (m/s)
38 C X5 Wind direction (radians)
39 C X6 Radiance (ignoring the atmosphere)
40 C
41 C #include "sensor_cmn.fin"
42  save
43  real pi
44  parameter(pi=3.1415926)
45  real X2, X3, X4, X5, acoss, RHO
46  real X1, Y1, RAD, X, Y2, OMEGA, BETA, ALPHA, SIGC, SIGU, Y3
47  real CHI, ALPHAP, SWIG, ETA, PROB, X6, EXPON
48  real Y4
49
50  rad(x) = x*pi/180.
51
52  y4 = max(x4,0.001)
53
54  y1 = rad(x1)
55  if (y1.eq.0.) y1 = 1.e-7
56  y2 = rad(x2)
57  if (y2.eq.0.) y2 = 1.e-7
58  y3 = rad(x3)
59  omega = acoss(cos(y1)*cos(y2)-sin(y1)*sin(y2)*cos(y3))/2.
60  if (omega.eq.0.) omega = 1.e-7
61  beta = acoss((cos(y1)+cos(y2))/(2.*cos(omega)))
62  if (beta.eq.0.) beta = 1.e-7
63  alpha = acoss((cos(beta)*cos(y2)-cos(omega))/(sin(beta)*sin(y2)))
64  if (sin(y3).lt.0.) alpha = -alpha
65
66 c Isotropic wind
67 c if (glintOn .ge. 2) then
68 c ! from Ebuchi & Kizu
69 c SIGC = SQRT(0.0048 + 0.00152*Y4)
70 c SIGU = SQRT(0.0053 + 0.000671*Y4)
71 c else
72 c ! from Cox & Munk
73  sigc = .04964*sqrt(y4)
74  sigu = .04964*sqrt(y4)
75 c endif
76  chi = x5
77  alphap = alpha+chi
78  swig = sin(alphap)*tan(beta)/sigc
79  eta = cos(alphap)*tan(beta)/sigu
80  expon = -(swig**2+eta**2)/2.
81  if (expon.lt.-30.) expon = -30. ! trap underflow
82  if (expon.gt.+30.) expon = +30. ! trap overflow
83  prob = exp(expon)/(2.*pi*sigu*sigc)
84  call reflec(omega,rho)
85
86 c Normal distribution
87  x6 = rho*prob/(4.*cos(y1)*cos(beta)**4)
88
89  return
90  end
91
92
93  subroutine getglint_iqu(X1,X2,X3,X4,X5,X6,X7,X8)
95 C Calculate sun gliter coefficient, including polarization
96 C
97 C X1(real) - Angle MU (sensor zenith angle)
98 C X2(real) - Angle PHI (solar zenith angle)
99 C X3(real) - Angle NU (sensor-sun azimuth)
100 C X4(real) - Wind speed (m/s)
101 C X5(real) - Wind direction
102 C X6(real) - Sun glitter coefficient
103 C X7(real) - Q/I for glitter
104 C X8(real) - U/I for glitter
105 C
106 C #include "sensor_cmn.fin"
107  save
108  real pi
109  parameter(pi=3.1415926)
110  real X2, X3, X4, X5, acoss, RHO_PLUS
111  real X1, Y1, RAD, X, Y2, OMEGA, BETA, ALPHA, SIGC, SIGU
112  real CHI, ALPHAP, SWIG, ETA, PROB, X6, EXPON
113  real X7, X8, RHO_MINUS, SR, CR, ROT_ANG, C2R, S2R, ASINN
114  real Y4
115
116  rad(x) = x*pi/180.
117
118  y4 = max(x4,0.001)
119
120  y1 = rad(x1)
121  if (y1.eq.0.) y1 = 1.e-7
122  y2 = rad(x2)
123  if (y2.eq.0.) y2 = 1.e-7
124  y3 = rad(x3)
125  omega = acoss(cos(y1)*cos(y2)-sin(y1)*sin(y2)*cos(y3))/2.
126  if (omega.eq.0.) omega = 1.e-7
127  beta = acoss((cos(y1)+cos(y2))/(2.*cos(omega)))
128  if (beta.eq.0.) beta = 1.e-7
129  alpha = acoss((cos(beta)*cos(y2)-cos(omega))/(sin(beta)*sin(y2)))
130  if (sin(y3).lt.0.) alpha = -alpha
131
132 c Isotropic wind
133 c if (glintOn .ge. 2) then
134 c ! from Ebuchi & Kizu
135 c SIGC = SQRT(0.0048 + 0.00152*Y4)
136 c SIGU = SQRT(0.0053 + 0.000671*Y4)
137 c else
138 c ! from Cox & Munk
139  sigc = .04964*sqrt(y4)
140  sigu = .04964*sqrt(y4)
141 c endif
142  chi = x5
143  alphap = alpha+chi
144  swig = sin(alphap)*tan(beta)/sigc
145  eta = cos(alphap)*tan(beta)/sigu
146  expon = -(swig**2+eta**2)/2.
147  if (expon.lt.-30.) expon = -30. ! trap underflow
148  if (expon.gt.+30.) expon = +30. ! trap overflow
149  prob = exp(expon)/(2.*pi*sigu*sigc)
150  call reflec_both(omega,rho_plus,rho_minus)
151
152 c Normal distribution
153  x6 = rho_plus*prob/(4.*cos(y1)*cos(beta)**4)
154
155 c Polarization components
156  if (omega .gt. .0001) then
157  cr = (cos(y2) - cos(2.*omega)*cos(y1))/(sin(2.*omega)*sin(y2))
158  sr = sin(y2)*sin(pi-y3) / sin(2.*omega)
159  rot_ang = sign(1., cr)*asinn(sr)
160  else
161  rot_ang = pi/2.
162  endif
163
164  c2r = cos(2.*rot_ang)
165  s2r = sin(2.*rot_ang)
166
167  x7 = c2r * rho_minus / rho_plus ! q_ov_i
168  x8 = -s2r * rho_minus / rho_plus ! u_ov_i
169
170  return
171  end
172
173
174  subroutine reflec (X1,X3)
175 C
176 C X1 Incident angle (radians)
177 C X3 Reflectance
178 C
179 C n1 sin(x1) = n2 sin(x2)
180 C
181 C tan(x1-x2)**2
182 C Refl(par ) = -------------
183 C tan(x1+x2)**2
184 C
185 C sin(x1-x2)**2
186 C Refl(perp) = -------------
187 C sin(x1+x2)**2
188 C
189 C Where:
190 C x1 Incident angle
191 C n1 Index refraction of Air
192 C x2 Refracted angle
193 C n2 Index refraction of Water
194 C
195  real X1, X2, X3, ref
196 * ! Index refraction of sea water
197  ref = 4./3.
198  if (x1.lt..00001) then
199  x3 = .0204078
200  else
201  x2 = asin(sin(x1)/ref)
202  x3 = (sin(x1-x2)/sin(x1+x2))**2+(tan(x1-x2)/tan(x1+x2))**2
203  x3 = x3/2.
204  end if
205  return
206  end
207
208
209  subroutine reflec_both (X1,X3,X4)
210 C
211 C X1 Incident angle (radians)
212 C X3 Reflectance sum
213 C X4 Reflectance difference
214 C
215 C n1 sin(x1) = n2 sin(x2)
216 C
217 C tan(x1-x2)**2
218 C Refl(par ) = -------------
219 C tan(x1+x2)**2
220 C
221 C sin(x1-x2)**2
222 C Refl(perp) = -------------
223 C sin(x1+x2)**2
224 C
225 C Where:
226 C x1 Incident angle
227 C n1 Index refraction of Air
228 C x2 Refracted angle
229 C n2 Index refraction of Water
230 C
231  real X1, X2, X3, X4, REF
232  real PERP, PAR
233 c ! Index refraction of sea water
234  ref = 4./3.
235  if (x1.lt..00001) then
236  x3 = .0204078
237  x4 = 0.
238  else
239  x2 = asin(sin(x1)/ref)
240  perp = (sin(x1-x2)/sin(x1+x2))**2
241  par = (tan(x1-x2)/tan(x1+x2))**2
242  x3 = perp+par
243  x3 = x3/2.
244  x4 = -perp+par
245  x4 = x4/2.
246  end if
247  return
248  end
249
250
251  function acoss (X1)
252  real acoss, x1
253 C
254 C This calculates ACOS(X) when X is near + or - 1
255 C
256  if (x1.ge.1.) then
257  acoss = 0.
258  else if (x1.le.-1.) then
259  acoss = 3.141592654
260  else
261  acoss = acos(x1)
262  end if
263  return
264  end
265
266
267  function asinn (X1)
268  real asinn, x1
269 C
270 C This calculates ASIN(X) when X is near + or - 1
271 C
272  if (x1.ge.1.) then
273  asinn = 3.141592654/2.0
274  else if (x1.le.-1.) then
275  asinn =-3.141592654/2.0
276  else
277  asinn = asin(x1)
278  end if
279  return
280  end
real function asinn(X1)
Definition: getglint.f:268
subroutine reflec_both(X1, X3, X4)
Definition: getglint.f:210
#define sign(x)
Definition: misc.h:95
real function acoss(X1)
Definition: getglint.f:252
subroutine reflec(X1, X3)
Definition: getglint.f:175
README for MOD_PR03(V6.1.0) 2. POINTS OF CONTACT it can be either SDP Toolkit or MODIS Packet for Terra input files The orbit validation configuration parameter(LUN 600281) must be either "TRUE" or "FALSE". It needs to be "FALSE" when running in Near Real Time mode
#define max(A, B)
Definition: main_biosmap.c:61
subroutine getglint(X1, X2, X3, X4, X5, X6)
Definition: getglint.f:31
subroutine getglint_iqu(X1, X2, X3, X4, X5, X6, X7, X8)
Definition: getglint.f:94