%	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.