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