Dim pApp As IApplication
    Set pApp = New AppRef
    Dim pMxDoc As IMxDocument
    Dim pMap As IMap
    Set pMxDoc = pApp.Document
    Set pMap = pMxDoc.FocusMap
    
    'Get the input raster from the first layer in ArcMap
    Dim pLayer As ILayer
    Set pLayer = pMap.Layer(0)
    If Not TypeOf pLayer Is IRasterLayer Then
        Exit Sub
    End If
    Dim pLayer2 As ILayer
    Dim RasterUID As UID
    Set RasterUID = New UID
    RasterUID = "{6CA416B1-E160-11D2-9F4E-00C04F6BC78E}"
    Set pLayer2 = FindLayer(RasterUID, "dem")
    If Not TypeOf pLayer Is IRasterLayer Then
        Exit Sub
    End If
    
    Dim pRasLayer As IRasterLayer
    Set pRasLayer = pLayer
    Dim SlopeRaster As IRaster
    Set SlopeRaster = pRasLayer.Raster
    Dim pRasLayer2 As IRasterLayer
    Set pRasLayer2 = pLayer2
    Dim DEMraster As IRaster
    Set DEMraster = pRasLayer2.Raster
    
    'Create a RasterModel object
    Dim pRasModel As IRasterModel
    Set pRasModel = New RasterModel
    
    'Set output workspace in the analysis environment
    Dim pEnv As IRasterAnalysisEnvironment
    Set pEnv = pRasModel
    
    Dim pWS As IWorkspace
    Dim pWSF As IWorkspaceFactory
    Set pWSF = New RasterWorkspaceFactory
    Set pWS = pWSF.OpenFromFile("c:\usertemp", 0)
    Set pEnv.OutWorkspace = pWS

    'Bind the input raster
    pRasModel.BindRaster SlopeRaster, "slope"
    pRasModel.BindRaster DEMraster, "dem"

    'Specify map algebra expression, use vbLf to create separate lines
    pRasModel.Script = "[out1] = [dem] / 12 * 39.37" + vbLf + _
                       "[out2] = ([slope] < 3.0 and [dem] < 400) > 0"
    'Execute map algebra expression(s)
    pRasModel.Execute
    
    'Get output rasters
    Dim pOutRaster1 As IRaster
    Set pOutRaster1 = pRasModel.BoundRaster("out1")
    Dim pOutRaster2 As IRaster
    Set pOutRaster2 = pRasModel.BoundRaster("out2")
    
    'Unbind input raster
    pRasModel.UnbindSymbol "input1"
    pRasModel.UnbindSymbol "input2"
    
    Dim pRasDS01 As IRasterDataset
'    Set pRasDS01 = OpenGridAsIRasterDataset(sGridPath, sGridName)
    'Add outputs into ArcMap as raster layers
    Dim pOutRasLayer1 As IRasterLayer
    Set pOutRasLayer1 = New RasterLayer
    pOutRasLayer1.CreateFromRaster pOutRaster1
    pMap.AddLayer pOutRasLayer1
    Dim pOutRasLayer2 As IRasterLayer
    Set pOutRasLayer2 = New RasterLayer
    pOutRasLayer2.CreateFromRaster pOutRaster2
    pMap.AddLayer pOutRasLayer2