source: issm/trunk-jpl/src/py3/boundaryconditions/PattynSMB.py@ 23670

import os
import numpy as  np
def PattynSMB(md,Tf):
        """
    PATTYNSMB- Compute SMB over Antarctica (from Pattyn 2006, pg. 18, "GRANTISM: An ExcelTM model for Greenland 
        and Antarctic ice-sheet response to climate changes")

    Usage:
      md=PattynSMB(md,Tf)

      where Tf is a background forcing temperature ("an anomalous temperature relative to the present conditions)


    See also: SETICESHELFBC, SETMARINEICESHEETBC
        """
 
        # Tma    : Mean annual surface temperature in [deg C]
        # Tms    : Mean summer temperature in [deg C]
        # h      : Surface/bedrock elevation (I assume in meters but paper does not specify)
        # phi    : Latitude in degrees SOUTH
        # lambda : Longitude in degrees WEST
        # Tf     : Background forcing temperature ("an anomalous temperature relative to the present conditions)
        # ACCdot : Accumulation rate in units of [m/a] ice equivalent
        # ABLdot : Surface ablation rate in [m/a] ice equivalent

        #Double check lat and long exist:
        if np.any(np.isnan(md.mesh.lat)): 
                raise IOError('PattynSMB error message: md.mesh.lat field required')

        # Calculate mean annual surface temperature, Eqn (11)
        # Here, -0.012 is the atmospheric Lapse rate from sea level in deg/m.
        # It is multiplied by surface elevation from sea level
        Tma = -15.15 - 0.012*md.geometry.surface
        

        # Calculate summer temperature, Eqn (12)
        # No melting at PIG in mean conditions - need about 6 degress Tf to start having a negative yearly SMB
        Tms = 16.81 - 0.00692*md.geometry.surface - 0.27937*np.abs(md.mesh.lat) + Tf
        Tms= Tms[0]

        # Calculate Accumulation perturbation with Tf forcing, Eqn (9)
        ACCdot = 2.5*2**((Tma+Tf)/10.) - 2.5*2**(Tma/10.)

        # Calculate Ablation, Eqn (10) (use for both Antarctica & Greenland), max melt is 10m/a
        ABLdot=0.*np.ones(md.mesh.numberofvertices)
        pos=np.nonzero(Tms>=0)
        ABLdot[pos]=np.minimum(1.4*Tms[pos],10)

        smb=ACCdot-ABLdot
        return smb[0]