ocssw V2020
seawater.c
Go to the documentation of this file.
1 #include "l12_proto.h"
2 #include <gsl/gsl_fit.h>
3 
4 // ============================================================================
5 // Compute refractive index of seawater given temperature and salinity
6 // ============================================================================
7 
8 float seawater_nsw(float wave, float sst, float sss, float *dnswds) {
9  float n0 = 1.31405;
10  float n1 = 1.779e-4;
11  float n2 = -1.05e-6;
12  float n3 = 1.6e-8;
13  float n4 = -2.02e-6;
14  float n5 = 15.868;
15  float n6 = 0.01155;
16  float n7 = -0.00423;
17  float n8 = -4382.0;
18  float n9 = 1.1455e6;
19 
20  float n_air;
21  float n_sw;
22  float S;
23  float T;
24  float T2;
25  float wave2;
26  float wave3;
27 
28  S = sss;
29  T = sst;
30  T2 = T*T;
31  wave2 = wave*wave;
32  wave3 = wave2*wave;
33 
34  // refractive index of air is from Ciddor (1996, Applied Optics)
35  n_air = 1.0 + (5792105.0 / (238.0185 - 1 / (wave2 / 1e6)) + 167917.0 / (57.362 - 1 / (wave2 / 1e6))) / 1e8;
36 
37  // refractive index of seawater is from Quan and Fry (1994, Applied Optics)
38  n_sw = (n0 + (n1 + n2 * T + n3 * T2) * S + n4 * T2 + (n5 + n6 * S + n7 * T) / wave + n8 / wave2 + n9 / wave3) * n_air;
39 
40  if (dnswds != NULL) {
41  *dnswds = (n1 + n2 * T + n3 * T2 + n6 / wave) * n_air;
42  }
43 
44  return (n_sw);
45 }
46 
47 // ============================================================================
48 // Compute seawater isothermal compressibility
49 // ============================================================================
50 
51 float seawater_betat(float sst, float sss) {
52  float S = sss;
53  float T = sst;
54  float T2, T3, T4;
55  float kw;
56  float a0, b0, Ks;
57  float betat;
58 
59  T2 = T*T;
60  T3 = T2*T;
61  T4 = T2*T2;
62 
63  // pure water secant bulk Millero (1980, Deep-sea Research)
64 
65  kw = 19652.21 + 148.4206 * T - 2.327105 * T2 + 1.360477e-2 * T3 - 5.155288e-5 * T4;
66 
67  // seawater secant bulk
68 
69  a0 = 54.6746 - 0.603459 * T + 1.09987e-2 * T2 - 6.167e-5 * T3;
70  b0 = 7.944e-2 + 1.6483e-2 * T - 5.3009e-4 * T2;
71 
72  Ks = kw + a0 * S + b0 * pow(S, 1.5);
73 
74  // calculate seawater isothermal compressibility from the secant bulk
75 
76  betat = 1 / Ks * 1e-5; // unit is Pa
77 
78  return (betat);
79 }
80 
81 
82 // ============================================================================
83 // Compute density of seawater, unit is kg/m3, from UNESCO,38,1981
84 // ============================================================================
85 
86 float seawater_density(float sst, float sss) {
87  float S = sss;
88  float T = sst;
89  float T2, T3, T4, T5;
90  float density_w, density_sw;
91 
92  float a0 = 8.24493e-1;
93  float a1 = -4.0899e-3;
94  float a2 = 7.6438e-5;
95  float a3 = -8.2467e-7;
96  float a4 = 5.3875e-9;
97  float a5 = -5.72466e-3;
98  float a6 = 1.0227e-4;
99  float a7 = -1.6546e-6;
100  float a8 = 4.8314e-4;
101  float b0 = 999.842594;
102  float b1 = 6.793952e-2;
103  float b2 = -9.09529e-3;
104  float b3 = 1.001685e-4;
105  float b4 = -1.120083e-6;
106  float b5 = 6.536332e-9;
107 
108  T2 = T*T;
109  T3 = T2*T;
110  T4 = T2*T2;
111  T5 = T3*T2;
112 
113  // density for pure water
114 
115  density_w = b0 + b1 * T + b2 * T2 + b3 * T3 + b4 * T4 + b5*T5;
116 
117  // density for pure seawater
118 
119  density_sw = density_w + ((a0 + a1 * T + a2 * T2 + a3 * T3 + a4 * T4) * S + (a5 + a6 * T + a7 * T2) * pow(S, 1.5) + a8 * S * S);
120 
121  return (density_sw);
122 }
123 
124 
125 // ============================================================================
126 // Compute activity wrt salinity
127 // ============================================================================
128 
129 float seawater_dlnaswds(float sst, float sss) {
130  float S = sss;
131  float T = sst;
132  float T2, T3;
133  float dlnaswds;
134 
135  T2 = T*T;
136  T3 = T2*T;
137 
138  // water activity data of seawater is from Millero and Leung (1976,American
139  // Journal of Science,276,1035-1077). Table 19 was reproduced using
140  // Eqs.(14,22,23,88,107) then were fitted to polynominal equation.
141  // dlnawds is partial derivative of natural logarithm of water activity
142  // w.r.t.salinity
143 
144  dlnaswds = (-5.58651e-4 + 2.40452e-7 * T - 3.12165e-9 * T2 + 2.40808e-11 * T3) +
145  1.5 * (1.79613e-5 - 9.9422e-8 * T + 2.08919e-9 * T2 - 1.39872e-11 * T3) * sqrt(S) +
146  2 * (-2.31065e-6 - 1.37674e-9 * T - 1.93316e-11 * T2) * S;
147 
148  return (dlnaswds);
149 }
150 
151 
152 // ============================================================================
153 // Compute derivative of nsw wrt density via PMH model
154 // ============================================================================
155 
156 float seawater_dnswdrho(float n_sw) {
157  float n_sw2;
158  float dnswdrho;
159 
160  n_sw2 = n_sw*n_sw;
161 
162  dnswdrho = (n_sw2 - 1)*(1 + 2. / 3. * (n_sw2 + 2) * pow((n_sw / 3 - 1. / 3. / n_sw), 2));
163 
164  return (dnswdrho);
165 }
166 
167 
168 
169 // ============================================================================
170 // Compute backscattering coeff for seawater at temperature and salinity
171 //
172 // Ref: X. Zhang, L. Hu, and M-X. He, "Scattering by pure seawater: Effect of salinity,"
173 // Optics Express 17(7), 5698-5710 (2009).
174 // ============================================================================
175 
176 float seawater_bb(float wave, float sst, float sss, double delta) {
177  double S = sss;
178  double T = sst;
179 
180  // constants
181 
182  double Na = 6.0221417930e23; // Avagadro
183  double Kbz = 1.3806503e-23; // Boltzmann
184  double M0 = 18e-3; // molecular weight of water in kg/mol
185  // double delta = 0.039; // depolarization ratio from Farinato and Roswell 1975
186 
187  double Tk;
188  float n_sw, dnswds;
189  float isocomp;
190  float density_sw;
191  float dlnawds;
192  float dfri;
193  float beta_df;
194  float flu_con;
195  float beta_cf;
196  float beta90sw;
197  float bbsw;
198 
199  Tk = T + 273.15;
200 
201  // n_sw: absolute refractive index of seawater
202  // dnswds: partial derivative of seawater refractive index w.r.t. salinity
203 
204  n_sw = seawater_nsw(wave, T, S, &dnswds);
205 
206  // isothermal compressibility is from Lepple & Millero 1971
207  // the error ~ 0.004e-6 bar-1
208 
209  isocomp = seawater_betat(T, S);
210 
211  // density of water and seawater, unit is kg/m3, from UNESCO,38,1981
212 
213  density_sw = seawater_density(T, S);
214  // derivative of natural logarithm of water activity w.r.t.salinity
215 
216  dlnawds = seawater_dlnaswds(T, S);
217 
218  // density derivative of refractive index from PMH model
219 
220  dfri = seawater_dnswdrho(n_sw);
221 
222  // volume scattering at 90 degree due to the density fluctuation
223 
224  beta_df = PI * PI / 2 * (pow(wave * 1e-9, -4)) * Kbz * Tk * isocomp * dfri * dfri * (6 + 6 * delta) / (6 - 7 * delta);
225 
226  // volume scattering at 90 degree due to the concentration fluctuation
227 
228  flu_con = S * M0 * dnswds * dnswds / density_sw / (-dlnawds) / Na;
229  beta_cf = 2 * PI * PI * (pow(wave * 1e-9, -4)) * n_sw * n_sw * (flu_con)*(6 + 6 * delta) / (6 - 7 * delta);
230 
231  // total volume scattering at 90 degrees
232 
233  beta90sw = beta_df + beta_cf;
234 
235  // total backscattering coefficient
236 
237  bbsw = 8 * PI / 3 * beta90sw * (2 + delta) / (1 + delta) * 0.5;
238 
239  return (bbsw);
240 }
241 
254 void get_bbws(l2str *l2rec, l2prodstr *p, float prod[]) {
255  int i, ip;
256  int32_t nbands = l2rec->l1rec->l1file->nbands;
257  float *fwave = l2rec->l1rec->l1file->fwave;
258 
259  int nvisbands = 0;
260 
261  double x[nbands];
262  double y[nbands];
263 
264  for (i = 0; i < nbands; i++)
265  if (fwave[i] > 400 && fwave[i] < 700.) {
266  x[i] = log(fwave[i]);
267  nvisbands++;
268  }
269 
270  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
271  prod[ip] = p->badData;
272 
273  double c0, c1, cov00, cov01, cov11, sumsq;
274  for (i = 0; i < nvisbands; i++) {
275  int ibp = ip * nbands + i;
276  y[i] = log(l2rec->l1rec->sw_bb[ibp]);
277  }
278  // gsl_fit_mul(x,1,y,1,nbands, &c1, &cov11,&sumsq);
279  gsl_fit_linear(x, 1, y, 1, nvisbands,
280  &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
281 
282  prod[ip] = (float) c1;
283  }
284 }
#define NULL
Definition: decode_rs.h:63
float seawater_dlnaswds(float sst, float sss)
Definition: seawater.c:129
float seawater_dnswdrho(float n_sw)
Definition: seawater.c:156
unsigned long long sumsq(signed short *in, int cnt)
float seawater_bb(float wave, float sst, float sss, double delta)
Definition: seawater.c:176
void get_bbws(l2str *l2rec, l2prodstr *p, float prod[])
Definition: seawater.c:254
#define PI
Definition: l3_get_org.c:6
float seawater_betat(float sst, float sss)
Definition: seawater.c:51
float seawater_nsw(float wave, float sst, float sss, float *dnswds)
Definition: seawater.c:8
const double delta
instead the metadata field ProcessingEnvinronment is filled in from the output of a call to the POSIX compliant function uname from within the L1B code A small bug in L1B_Tables an incorrect comparison of RVS coefficients for TEBs to RVS coefficients for RSBs was being made This was replaced with a comparison between TEB coefficients This error never resulted in an incorrect RVS correction but did lead to recalculating the coefficients for each detector in a thermal band even if the coefficients were the same for all detectors To reduce to overall size of the reflective LUT HDF fill values were eliminated from all LUTs previously dimensioned where and where NUM_TIMES is the number of time dependent table pieces In Preprocess a small error where the trailing dropped scan counter was incremented when the leading dropped scan counter should have been was fixed This counter is internal only and is not yet used for any chiefly to casting of were added to make it LINUX compatible Output of code run on LINUX machines displays differences of at most scaled sector incalculable values of the Emissive calibration factor b1
Definition: HISTORY.txt:576
int32_t nbands
float seawater_density(float sst, float sss)
Definition: seawater.c:86
int i
Definition: decode_rs.h:71
float p[MODELMAX]
Definition: atrem_corl1.h:131