REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_rho_eos.cpp
Go to the documentation of this file.
1#include <REMORA.H>
2
3using namespace amrex;
4
5/**
6 * @param[in ] bx box for calculation
7 * @param[in ] state state holds temp, salt
8 * @param[ out] rho density
9 * @param[ out] rhoA vertically-averaged density
10 * @param[ out] rhoS density perturbation
11 * @param[ out] bvf Brunt-Vaisala frequency
12 * @param[ out] alpha thermal expansion coefficient
13 * @param[ out] beta saline contraction coefficient
14 * @param[in ] Hz vertical cell height
15 * @param[in ] z_r z coordinates at rho points
16 * @param[in ] z_w z coordinates at w points
17 * @param[in ] h bathymetry
18 * @param[in ] mskr land-sea mask on rho points
19 * @param[in ] N number of vertical levels
20 */
21void
22REMORA::rho_eos (const Box& bx,
23 const Array4<Real const>& state,
24 const Array4<Real >& rho,
25 const Array4<Real >& rhoA,
26 const Array4<Real >& rhoS,
27 const Array4<Real >& bvf,
28 const Array4<Real >& alpha,
29 const Array4<Real >& beta,
30 const Array4<Real const>& Hz,
31 const Array4<Real const>& z_w,
32 const Array4<Real const>& z_r,
33 const Array4<Real const>& h,
34 const Array4<Real const>& mskr,
35 const int N)
36{
37 if (solverChoice.eos_type == EOSType::linear) {
38 lin_eos(bx, state, rho, rhoA, rhoS, bvf, Hz, z_w, z_r, h, mskr, N);
39 } else if (solverChoice.eos_type == EOSType::nonlinear) {
40 nonlin_eos(bx, state, rho, rhoA, rhoS, bvf, alpha, beta, Hz, z_w, z_r, h, mskr, N);
41 } else {
42 Abort("Unknown EOS type in rho_eos");
43 }
44}
45
46/**
47 * @param[in ] bx box for calculation
48 * @param[in ] state state holds temp, salt
49 * @param[ out] rho density
50 * @param[ out] rhoA vertically-averaged density
51 * @param[ out] rhoS density perturbation
52 * @param[ out] bvf Brunt-Vaisala frequency
53 * @param[in ] Hz vertical cell height
54 * @param[in ] z_r z coordinates at rho points
55 * @param[in ] z_w z coordinates at w points
56 * @param[in ] h bathymetry
57 * @param[in ] mskr land-sea mask on rho points
58 * @param[in ] N number of vertical levels
59 */
60void
61REMORA::lin_eos (const Box& bx,
62 const Array4<Real const>& state,
63 const Array4<Real >& rho,
64 const Array4<Real >& rhoA,
65 const Array4<Real >& rhoS,
66 const Array4<Real >& bvf,
67 const Array4<Real const>& Hz,
68 const Array4<Real const>& z_w,
69 const Array4<Real const>& z_r,
70 const Array4<Real const>& h,
71 const Array4<Real const>& mskr,
72 const int N)
73{
74 BL_PROFILE("REMORA::lin_eos()");
75//
76 AMREX_ASSERT(bx.smallEnd(2) == 0 && bx.bigEnd(2) == N);
77//
78//=======================================================================
79// Linear equation of state.
80//=======================================================================
81//
82//
83//-----------------------------------------------------------------------
84// Compute "in situ" density anomaly (kg/m3 - 1000) using the linear
85// equation of state.
86//-----------------------------------------------------------------------
87//
88 Real R0 = solverChoice.R0;
89 Real S0 = solverChoice.S0;
90 Real T0 = solverChoice.T0;
91 Real Tcoef = solverChoice.Tcoef;
92 Real Scoef = solverChoice.Scoef;
93
94 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
95 {
96 rho(i,j,k) = (R0 - R0*Tcoef*(state(i,j,k,Temp_comp)-T0)
97 + R0*Scoef*(state(i,j,k,Salt_comp)-S0)
98 - Real(1000.0)) * mskr(i,j,0);
99 });
100
101//
102//-----------------------------------------------------------------------
103// Compute vertical averaged density (rhoA) and density perturbation
104// used (rhoS) in barotropic pressure gradient.
105//-----------------------------------------------------------------------
106//
107 Real cff2 =one/solverChoice.rho0;
108
109 ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
110 {
111 Real cff0 = rho(i,j,N)*Hz(i,j,N);
112 rhoS(i,j,0) = Real(0.5)*cff0*Hz(i,j,N);
113 rhoA(i,j,0) = cff0;
114
115 for (int k = 1; k <= N; ++k) {
116 Real cff1=rho(i,j,N-k)*Hz(i,j,N-k);
117 rhoS(i,j,0) += Hz(i,j,N-k)*(rhoA(i,j,0)+Real(0.5)*cff1);
118 rhoA(i,j,0) += cff1;
119 }
120
121 Real cff11 =one/(z_w(i,j,N+1)+h(i,j,0,0));
122
123 rhoA(i,j,0) *= cff2*cff11;
124
125 rhoS(i,j,0) *= two*cff11*cff11*cff2;
126 });
127
128 // Compute Brunt-Vaisala frequency (1/s2)
129 Real gorho0 = g / solverChoice.rho0;
130 // Really want enclosed nodes or something similar
131 Box box_w = bx;
132 box_w.surroundingNodes(2);
133 box_w.grow(IntVect(0,0,-1));
134
135 ParallelFor(box_w, [=] AMREX_GPU_DEVICE (int i, int j, int k)
136 {
137 bvf(i,j,k) = -gorho0 * (rho(i,j,k)-rho(i,j,k-1)) / (z_r(i,j,k)-z_r(i,j,k-1));
138 });
139}
140
141/**
142 * @param[in ] bx box for calculation
143 * @param[in ] state state holds temp, salt
144 * @param[ out] rho density
145 * @param[ out] rhoA vertically-averaged density
146 * @param[ out] rhoS density perturbation
147 * @param[ out] bvf Brunt-Vaisala frequency
148 * @param[ out] alpha thermal expansion coefficient
149 * @param[ out] beta saline contraction coefficient
150 * @param[in ] Hz vertical cell height
151 * @param[in ] z_r z coordinates at rho points
152 * @param[in ] z_w z coordinates at w points
153 * @param[in ] h bathymetry
154 * @param[in ] mskr land-sea mask on rho points
155 * @param[in ] N number of vertical levels
156 */
157void
158REMORA::nonlin_eos (const Box& bx,
159 const Array4<Real const>& state,
160 const Array4<Real >& rho,
161 const Array4<Real >& rhoA,
162 const Array4<Real >& rhoS,
163 const Array4<Real >& bvf,
164 const Array4<Real >& alpha,
165 const Array4<Real >& beta,
166 const Array4<Real const>& Hz,
167 const Array4<Real const>& z_w,
168 const Array4<Real const>& z_r,
169 const Array4<Real const>& h,
170 const Array4<Real const>& mskr,
171 const int N)
172{
173 BL_PROFILE("REMORA::nonlin_eos()");
174//
175 AMREX_ASSERT(bx.smallEnd(2) == 0 && bx.bigEnd(2) == N);
176
177 FArrayBox fab_den1 (bx,1,amrex::The_Async_Arena()); auto den1 = fab_den1.array();
178 FArrayBox fab_den (bx,1,amrex::The_Async_Arena()); auto den = fab_den.array();
179 FArrayBox fab_bulk (bx,1,amrex::The_Async_Arena()); auto bulk = fab_bulk.array();
180 FArrayBox fab_bulk0(bx,1,amrex::The_Async_Arena()); auto bulk0 = fab_bulk0.array();
181 FArrayBox fab_bulk1(bx,1,amrex::The_Async_Arena()); auto bulk1 = fab_bulk1.array();
182 FArrayBox fab_bulk2(bx,1,amrex::The_Async_Arena()); auto bulk2 = fab_bulk2.array();
183 FArrayBox fab_DbulkDT(bx,1,amrex::The_Async_Arena()); auto DbulkDT = fab_DbulkDT.array();
184 FArrayBox fab_Dden1DT(bx,1,amrex::The_Async_Arena()); auto Dden1DT = fab_Dden1DT.array();
185 FArrayBox fab_DbulkDS(bx,1,amrex::The_Async_Arena()); auto DbulkDS = fab_DbulkDS.array();
186 FArrayBox fab_Dden1DS(bx,1,amrex::The_Async_Arena()); auto Dden1DS = fab_Dden1DS.array();
187
188 bool bulk_fluxes = solverChoice.bulk_fluxes;
189//
190//=======================================================================
191// Non-linear equation of state.
192//=======================================================================
193 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
194 {
195 Real Tt = std::max(Real(-2.5), state(i,j,k,Temp_comp));
196 Tt = std::min(Real(40.0), Tt);
197 Real Ts = std::max(Real(0.0), state(i,j,k,Salt_comp));
198 Ts = std::min(Real(100.0), Ts);
199 Real sqrtTs = std::sqrt(Ts);
200
201 Real Tp = z_r(i,j,k);
202 Real Tpr10 = Real(0.1)*Tp;
203
204 Real C0 = Q00+Tt*(Q01+Tt*(Q02+Tt*(Q03+Tt*(Q04+Tt*Q05))));
205 Real C1 = U00+Tt*(U01+Tt*(U02+Tt*(U03+Tt*U04)));
206 Real C2 = V00+Tt*(V01+Tt*V02);
207
208 if (bulk_fluxes) {
209 Real dCdT0=Q01+Tt*(two*Q02+Tt*(Real(3.0)*Q03+Tt*(Real(4.0)*Q04+
210 Tt*Real(5.0)*Q05)));
211 Real dCdT1=U01+Tt*(two*U02+Tt*(Real(3.0)*U03+Tt*Real(4.0)*U04));
212 Real dCdT2=V01+Tt*two*V02;
213 // Compute d(den1)/d(S) and d(den1)/d(T) derivatives used in the
214 // computation of thermal expansion and saline contraction
215 // coefficients.
216
217 Dden1DS(i,j,k)=C1+Real(1.5)*C2*sqrtTs+two*W00*Ts;
218 Dden1DT(i,j,k)=dCdT0+Ts*(dCdT1+sqrtTs*dCdT2);
219 }
220
221 den1(i,j,k) = C0 + Ts*(C1+sqrtTs*C2+Ts*W00);
222
223 //-----------------------------------------------------------------------
224 // Compute secant bulk modulus.
225 //-----------------------------------------------------------------------
226 Real C3=A00+Tt*(A01+Tt*(A02+Tt*(A03+Tt*A04)));
227 Real C4=B00+Tt*(B01+Tt*(B02+Tt*B03));
228 Real C5=D00+Tt*(D01+Tt*D02);
229 Real C6=E00+Tt*(E01+Tt*(E02+Tt*E03));
230 Real C7=F00+Tt*(F01+Tt*F02);
231 Real C8=G01+Tt*(G02+Tt*G03);
232 Real C9=H00+Tt*(H01+Tt*H02);
233
234 bulk0(i,j,k)=C3+Ts*(C4+sqrtTs*C5);
235 bulk1(i,j,k)=C6+Ts*(C7+sqrtTs*G00);
236 bulk2(i,j,k)=C8+Ts*C9;
237 bulk (i,j,k)=bulk0(i,j,k)-Tp*(bulk1(i,j,k)-Tp*bulk2(i,j,k));
238
239 Real cff = one / (bulk(i,j,k) + Tpr10);
240 den(i,j,k) = (den1(i,j,k)*bulk(i,j,k)*cff - Real(1000.0)) * mskr(i,j,0);
241 // This line may need to move once bluk fluxes are added
242 rho(i,j,k) = den(i,j,k);
243
244 if (bulk_fluxes) {
245 Real dCdT3=A01+Tt*(two*A02+Tt*(Real(3.0)*A03+Tt*Real(4.0)*A04));
246 Real dCdT4=B01+Tt*(two*B02+Tt*Real(3.0)*B03);
247 Real dCdT5=D01+Tt*two*D02;
248 Real dCdT6=E01+Tt*(two*E02+Tt*Real(3.0)*E03);
249 Real dCdT7=F01+Tt*two*F02;
250 Real dCdT8=G02+Tt*two*G03;
251 Real dCdT9=H01+Tt*two*H02;
252
253 DbulkDS(i,j,k)=C4+sqrtTs*Real(1.5)*C5-
254 Tp*(C7+sqrtTs*Real(1.5)*G00-Tp*C9);
255 DbulkDT(i,j,k)=dCdT3+Ts*(dCdT4+sqrtTs*dCdT5)-
256 Tp*(dCdT6+Ts*dCdT7-
257 Tp*(dCdT8+Ts*dCdT9));
258 }
259 });
260
261//
262//-----------------------------------------------------------------------
263// Compute vertical averaged density (rhoA) and density perturbation
264// used (rhoS) in barotropic pressure gradient.
265//-----------------------------------------------------------------------
266//
267 Real cff_rho =one/solverChoice.rho0;
268
269 ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
270 {
271 Real cff0 = den(i,j,N)*Hz(i,j,N);
272 rhoS(i,j,0) = Real(0.5)*cff0*Hz(i,j,N);
273 rhoA(i,j,0) = cff0;
274
275 for (int k = 1; k <= N; ++k) {
276 Real cff1=den(i,j,N-k)*Hz(i,j,N-k);
277 rhoS(i,j,0) += Hz(i,j,N-k)*(rhoA(i,j,0)+Real(0.5)*cff1);
278 rhoA(i,j,0) += cff1;
279 }
280 Real cff11 =one/(z_w(i,j,N+1)+h(i,j,0,0));
281 rhoA(i,j,0) *= cff_rho*cff11;
282 rhoS(i,j,0) *= two*cff11*cff11*cff_rho;
283 });
284
285 // Compute Brunt-Vaisala frequency (1/s2)
286 Box bxD = bx; bxD.makeSlab(2,0);
287
288 ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
289 {
290 bvf(i,j,0) = Real(0.0);
291 bvf(i,j,N+1) = Real(0.0);
292 for (int k=0; k<=N-1; k++) {
293 Real bulk_up = bulk0(i,j,k+1) - (z_w(i,j,k+1) * (bulk1(i,j,k+1) - bulk2(i,j,k+1)*z_w(i,j,k+1)));
294 Real bulk_dn = bulk0(i,j,k ) - (z_w(i,j,k+1) * (bulk1(i,j,k ) - bulk2(i,j,k )*z_w(i,j,k+1)));
295 Real cff1 = one / (bulk_up + Real(0.1) * z_w(i,j,k+1));
296 Real cff2 = one / (bulk_dn + Real(0.1) * z_w(i,j,k+1));
297 Real den_up = cff1 * (den1(i,j,k+1) * bulk_up);
298 Real den_dn = cff2 * (den1(i,j,k ) * bulk_dn);
299 bvf(i,j,k+1) = -g * (den_up - den_dn) / (Real(0.5) * (den_up+den_dn) * (z_r(i,j,k+1) - z_r(i,j,k)));
300 }
301 });
302
303 if (solverChoice.bulk_fluxes) {
304 ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
305 {
306 Real Tp = z_r(i,j,N);
307 Real Tpr10 = Real(0.1)*Tp;
308
309 // Compute thermal expansion and saline contraction coefficients
310 Real cff = bulk(i,j,N) + Tpr10;
311 Real cff1 = Tpr10 * den1(i,j,N);
312 Real cff2 = bulk(i,j,N) * cff;
313 Real wrk = (den(i,j,N) + Real(1000.0)) * cff * cff;
314 Real Tcof = -(DbulkDT(i,j,N) * cff1 + Dden1DT(i,j,N) * cff2);
315 Real Scof = (DbulkDS(i,j,N) * cff1 + Dden1DS(i,j,N) * cff2);
316 alpha(i,j,0) = Tcof / wrk;
317 beta(i,j,0) = Scof / wrk;
318 });
319 }
320}
A01
constexpr amrex::Real A01
Definition REMORA_Constants.H:75
U01
constexpr amrex::Real U01
Definition REMORA_Constants.H:107
E02
constexpr amrex::Real E02
Definition REMORA_Constants.H:88
W00
constexpr amrex::Real W00
Definition REMORA_Constants.H:114
E03
constexpr amrex::Real E03
Definition REMORA_Constants.H:89
H00
constexpr amrex::Real H00
Definition REMORA_Constants.H:97
Q04
constexpr amrex::Real Q04
Definition REMORA_Constants.H:104
two
constexpr amrex::Real two
Definition REMORA_Constants.H:8
V02
constexpr amrex::Real V02
Definition REMORA_Constants.H:113
B02
constexpr amrex::Real B02
Definition REMORA_Constants.H:81
U03
constexpr amrex::Real U03
Definition REMORA_Constants.H:109
H01
constexpr amrex::Real H01
Definition REMORA_Constants.H:98
one
constexpr amrex::Real one
Definition REMORA_Constants.H:7
G00
constexpr amrex::Real G00
Definition REMORA_Constants.H:93
D02
constexpr amrex::Real D02
Definition REMORA_Constants.H:85
A04
constexpr amrex::Real A04
Definition REMORA_Constants.H:78
B00
constexpr amrex::Real B00
Definition REMORA_Constants.H:79
A02
constexpr amrex::Real A02
Definition REMORA_Constants.H:76
B01
constexpr amrex::Real B01
Definition REMORA_Constants.H:80
F01
constexpr amrex::Real F01
Definition REMORA_Constants.H:91
g
constexpr amrex::Real g
Definition REMORA_Constants.H:38
Q00
constexpr amrex::Real Q00
Definition REMORA_Constants.H:100
G01
constexpr amrex::Real G01
Definition REMORA_Constants.H:94
G03
constexpr amrex::Real G03
Definition REMORA_Constants.H:96
Q02
constexpr amrex::Real Q02
Definition REMORA_Constants.H:102
H02
constexpr amrex::Real H02
Definition REMORA_Constants.H:99
Q01
constexpr amrex::Real Q01
Definition REMORA_Constants.H:101
Q03
constexpr amrex::Real Q03
Definition REMORA_Constants.H:103
U04
constexpr amrex::Real U04
Definition REMORA_Constants.H:110
F00
constexpr amrex::Real F00
Definition REMORA_Constants.H:90
B03
constexpr amrex::Real B03
Definition REMORA_Constants.H:82
E00
constexpr amrex::Real E00
Definition REMORA_Constants.H:86
Q05
constexpr amrex::Real Q05
Definition REMORA_Constants.H:105
G02
constexpr amrex::Real G02
Definition REMORA_Constants.H:95
V00
constexpr amrex::Real V00
Definition REMORA_Constants.H:111
U02
constexpr amrex::Real U02
Definition REMORA_Constants.H:108
F02
constexpr amrex::Real F02
Definition REMORA_Constants.H:92
E01
constexpr amrex::Real E01
Definition REMORA_Constants.H:87
U00
constexpr amrex::Real U00
Definition REMORA_Constants.H:106
A00
constexpr amrex::Real A00
Definition REMORA_Constants.H:74
V01
constexpr amrex::Real V01
Definition REMORA_Constants.H:112
A03
constexpr amrex::Real A03
Definition REMORA_Constants.H:77
D01
constexpr amrex::Real D01
Definition REMORA_Constants.H:84
D00
constexpr amrex::Real D00
Definition REMORA_Constants.H:83
EOSType::nonlinear
@ nonlinear
EOSType::linear
@ linear
Temp_comp
#define Temp_comp
Definition REMORA_IndexDefines.H:8
Salt_comp
#define Salt_comp
Definition REMORA_IndexDefines.H:9
T0
auto T0
Definition REMORA_InitAnalyticProb_Advection.H:11
S0
auto S0
Definition REMORA_InitAnalyticProb_Advection.H:12
REMORA::rho_eos
void rho_eos(const amrex::Box &bx, const amrex::Array4< amrex::Real const > &state, const amrex::Array4< amrex::Real > &rho, const amrex::Array4< amrex::Real > &rhoA, const amrex::Array4< amrex::Real > &rhoS, const amrex::Array4< amrex::Real > &bvf, const amrex::Array4< amrex::Real > &alpha, const amrex::Array4< amrex::Real > &beta, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &z_w, const amrex::Array4< amrex::Real const > &z_r, const amrex::Array4< amrex::Real const > &h, const amrex::Array4< amrex::Real const > &mskr, const int N)
Wrapper around equation of state calculation.
Definition REMORA_rho_eos.cpp:22
REMORA::lin_eos
void lin_eos(const amrex::Box &bx, const amrex::Array4< amrex::Real const > &state, const amrex::Array4< amrex::Real > &rho, const amrex::Array4< amrex::Real > &rhoA, const amrex::Array4< amrex::Real > &rhoS, const amrex::Array4< amrex::Real > &bvf, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &z_w, const amrex::Array4< amrex::Real const > &z_r, const amrex::Array4< amrex::Real const > &h, const amrex::Array4< amrex::Real const > &mskr, const int N)
Calculate density and related quantities from linear equation of state.
Definition REMORA_rho_eos.cpp:61
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1603
REMORA::nonlin_eos
void nonlin_eos(const amrex::Box &bx, const amrex::Array4< amrex::Real const > &state, const amrex::Array4< amrex::Real > &rho, const amrex::Array4< amrex::Real > &rhoA, const amrex::Array4< amrex::Real > &rhoS, const amrex::Array4< amrex::Real > &bvf, const amrex::Array4< amrex::Real > &alpha, const amrex::Array4< amrex::Real > &beta, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &z_w, const amrex::Array4< amrex::Real const > &z_r, const amrex::Array4< amrex::Real const > &h, const amrex::Array4< amrex::Real const > &mskr, const int N)
Calculate density and related quantities from nonlinear equation of state.
Definition REMORA_rho_eos.cpp:158
REMORA::bulk_fluxes
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.
Definition REMORA_bulk_flux.cpp:25
SolverChoice::Tcoef
amrex::Real Tcoef
Definition REMORA_DataStruct.H:859
SolverChoice::Scoef
amrex::Real Scoef
Definition REMORA_DataStruct.H:860