% dyn.s
%
% Purpose: HIRLAM DYN specification
% Copyright: R.A. van Engelen, Leiden University, 1997
include vecalg, arakawa.
SomeDim := _G2510.
Time := "s".
Temperature := "K".
Distance := "m".
Pressure := "Pa".
DimLess := "1".
[Velocity] := [Distance] / [Time].
[Accelleration] := [Distance] / [Time] ^ 2.
[PressureTendency] := [Pressure] / [Time].
[TemperatureTendency] := [Temperature] / [Time].
r :: real ~ Distance := 6371 ~ "km".
% Gas constant for dry air
R_d :: real ~ "J/kg/K" := 287.04 ~ "J/kg/K".
% Gas constant for water vapour
R_v :: real ~ "J/kg/K" := 461.51 ~ "J/kg/K".
% Specific heat for dry air at constant pressure
c_pd :: real ~ "J/kg/K" := 1004.64 ~ "J/kg/K".
% Specific heat for moist air at constant pressure
c_pv :: real ~ "J/kg/K" := 1869.46 ~ "J/kg/K".
delta :: real := c_pv / c_pd.
epsilon :: real := R_d / R_v.
kappa :: real := R_d / c_pd.
natural := integer(1 .. infinity).
nlon :: natural.
% Grid size in latitudinal y-direction
nlat :: natural.
% Grid size in vertical z-direction
nlev :: integer(2 .. infinity).
longitudinal := i = 1 .. nlon.
latitudinal := j = 1 .. nlat.
horizontal := longitudinal by latitudinal.
vertical := k = 1 .. nlev.
global := horizontal by vertical.
atmosphere := k = 2 .. nlev.
h_x is abstraction_of(grid, choices = (hxv, hxu, hxt)).
h_x :: real field (x, y).
% Metric coefficient in the latitudinal y-direction
h_y is abstraction_of(grid, choices = (hyu, hyv, hyt)).
h_y :: real field (x, y).
H := [h_x, h_y].
hxv :: real field (x(grid), y(half)) on horizontal.
hxu :: real field (x(half), y(grid)) on horizontal.
hxt :: real field (x(grid), y(grid)) on horizontal.
hyu :: real field (x(half), y(grid)) on horizontal.
hyv :: real field (x(grid), y(half)) on horizontal.
hyt :: real field (x(grid), y(grid)) on horizontal.
include diff.
mode(diff, chain).
coordinates := [x, y].
coefficients := H.
eta := z ~ DimLess.
eta(k) := k.
rdlam :: real.
rdth :: real.
% Grid point distances in the x, y, and z-directions
delta x := r / rdlam.
delta y := r / rdth.
delta z := 1.
A :: real(0 .. 107000) ~ Pressure field z(grid) on vertical.
B :: real(0 .. 1) ~ DimLess field z(grid) on vertical.
% Surface pressure
p_s :: real(10000 .. 107000) ~ Pressure field ( x(grid)
, y(grid)
) on horizontal.
% Pressure
p :: real(0 .. 107000) ~ Pressure field ( x(grid)
, y(grid)
, z(grid)
) monotonic k(+) on global.
p = A + B * p_s.
alpha := ((1 - ( shift(p, z, - 1 / 2)
* d log(p) / d eta / (d p / d eta)
) if
on atmosphere
) \\
log(2) otherwise
).
beta := (d log(p) / d eta - alpha if on atmosphere \\
- log(2) otherwise
).
ln_p :: real ~ DimLess field ( x(grid)
, y(grid)
, z(half)
) on global.
ln_p = (d (p_bc * log(p_bc)) / d eta / (d p / d eta) where
p_bc = 0 if before atmosphere \\ p otherwise
).
u :: real ~ Velocity field ( x(half)
, y(grid)
, z(half)
) on global.
v :: real ~ Velocity field ( x(grid)
, y(half)
, z(half)
) on global.
V_h := [u, v].
u_aux :: (real ~ (Pressure * Velocity) field ( x(half)
, y(grid)
, z(half)
) on
global
).
v_aux :: (real ~ (Pressure * Velocity) field ( x(grid)
, y(half)
, z(half)
) on
global
).
V := [u_aux, v_aux].
V = V_h * d p / d eta.
f :: real ~ (1 / Time) field (x(half), y(half)) on horizontal.
q :: real field (x(grid), y(grid), z(half)) on global.
ltvir := true.
T :: (real(0 .. infinity) ~ Temperature field ( x(grid)
, y(grid)
, z(half)
) on
global
).
T_v := T * (1 + (1 / epsilon - 1) * q) if ltvir \\ T otherwise.
phi_s :: real ~ "m2/s2" field (x(grid), y(grid)) on horizontal.
phi :: real ~ "m2/s2" field ( x(grid)
, y(grid)
, z(half)
) on global.
phi = ( phi_s
+ alpha * R_d * T_v
+ sum((alpha + beta) * R_d * T_v, k = k + 1 .. nlev)
).
E :: real ~ "m2/s2" field ( x(grid)
, y(grid)
, z(half)
) on global.
E = ( 1 / 2
* ( 1 / h_y * ave(u ^ 2 * h_y, x)
+ 1 / h_x * ave(v ^ 2 * h_x, y)
)
).
Z :: real ~ "1/Pa/s" field ( x(half)
, y(half)
, z(half)
) on global.
Z = ( 1 / ave(ave(h_x * h_y * d p / d eta, x), y)
* ( (f * ave(ave(h_x * h_y, x), y) + d (h_y * v) / d x)
- d (h_x * u) / d y
)
).
p_s_t :: (real ~ PressureTendency field ( x(grid)
, y(grid)
, z(grid)
) on
horizontal
).
p_s_t = - integrate(nabla .* V, eta = eta(1) .. eta(nlev + 1)).
eta_p :: (real ~ PressureTendency field ( x(grid)
, y(grid)
, z(grid)
) on
global
).
eta_p = (1 - B) * p_s_t + integrate(nabla .* V, eta = eta .. eta(nlev + 1)).
% Definition of its LaTeX representation
eta_p latex "\dot\eta\frac{\partial p}{\partial\eta}".
advection(X) := adv_h(X) + adv_v(X).
(adv_h(X) :=
( ( ave(u_aux * h_y * d X / d x, x)
+ ave(v_aux * h_x * d X / d y, y)
)
/ (h_x * h_y * d p / d eta)
))
.
(adv_v(X) :=
(ave(eta_p_bc * d X / d eta, z) / (d p / d eta) where
eta_p_bc = eta_p if on atmosphere \\ 0 otherwise
))
.
omega := ( ( (kappa / (1 + (delta - 1) * q) if ltvir \\ kappa otherwise)
/ (d p / d eta)
)
* ( ( (d log(p) / d eta * (p_s_t + protect(integrate(nabla .* V, eta = eta(k + 1) .. eta(nlev + 1)), z)) if on atmosphere \\
0 otherwise
)
+ beta * (nabla .* V)
) * T_v
+ ( ( ave(u_aux * ave(T_v, x) * h_y * d ln_p / d x, x)
+ ave(v_aux * ave(T_v, y) * h_x * d ln_p / d y, y)
)
/ (h_x * h_y)
)
)
).
u_t :: real ~ Accelleration field ( x(half)
, y(grid)
, z(half)
) on global.
v_t :: real ~ Accelleration field ( x(grid)
, y(half)
, z(half)
) on global.
V_h_t := subscript(V_h, t).
% Note: this is a trick to get V_h_t := [u_t, v_t]
V_h_t = ( protect(Z) / H * ave(ave([[0, 1], [-1, 0]] &* V * H, x), y)
- (nabla * (phi + E) + R_d * T_v * nabla * ln_p)
- adv_v(V_h)
).
T_t :: (real ~ TemperatureTendency field ( x(grid)
, y(grid)
, z(half)
) on
global
).
T_t = - advection(T) + omega.
q_t :: (real ~ (DimLess / Time) field ( x(grid)
, y(grid)
, z(half)
) on
global
).
q_t = - advection(q).
subroutine(dyn, [V_h, T, q], [p_s_t, V_h_t, T_t, q_t]).
A has fortran_name(AHYB).
B has fortran_name(BHYB).
p_s has fortran_name(PS).
u has fortran_name(PUZ).
v has fortran_name(PVZ).
T has fortran_name(PTZ).
q has fortran_name(PQZ).
p_s_t has fortran_name(PDPSDT).
u_t has fortran_name(PDUDT).
v_t has fortran_name(PDVDT).
T_t has fortran_name(PDTDT).
q_t has fortran_name(PDQDT).
df_g(Expr, z) latex precedence(^); "\Delta^+_z", Expr.
df_h(Expr, z) latex precedence(^); "\Delta^-_z", Expr.