REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_bulk_flux.cpp
Go to the documentation of this file.
1#include <REMORA.H>
2
3using namespace amrex;
4
5/**
6 * @param[in ] lev level to operate on
7 * @param[in ] mf_cons scalar data: temperature, salinity, passsive scalar, etc
8 * @param[in ] mf_uwind u-direction wind dvelocity
9 * @param[in ] mf_vwind v-direction wind dvelocity
10 * @param[in ] mf_Tair air temperature [°C]
11 * @param[in ] mf_qair specific humidity [kg/kg]
12 * @param[in ] mf_Pair air pressure [mb]
13 * @param[in ] mf_srflx shortwave radiation flux [W/m²]
14 * @param[in ] mf_longwave_down external longwave radiation flux [W/m²]
15 * @param[inout] mf_evap evaporation rate
16 * @param[ out] mf_sustr u-direction surface momentum stress
17 * @param[ out] mf_svstr v-direction surface momentum stress
18 * @param[ out] mf_stflux surface scalar flux (temperature, salinity)
19 * @param[ out] mf_lrflx longwave radiation flux
20 * @param[inout] mf_lhflx latent heat flux
21 * @param[inout] mf_shflx sensible heat flux
22 * @param[in ] N number of vertical levels
23 */
24void
25REMORA::bulk_fluxes (int lev, MultiFab* mf_cons, MultiFab* mf_uwind, MultiFab* mf_vwind,
26 MultiFab* mf_Tair, MultiFab* mf_qair, MultiFab* mf_Pair,
27 MultiFab* mf_srflx,
28 MultiFab* mf_longwave_down,
29 MultiFab* mf_evap, MultiFab* mf_sustr, MultiFab* mf_svstr,
30 MultiFab* mf_stflux, MultiFab* mf_lrflx, MultiFab* mf_lhflx,
31 MultiFab* mf_shflx,
32 const int N)
33{
34 BL_PROFILE("REMORA::bulk_fluxes()");
35 const int IterMax = 3;
36 const BoxArray& ba = mf_cons->boxArray();
37 const DistributionMapping& dm = mf_cons->DistributionMap();
38 MultiFab mf_Taux(ba, dm, 1, IntVect(NGROW,NGROW,0));
39 MultiFab mf_Tauy(ba, dm, 1, IntVect(NGROW,NGROW,0));
40
41 // temps: Taux, Tauy,
42 for ( MFIter mfi(*mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
43 Array4<Real const> const& uwind = mf_uwind->const_array(mfi);
44 Array4<Real const> const& vwind = mf_vwind->const_array(mfi);
45 Array4<Real const> const& Tair_arr = mf_Tair->const_array(mfi);
46 Array4<Real const> const& qair_arr = mf_qair->const_array(mfi);
47 Array4<Real const> const& Pair_arr = mf_Pair->const_array(mfi);
48 Array4<Real const> const& srflx_arr = mf_srflx->const_array(mfi);
50 if (mf_longwave_down != nullptr) {
52 }
53 Array4<Real const> const& cons = mf_cons->const_array(mfi);
54 Array4<Real> const& sustr = mf_sustr->array(mfi);
55 Array4<Real> const& svstr = mf_svstr->array(mfi);
56 Array4<Real> const& stflux = mf_stflux->array(mfi);
57 Array4<Real> const& lrflx = mf_lrflx->array(mfi);
58 Array4<Real> const& lhflx = mf_lhflx->array(mfi);
59 Array4<Real> const& shflx = mf_shflx->array(mfi);
60 Array4<Real> const& evap = mf_evap->array(mfi);
61 Array4<Real> const& Taux = mf_Taux.array(mfi);
62 Array4<Real> const& Tauy = mf_Tauy.array(mfi);
63
64 Array4<const Real> const& mskr = vec_mskr[lev]->const_array(mfi);
65 Array4<const Real> const& msku = vec_msku[lev]->const_array(mfi);
66 Array4<const Real> const& mskv = vec_mskv[lev]->const_array(mfi);
67 Array4<const Real> const& rain = vec_rain[lev]->const_array(mfi);
68 Array4<const Real> const& EminusP = vec_EminusP[lev]->const_array(mfi);
69 Array4<const Real> const& cloud_arr = vec_cloud[lev]->const_array(mfi);
70
71 Real Hscale = solverChoice.rho0 * Cp;
72 Real Hscale2 = one / (solverChoice.rho0 * Cp);
73 Real blk_ZQ = solverChoice.blk_ZQ;
74 Real blk_ZT = solverChoice.blk_ZT;
75 Real blk_ZW = solverChoice.blk_ZW;
76 Real l_g = solverChoice.g;
77
79 bool longwave_is_net = solverChoice.longwave_is_net;
80 bool have_external_longwave = (mf_longwave_down != nullptr);
83
84 Real eps = Real(1e-20);
85
86 Box bx = mfi.tilebox();
87 Box ubx = mfi.grownnodaltilebox(0,IntVect(NGROW-1,NGROW-1,0));
88 Box vbx = mfi.grownnodaltilebox(1,IntVect(NGROW-1,NGROW-1,0));
89 Box gbx1 = bx; gbx1.grow(IntVect(NGROW,NGROW,0));
90
91 ParallelFor(makeSlab(gbx1,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) {
92 // Get spatially-varying atmospheric forcing from input arrays
93 Real PairM = Pair_arr(i,j,0); // Air pressure [mb]
94 Real TairC = Tair_arr(i,j,0); // Air temperature [°C]
95 Real TairK = TairC + Real(273.16); // Air temperature [K]
96 Real Hair = qair_arr(i,j,0); // Specific humidity [kg/kg] or RH [fraction]
97 Real RH = Hair;
98 Real srflux = srflx_arr(i,j,0); // Shortwave radiation flux [W/m²]
99 Real cloud = cloud_arr(i,j,0); // Cloud cover fraction [0-1]
100
101 // Input bulk parametrization fields
102 Real wind_mag = std::sqrt(uwind(i,j,0)*uwind(i,j,0) + vwind(i,j,0) * vwind(i,j,0)) + eps;
103 Real TseaK = cons(i,j,N,Temp_comp) + Real(273.16);
104
105 // Initialize
106 Real delTc = zero;
107 Real delQc = zero;
108 Real cff = zero;
109
110 Real LHeat = lhflx(i,j,0) * Hscale;
111 Real SHeat = shflx(i,j,0) * Hscale;
112 Real Taur = zero;
113 Taux(i,j,0) = zero;
114 Tauy(i,j,0) = zero;
115 Real LRad;
116
117 /*-----------------------------------------------------------------------
118 Compute outward or net longwave radiation (W/m2), LRad.
119 -----------------------------------------------------------------------
120 If external longwave radiation is net, use it directly. If it is
121 downward, compute net longwave radiation as Ldown - Lemit, where
122 Lemit is computed from the model SST and an emissivity. Otherwise,
123 use Berliand (1952) formula to calculate net longwave radiation.
124 The equation for saturation vapor pressure is from Gill (Atmosphere-
125 Ocean Dynamics, pp 606). Here the coefficient in the cloud term
126 is assumed constant, but it is a function of latitude varying from
127 1.0 at poles to 0.5 at the Equator).
128
129 */
130 if (have_external_longwave && longwave_is_net) {
131 // External forcing provides net longwave directly (W/m2).
134 Real Ldown = longwave_down_arr(i,j,0);
135 Real Lemit = emmiss * StefBo * std::pow(TseaK,4);
136 LRad = Ldown - Lemit;
137 } else {
138 // Original Berliand parameterization
139 cff=(Real(0.7859)+Real(0.03477)*TairC)/(one+Real(0.00412)*TairC);
140 Real e_sat=std::pow(Real(10.0),cff);
141 Real vap_p=e_sat*RH;
142 Real cff2=TairK*TairK*TairK;
143 Real cff1=cff2*TairK;
144
146 (cff1*(Real(0.39)-Real(0.05)*std::sqrt(vap_p))*
147 (one-Real(0.6823)*cloud*cloud)+
148 cff2*Real(4.0)*(TseaK-TairK));
149 }
150 /*
151 -----------------------------------------------------------------------
152 Compute specific humidities (kg/kg).
153
154 note that Qair is the saturation specific humidity at Tair
155 Q is the actual specific humidity
156 Qsea is the saturation specific humidity at Tsea
157
158 Saturation vapor pressure in mb is first computed and then
159 converted to specific humidity in kg/kg
160
161 The saturation vapor pressure is computed from Teten formula
162 using the approach of Buck (1981):
163
164 Esat(mb) = (Real(1.0007)+3.46E-Real(6)*PairM(mb))*Real(6.1121)*
165 EXP(Real(17.502)*TairC(C)/(Real(240.97)+TairC(C)))
166
167 The ambient vapor is found from the definition of the
168 Relative humidity:
169
170 RH = W/Ws*100 ~ E/Esat*100 E = RH/100*Esat if RH is in %
171 E = RH*Esat if RH fractional
172
173 The specific humidity is then found using the relationship:
174
175 Q = 0.622 E/(P + (0.622-1)e)
176
177 Q(kg/kg) = Real(0.62197)*(E(mb)/(PairM(mb)Real(-0.378)*E(mb)))
178
179 -----------------------------------------------------------------------
180 */
181
182 // Compute air saturation vapor pressure (mb), using Teten formula.
183
184 Real cff_saturation_air=(Real(1.0007)+Real(3.46e-6)*PairM)*Real(6.1121)*
185 std::exp(Real(17.502)*TairC/(Real(240.97)+TairC));
186
187 // Compute specific humidity at Saturation, Qair (kg/kg).
188
189 Real Qair = Real(0.62197)*(cff_saturation_air/(PairM-Real(0.378)*cff_saturation_air+eps));
190
191 // Compute specific humidity, Q (kg/kg).
192 Real Q;
193 if (RH < 2.0) {
194 Real cff_Q = cff_saturation_air*RH; //Vapor pressure (mb)
195 Q=Real(0.62197)*(cff_Q/(PairM-Real(0.378)*cff_Q+eps)); //Spec hum (kg/kg)
196 } else { // RH input was actually specific humidity in g/kg
197 Q=RH/Real(1000.0); //!Spec Hum (kg/kg)
198 }
199
200 // Compute water saturation vapor pressure (mb), using Teten formula.
201
202 Real cff_saturation_water=(Real(1.0007)+Real(3.46e-6)*PairM)*Real(6.1121)*
203 std::exp(Real(17.502)*cons(i,j,N,Temp_comp)/(Real(240.97)+cons(i,j,N,Temp_comp)));
204
205 // Compute water saturation vapor pressure (mb), using Teten formula.
206 // Vapor Pressure reduced for salinity (Kraus and Businger, 1994, pp42).
207 Real cff_vp=cff_saturation_water*Real(0.98);
208
209 // Compute Qsea (kg/kg) from vapor pressure.
210 // NOTE: ROMS does not have the small-value guard here, but does for
211 // Q and Qair
212
213 Real Qsea=Real(0.62197)*(cff_vp/(PairM-Real(0.378)*cff_vp+eps));
214 //
215 // -----------------------------------------------------------------------
216 // Compute Monin-Obukhov similarity parameters for wind (Wstar),
217 // heat (Tstar), and moisture (Qstar), Liu et al. (1979).
218 // -----------------------------------------------------------------------
219 //
220 // Moist air density (kg/m3).
221
222 Real rhoAir=PairM*Real(100.0)/(blk_Rgas*TairK*(one+Real(0.61)*Q));
223
224 // Kinematic viscosity of dry air (m2/s), Andreas (1989).
225
226 Real VisAir=Real(1.326e-5)*(one+TairC*(Real(6.542e-3)+TairC*
227 (Real(8.301e-6)-Real(4.84e-9)*TairC)));
228
229
230 // Compute latent heat of vaporization (J/kg) at sea surface, Hlv.
231
232 Real Hlv = (Real(2.501)-Real(0.00237)*cons(i,j,N,Temp_comp))*Real(1.0e6);
233
234 // Assume that wind is measured relative to sea surface and include
235 // gustiness.
236
237 Real Wgus=Real(0.5);
238 Real delW=std::sqrt(wind_mag*wind_mag+Wgus*Wgus);
239 Real delQ=Qsea-Q;
240 Real delT=cons(i,j,N,Temp_comp)-TairC;
241
242 // Neutral coefficients.
243 Real ZoW=Real(0.0001);
244 Real u10=delW*std::log(Real(10.0)/ZoW)/std::log(blk_ZW/ZoW);
245 Real Wstar=Real(0.035) * u10;
246 Real Zo10=Real(0.011)*Wstar*Wstar/l_g+Real(0.11)*VisAir/Wstar;
247 Real Cd10 =(vonKar/std::log(Real(10.0)/Zo10));
248 Cd10 = Cd10 * Cd10;
249 Real Ch10 =Real(0.00115);
250 Real Ct10 = Ch10/std::sqrt(Cd10);
251 Real ZoT10=Real(10.0)/std::exp(vonKar/Ct10);
252 Real Cd=(vonKar/std::log(blk_ZW/Zo10));
253 Cd = Cd * Cd;
254
255 // Compute Richardson number.
256 Real Ct=vonKar/std::log(blk_ZT/ZoT10); // T transfer coefficient
257 Real CC=vonKar*Ct/Cd;
258
259 Real Ribcu = -blk_ZW/(blk_Zabl*Real(0.004)*blk_beta*blk_beta*blk_beta);
260 Real Ri = -l_g*blk_ZW*((delT-delTc)+Real(0.61)*TairK*delQ)/
261 (TairK*delW*delW+eps);
262 Real Zetu;
263 if (Ri < zero) {
264 Zetu=CC*Ri/(one+Ri/Ribcu); // Unstable
265 } else {
266 Zetu=CC*Ri/(one+Real(3.0)*Ri/CC); // Stable
267 }
268 Real L10 = blk_ZW/Zetu;
269
270 // First guesses for Monon-Obukhov similarity scales.
271 Wstar=delW*vonKar/(std::log(blk_ZW/Zo10)-
272 bulk_psiu(blk_ZW/L10));
273 Real Tstar=-(delT-delTc)*vonKar/(std::log(blk_ZT/ZoT10)-
274 bulk_psit(blk_ZT/L10));
275 Real Qstar=-(delQ-delQc)*vonKar/(std::log(blk_ZQ/ZoT10)-
276 bulk_psit(blk_ZQ/L10));
277
278 // Modify Charnock for high wind speeds. The 0.125 factor below is for
279 // 1.0/(18.0-10.0).
280
281 Real charn;
282 if (delW > Real(18.0)) {
283 charn=Real(0.018);
284 } else if ((Real(10.0) < delW) and (delW <= Real(18.0))) {
285 charn=Real(0.011)+Real(0.125)*(Real(0.018)-Real(0.011))*(delW-Real(10.0));
286 } else {
287 charn=Real(0.011);
288 }
289
290 // Iterate until convergence. It usually converges within 3 iterations.
291 for (int it=0; it<IterMax; it++) {
292 ZoW=charn*Wstar*Wstar/l_g+Real(0.11)*VisAir/(Wstar+eps);
293 Real Rr=ZoW*Wstar/VisAir;
294 // Compute Monin-Obukhov stability parameter, Z/L.
295 Real ZoQ=std::min(Real(1.15e-4),Real(5.5e-5)/std::pow(Rr,Real(0.6)));
296 Real ZoT=ZoQ;
297 Real ZoL=vonKar*l_g*blk_ZW*(Tstar*(one+Real(0.61)*Q)+
298 Real(0.61)*TairK*Qstar)/
299 (TairK*Wstar*Wstar*(one+Real(0.61)*Q)+eps);
300 Real L=blk_ZW/(ZoL+eps);
301
302 // Evaluate stability functions at Z/L.
303 Real Wpsi=bulk_psiu(ZoL);
304 Real Tpsi=bulk_psit(blk_ZT/L);
305 Real Qpsi=bulk_psit(blk_ZQ/L);
306
307 // Compute wind scaling parameters, Wstar.
308 Wstar=std::max(eps,delW*vonKar/(std::log(blk_ZW/ZoW)-Wpsi));
309 Tstar=-(delT-delTc)*vonKar/(std::log(blk_ZT/ZoT)-Tpsi);
310 Qstar=-(delQ-delQc)*vonKar/(std::log(blk_ZQ/ZoQ)-Qpsi);
311
312 // Compute gustiness in wind speed.
313 Real Bf=-l_g/TairK*Wstar*(Tstar+Real(0.61)*TairK*Qstar);
314 if (Bf>zero) {
315 Wgus=blk_beta*std::pow(Bf*blk_Zabl,one/Real(3.0));
316 } else {
317 Wgus=Real(0.2);
318 }
319 delW=std::sqrt(wind_mag*wind_mag+Wgus*Wgus);
320 }
321
322 // Compute transfer coefficients for momentum (Cd).
323 Real Wspeed=std::sqrt(wind_mag*wind_mag+Wgus*Wgus);
325
326 // Compute turbulent sensible heat flux (W/m2), Hs.
328
329 // Compute sensible heat flux (W/m2) due to rainfall (kg/m2/s), Hsr.
330 Real diffw=Real(2.11e-5)*std::pow(TairK/Real(273.16),Real(1.94));
331 Real diffh=Real(0.02411)*(one+TairC*
332 (Real(3.309e-3)-Real(1.44e-6)*TairC))/
334 cff=Qair*Hlv/(blk_Rgas*TairK*TairK);
335 Real wet_bulb=one/(one+Real(0.622)*(cff*Hlv*diffw)/
336 (blk_Cpa*diffh));
337 Real Hsr=rain(i,j,0)*wet_bulb*blk_Cpw*
338 ((cons(i,j,N,Temp_comp)-TairC)+(Qsea-Q)*Hlv/blk_Cpa);
339 SHeat=(Hs+Hsr) * mskr(i,j,0);
340
341 // Compute turbulent latent heat flux (W/m2), Hl.
342
343 Real Hl=-Hlv*rhoAir*Wstar*Qstar;
344
345 // Compute Webb correction (Webb effect) to latent heat flux, Hlw.
346 Real upvel=Real(-1.61)*Wstar*Qstar-
347 (one+Real(1.61)*Q)*Wstar*Tstar/TairK;
348 Real Hlw=rhoAir*Hlv*upvel*Q;
349 LHeat=(Hl+Hlw) * mskr(i,j,0);
350
351 // Compute momentum flux (N/m2) due to rainfall (kg/m2/s).
352 Taur=Real(0.85)*rain(i,j,0)*wind_mag;
353
354 // Compute wind stress components (N/m2), Tau.
356 // amrex::Print() << "rhoAir: " << rhoAir << " Cd: " << Cd << " Wspeed: " << Wspeed << " cff: " << cff << "\n";
357 Real sign_u = (uwind(i,j,0) >= zero) ? 1 : -1;
358 Real sign_v = (vwind(i,j,0) >= zero) ? 1 : -1;
359 Taux(i,j,0)=(cff*uwind(i,j,0)+Taur*sign_u) * mskr(i,j,0);
360 Tauy(i,j,0)=(cff*vwind(i,j,0)+Taur*sign_v) * mskr(i,j,0);
361 // amrex::Print() << "Taux: " << Taux(i,j,0) << " Tauy: " << Tauy(i,j,0) << "\n";
362
363 //=======================================================================
364 // Compute surface net heat flux and surface wind stress.
365 //=======================================================================
366 //
367 // Compute kinematic, surface, net heat flux (degC m/s). Notice that
368 // the signs of latent and sensible fluxes are reversed because fluxes
369 // calculated from the bulk formulations above are positive out of the
370 // ocean.
371 //
372 // For EMINUSP option, EVAP = LHeat (W/m2) / Hlv (J/kg) = kg/m2/s
373 // PREC = rain = kg/m2/s
374 //
375 // To convert these rates to m/s divide by freshwater density, rhow.
376 //
377 // Note that when the air is undersaturated in water vapor (Q < Qsea)
378 // the model will evaporate and LHeat > 0:
379 //
380 // LHeat positive out of the ocean
381 // evap positive out of the ocean
382 //
383 // Note that if evaporating, the salt flux is positive
384 // and if raining, the salt flux is negative
385 //
386 // Note that stflux(:,:,isalt) is the E-P flux. The fresh water flux
387 // is positive out of the ocean and the salt flux is positive into the
388 // ocean. It is multiplied by surface salinity when computing state
389 // variable stflx(:,:,isalt) in "set_vbc.F".
390
391// Real one_over_rhow=one/rhow;
392 lrflx(i,j,0) = LRad*Hscale2;
393 lhflx(i,j,0) = -LHeat*Hscale2;
394 shflx(i,j,0) = -SHeat*Hscale2;
395 // Note: srflx from NetCDF is in W/m², convert to degC m/s by multiplying by Hscale2
396 stflux(i,j,0,Temp_comp)=(srflux*Hscale2 + lrflx(i,j,0) + lhflx(i,j,0) + shflx(i,j,0)) * mskr(i,j,0);
397 evap(i,j,0) = (LHeat / Hlv+eps) * mskr(i,j,0);
399 // Use prescribed E-P directly
400 stflux(i,j,0,Salt_comp) = mskr(i,j,0) * EminusP(i,j,0);
401 } else {
402 stflux(i,j,0,Salt_comp) = mskr(i,j,0) * (evap(i,j,0)-rain(i,j,0)) / rhow;
403 }
404 });
405
406 Real cff_rho = Real(0.5) / solverChoice.rho0;
407 ParallelFor(makeSlab(ubx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) {
408 sustr(i,j,0) = cff_rho*(Taux(i-1,j,0) + Taux(i,j,0)) * msku(i,j,0);
409 });
410 ParallelFor(makeSlab(vbx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) {
411 svstr(i,j,0) = cff_rho*(Tauy(i,j-1,0) + Tauy(i,j,0)) * mskv(i,j,0);
412 });
413 }
414
415}
constexpr amrex::Real one
constexpr amrex::Real zero
constexpr amrex::Real blk_Zabl
constexpr amrex::Real vonKar
constexpr amrex::Real rhow
constexpr amrex::Real blk_Rgas
constexpr amrex::Real blk_beta
constexpr amrex::Real emmiss
constexpr amrex::Real Cp
constexpr amrex::Real StefBo
constexpr amrex::Real blk_Cpa
constexpr amrex::Real blk_Cpw
#define NGROW
#define Temp_comp
#define Salt_comp
mf_h setVal(geomdata.ProbHi(2))
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_EminusP
evaporation minus precipitation [kg/m^2/s], defined at rho-points
Definition REMORA.H:521
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskr
land/sea mask at cell centers (2D)
Definition REMORA.H:584
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_rain
precipitation rate [kg/m^2/s]
Definition REMORA.H:515
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_msku
land/sea mask at x-faces (2D)
Definition REMORA.H:586
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskv
land/sea mask at y-faces (2D)
Definition REMORA.H:588
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real bulk_psiu(amrex::Real ZoL)
Evaluate stability function psi for wind speed.
Definition REMORA.H:2192
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1949
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_cloud
cloud cover fraction [0-1], defined at rho-points
Definition REMORA.H:519
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real bulk_psit(amrex::Real ZoL)
Evaluate stability function psi for moisture and heat.
Definition REMORA.H:2224
void bulk_fluxes(int lev, amrex::MultiFab *mf_cons, amrex::MultiFab *mf_uwind, amrex::MultiFab *mf_vwind, amrex::MultiFab *mf_Tair, amrex::MultiFab *mf_qair, amrex::MultiFab *mf_Pair, amrex::MultiFab *mf_srflx, amrex::MultiFab *mf_longwave_down, amrex::MultiFab *mf_evap, amrex::MultiFab *mf_sustr, amrex::MultiFab *mf_svstr, amrex::MultiFab *mf_stflux, amrex::MultiFab *mf_lrflx, amrex::MultiFab *mf_lhflx, amrex::MultiFab *mf_shflx, const int N)
Calculate bulk temperature, salinity, wind fluxes.
@ EminusP
evaporation minus precipitation [m/s]
std::array< BulkForcingType, BulkFlux::NumTypes > bulk_flux_type