src/verif/tests/exactTestO.c

Go to the documentation of this file.

/*
   Copyright (C) 2011 Ed Bueler
  
   This file is part of PISM.
  
   PISM is free software; you can redistribute it and/or modify it under the
   terms of the GNU General Public License as published by the Free Software
   Foundation; either version 2 of the License, or (at your option) any later
   version.
  
   PISM is distributed in the hope that it will be useful, but WITHOUT ANY
   WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
   FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
   details.
  
   You should have received a copy of the GNU General Public License
   along with PISM; if not, write to the Free Software
   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
*/

#include "exactTestO.h"

#define beta_CC        7.9e-8       /* K Pa-1; Clausius-Clapeyron constant [\\ref Luethi2002] */
#define T_triple       273.15       /* K; triple point of pure water */
#define L              3.34e5       /* J kg-1; latent heat of fusion for water [\\ref AschwandenBlatter] */
#define grav           9.81         /* m/s^2; accel of gravity */

#define rho_ICE        910.0        /* kg/(m^3)  density of ice */
#define k_ICE          2.10         /* J/(m K s) = W/(m K)  thermal conductivity of ice */

#define k_BEDROCK      3.0          /* J/(m K s) = W/(m K)  thermal conductivity of bedrock */

#define H0             3000.0       /* m */
#define B0             1000.0       /* m */
#define Ts             223.15       /* K */
#define G              0.042        /* W/(m^2) = J m-2 s-1 */


int exactO(const double z, double *TT, double *Tm, double *qice, double *qbed, double *bmelt) {

  double P_base;

  P_base = rho_ICE * grav * H0;     /* Pa; hydrostatic approximation to pressure
                                       at base */

  *Tm = T_triple - beta_CC * P_base;/* K; pressure-melting point at base */

  *qice = - k_ICE * (Ts - *Tm) / H0;/* J m-1 K-1 s-1 K m-1 = J m-2 s-1 = W m-2;
                                       equilibrium heat flux everywhere in ice */

  *qbed = G;                        /* J m-2 s-1 = W m-2; equilibrium heat flux
                                       everywhere in bedrock */

  *bmelt = (G - *qice) / (L * rho_ICE);/* J m-2 s-1 / (J kg-1 kg m-3) = m s-1;
                                          ice-equivalent basal melt rate */

  if (z > H0) {
    *TT = Ts;                       /* K; in air above ice */
  } else if (z >= 0.0) {
    *TT = *Tm + (Ts - *Tm) * (z / H0); /* in ice */
  } else {
    *TT = *Tm - (G / k_BEDROCK) * z;   /* in bedrock */
  }
  return 0;
}