tneumann_bc.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
       ---
       tneumann_bc.mac (1223B)
       ---
            1 load("generic_equation.mac")$
            2 
            3 /* ethalpy equation for the positive and negative w[k] cases */
            4 eq[pm] := E[k-1]*l[pm] + E[k]*d[pm] + E[k+1]*u[pm] = b[pm];
            5 
            6 /* Additional equation implementing the Neumann B.C. */
            7 neumann[k] := (E[k+1] - E[k-1]) / (2 * dz) = G[k];
            8 
            9 /* Basal equation */
           10 E_base : solve(neumann[0], E[-1])[1];
           11 eq_base[pm] := facsum(subst([k=0, E_base], eq[pm]), E[0], E[1]);
           12 
           13 /* Matrix entries for the basal equation */
           14 d_base[pm] := facsum(coeff(lhs(eq_base[pm]), E[0]), w[0]);
           15 u_base[pm] := facsum(coeff(lhs(eq_base[pm]), E[1]), w[0]);
           16 
           17 b_base[pm] :=
           18 expandwrt(
           19   facsum(
           20     expand(
           21       rhs(eq_base[pm]) - (lhs(eq_base[pm]) - E[0]*d_base[pm] - E[1]*u_base[pm])),
           22     dz),
           23   rho);
           24 
           25 /* Surface equation */
           26 eq_surface[pm] := facsum(subst([k=k_s, E_surface], eq[pm]), E[k_s-1], E[k_s]);
           27 E_surface : solve(neumann[k_s], E[k_s+1])[1];
           28 
           29 /* Matrix entries for the surface equation */
           30 l_surface[pm] := facsum(coeff(lhs(eq_surface[pm]), E[k_s-1]), w[k_s]);
           31 d_surface[pm] := facsum(coeff(lhs(eq_surface[pm]), E[k_s]), w[k_s]);
           32 
           33 b_surface[pm] :=
           34 expandwrt(
           35   facsum(
           36     expand(
           37       rhs(eq_surface[pm]) - (lhs(eq_surface[pm]) - E[k_s-1]*l_surface[pm] - E[k_s]*d_surface[pm])),
           38     dz),
           39   rho);
           40 
           41 G_eq : G = phi / K;