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