tflux_balance.mac - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tflux_balance.mac (1038B)
---
1 /* heat flux balance: LHS is ice, RHS is ocean */
2 heat_flux_balance : Qtb + Qti(grad_T) = Qtm$
3
4 /* salt flux balance */
5 salt_flux_balance : Qsb + Qsi = Qsm$
6
7 /* Solve the salt flux balance equation for the melt rate, then
8 substitute the melt rate and linearized equations for T_B and Theta_B
9 into the heat flux balance */
10 eq : heat_flux_balance, solve(salt_flux_balance, melt_rate), T_B = T_B(S_B, ice_thickness), Theta_B = Theta_B(S_B, ice_thickness)$
11 eq : eq - rhs(eq)$
12
13 /* expand the left hand side: */
14 eq_lhs : lhs(expand(eq*S_B/ocean_density))$
15
16 /* eq is a quadratic equation like this one: */
17 eq_quadratic : A*x^2 + B*x + C = 0$
18
19 /* with these coefficients */
20 A : factorsum(coeff(eq_lhs, S_B, 2))$
21 B : factorsum(coeff(eq_lhs, S_B, 1))$
22 C : factorsum(coeff(eq_lhs, S_B, 0))$
23
24 eq_melt_rate : factor(solve(salt_flux_balance, melt_rate)[1])$
25
26 tex('T_B(S_B,p) = T_B(S_B,p));
27 tex('Theta_B(S_B,p) = Theta_B(S_B,p));
28 tex(eq_quadratic);
29 tex('A = A);
30 tex('B = B);
31 tex('C = C);
32
33 tex(eq_melt_rate);
34
35 grind('A = A);
36 grind('B = B);
37