%	smalldyn.s
%
%	Purpose:	Small part of HIRLAM DYN specification for testing
%	Copyright:	R.A. van Engelen, Leiden University, 1997

include vecalg, arakawa.

% Dimensions

SomeDim := _G207699.
Time := "s".
Temperature := "K".
Distance := "m".
Pressure := "Pa".
DimLess := "1".
Velocity := Distance / Time.
Accelleration := Distance / Time ^ 2.
PressureTendency := Pressure / Time.
TemperatureTendency := Temperature / Time.

% Constants and parameters

r :: real ~ Distance := 6371 ~ "km".
	% radius of earth
R_d :: real ~ "J/kg/K" := 287.04 ~ "J/kg/K".
% gas constant for dry air
R_v :: real ~ "J/kg/K" := 461.51 ~ "J/kg/K".
% gas constant for water vapour
c_pd :: real ~ "J/kg/K" := 1004.64 ~ "J/kg/K".
% specific heat for dry air at constant pressure
c_pv :: real ~ "J/kg/K" := 1869.46 ~ "J/kg/K".
% specific heat for moist air at constant pressure

[delta, epsilon, kappa] isa parameter_var.

delta :: real := c_pv / c_pd.
epsilon :: real := R_d / R_v.
kappa :: real := R_d / c_pd.

% Grid size

natural := integer(1 .. infinity).

[nlon, nlat, nlev] isa parameter_var.

nlon :: natural.
nlat :: natural.
nlev :: integer(2 .. infinity).

% Vertical grid consists of at least 2 grid points.
% This constraint allows stronger optimization of vertical boundary conditions
% imposed on the eta_p field

% Computational domains

longitudinal := i = 1 .. nlon.
latitudinal := j = 1 .. nlat.
horizontal := longitudinal by latitudinal.
vertical := k = 1 .. nlev.
global := horizontal by vertical.
atmosphere := k = 2 .. nlev.

% Metric coefficients for the HIRLAM coordinate system

[h_x, h_y] isa parameter_var.
h_x :: real field (x, y).
h_y :: real field (x, y).

h_x is abstraction_of(grid, choices = (hxv, hxu, hxt)).
h_y is abstraction_of(grid, choices = (hyu, hyv, hyt)).

[hxu, hxv, hxt] isa parameter_var.
hxu :: real field (x(half), y(grid)) on horizontal.
hxv :: real field (x(grid), y(half)) on horizontal.
hxt :: real field (x(grid), y(grid)) on horizontal.
hxu has latex_name("h_{x,u}").
hxv has latex_name("h_{x,v}").
hxt has latex_name("h_{x,t}").

[hyu, hyv, hyt] isa parameter_var.
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.
hyu has latex_name("h_{y,u}").
hyv has latex_name("h_{y,v}").
hyt has latex_name("h_{y,t}").

H := [h_x, h_y].

% Use the chain-rule for symbolic differentiation and integration

include diff.

mode(diff, chain).

coordinates := [x, y].
coefficients := H.

eta := z ~ "1".
eta(N) := N.

% Grid-point distances for finite differences

rdlam isa parameter_var.
rdth isa parameter_var.
rdlam :: real(0 .. 1).
rdth :: real(0 .. 1).

delta x :: real := r / rdlam.
delta y :: real := r / rdth.
delta z :: real := 1.

% Override LaTeX descriptions of vertical finite differences

df_g(Expr, z) latex precedence(300); "\Delta^+_z", Expr.
df_h(Expr, z) latex precedence(300); "\Delta^-_z", Expr.

% Pressure

A isa parameter_var.
B isa parameter_var.
A :: real(0 .. 107000) ~ Pressure field z(grid) on vertical.
B :: real(0 .. 1) ~ DimLess field z(grid) on vertical.
p_s :: real(10000 .. 107000) ~ Pressure field ( x(grid)
                                              , y(grid)
                                              ) on horizontal.
p :: (real(0 .. 107000) ~ Pressure field ( x(grid)
                                         , y(grid)
                                         , z(grid)
                                         ) monotonic k(+) on global :=
      A + B * p_s
     ).


% Horizontal wind velocity vector

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

% Auxiliary horizontal wind velocity vector

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.

% Temperature (virtual)

T :: (real(0 .. infinity) ~ Temperature field ( x(grid)
                                              , y(grid)
                                              , z(half)
                                              ) on
      global
     ).

ltvir :: boolean.

ltvir := false.

T_v := T * (1 + (1 / epsilon - 1) * q) if ltvir \\ T otherwise.

% Surface pressure tendency

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

% Subroutine declaration

subroutine(smalldyn, [V_h, T], [p_s_t]).

% Fortran stuff

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).
T_t has fortran_name(PDTDT).