26 MultiFab* mf_Tair, MultiFab* mf_qair, MultiFab* mf_Pair,
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,
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));
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);
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);
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);
79 bool have_external_longwave = (mf_longwave_down !=
nullptr);
83 Real eps = Real(1e-20);
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));
90 ParallelFor(makeSlab(gbx1,2,0), [=] AMREX_GPU_DEVICE (
int i,
int j,
int ) {
92 Real PairM = Pair_arr(i,j,0);
93 Real TairC = Tair_arr(i,j,0);
94 Real TairK = TairC + Real(273.16);
95 Real Hair = qair_arr(i,j,0);
97 Real srflux = srflx_arr(i,j,0);
98 Real cloud = cloud_arr(i,j,0);
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);
109 Real LHeat = lhflx(i,j,0) * Hscale;
110 Real SHeat = shflx(i,j,0) * Hscale;
129 if (have_external_longwave && longwave_is_net) {
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);
135 LRad = Ldown - Lemit;
138 cff=(Real(0.7859)+Real(0.03477)*TairC)/(
one+Real(0.00412)*TairC);
139 Real e_sat=std::pow(Real(10.0),cff);
141 Real cff2=TairK*TairK*TairK;
142 Real cff1=cff2*TairK;
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));
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));
188 Real Qair = Real(0.62197)*(cff_saturation_air/(PairM-Real(0.378)*cff_saturation_air+eps));
193 Real cff_Q = cff_saturation_air*RH;
194 Q=Real(0.62197)*(cff_Q/(PairM-Real(0.378)*cff_Q+eps));
201 Real cff_saturation_water=(Real(1.0007)+Real(3.46e-6)*PairM)*Real(6.1121)*
206 Real cff_vp=cff_saturation_water*Real(0.98);
212 Real Qsea=Real(0.62197)*(cff_vp/(PairM-Real(0.378)*cff_vp+eps));
221 Real rhoAir=PairM*Real(100.0)/(
blk_Rgas*TairK*(
one+Real(0.61)*Q));
225 Real VisAir=Real(1.326e-5)*(
one+TairC*(Real(6.542e-3)+TairC*
226 (Real(8.301e-6)-Real(4.84e-9)*TairC)));
231 Real Hlv = (Real(2.501)-Real(0.00237)*cons(i,j,N,
Temp_comp))*Real(1.0e6);
237 Real delW=std::sqrt(wind_mag*wind_mag+Wgus*Wgus);
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));
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));
255 Real Ct=
vonKar/std::log(blk_ZT/ZoT10);
259 Real Ri = -
g*blk_ZW*((delT-delTc)+Real(0.61)*TairK*delQ)/
260 (TairK*delW*delW+eps);
263 Zetu=CC*Ri/(
one+Ri/Ribcu);
265 Zetu=CC*Ri/(
one+Real(3.0)*Ri/CC);
267 Real L10 = blk_ZW/Zetu;
270 Wstar=delW*
vonKar/(std::log(blk_ZW/Zo10)-
272 Real Tstar=-(delT-delTc)*
vonKar/(std::log(blk_ZT/ZoT10)-
274 Real Qstar=-(delQ-delQc)*
vonKar/(std::log(blk_ZQ/ZoT10)-
281 if (delW > Real(18.0)) {
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));
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;
294 Real ZoQ=std::min(Real(1.15e-4),Real(5.5e-5)/std::pow(Rr,Real(0.6)));
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);
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);
312 Real Bf=-
g/TairK*Wstar*(Tstar+Real(0.61)*TairK*Qstar);
318 delW=std::sqrt(wind_mag*wind_mag+Wgus*Wgus);
322 Real Wspeed=std::sqrt(wind_mag*wind_mag+Wgus*Wgus);
323 Cd=Wstar*Wstar/(Wspeed*Wspeed+eps);
326 Real Hs=-
blk_Cpa*rhoAir*Wstar*Tstar;
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))/
333 cff=Qair*Hlv/(
blk_Rgas*TairK*TairK);
334 Real wet_bulb=
one/(
one+Real(0.622)*(cff*Hlv*diffw)/
336 Real Hsr=rain(i,j,0)*wet_bulb*
blk_Cpw*
338 SHeat=(Hs+Hsr) * mskr(i,j,0);
342 Real Hl=-Hlv*rhoAir*Wstar*Qstar;
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);
351 Taur=Real(0.85)*rain(i,j,0)*wind_mag;
354 cff=rhoAir*Cd*Wspeed;
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);
391 lrflx(i,j,0) = LRad*Hscale2;
392 lhflx(i,j,0) = -LHeat*Hscale2;
393 shflx(i,j,0) = -SHeat*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) {
399 stflux(i,j,0,
Salt_comp) = mskr(i,j,0) * EminusP(i,j,0);
401 stflux(i,j,0,
Salt_comp) = mskr(i,j,0) * (evap(i,j,0)-rain(i,j,0)) /
rhow;
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);
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);