f = addfile("mslp_jan_2017.nc","r")
 p = f->slp           ; ordered North-to-South

 if (p&lat(0) .gt. p&lat(1)) then
     p = p(:,::-1,:)      ; ordered South-to-North [required by spherical harmonics]
     print("pressure array reordered South-to-North")
     print("----------")
 end if

 nmsg = num(ismissing(p))   ; count missing values 
 if (nmsg.gt.0) then
     print("EXIT: missing values encountered: nmsg="+nmsg)
     exit
 end if

 printVarSummary(p)
 printMinMax(p,0)
 print("----------")

; Compute Laplacian

;lap = lapsF(p)       ; Laplacian function
;copy_VarCoords(p,lap)
;lap@long_name = "Laplacian(p)"
;lap@units     = "..."
;printVarSummary(lap)
;printMinMax(lap,0)
;print("----------")

; gradsf is a procedure which requires that the return arrays be preallocated
; create arrays to hold output, same size and type as input

 dimp = dimsizes(p)
 dpdx = new(dimp, typeof(p), p@_FillValue) ; dp/dlon => dp/dx
 dpdy = new(dimp, typeof(p), p@_FillValue) ; dp/dlat => dp/dy 

; procedure gradsf will return dp/dx, dp/dy in the preallocated arrays

 gradsf (p,  dpdx, dpdy) 

 work = new(dimp, typeof(p), p@_FillValue) ; placeholder for temporary values
 px2  = new(dimp, typeof(p), p@_FillValue) ; d2(px)/dx2 
 py2  = new(dimp, typeof(p), p@_FillValue) ; d2(py)/dy2 

 gradsf (dpdx,  px2, work) 
 gradsf (dpdy, work, py2 ) 

 px2@long_name = "d2(p)/d2lon"
 px2@units = "..."
 copy_VarCoords(p, px2)

 py2@long_name = "d2(p)/d2lat"
 py2@units = "..."
 copy_VarCoords(p, py2)
 print("----------")

 printVarSummary(px2)
 print("----------")
 printVarSummary(py2)
 print("----------")
 
;=======================

 px2py2 = px2*py2
 px2py2@long_name = "product: [d2(p)/d2lon]*[d2(p)/d2lat]"
 px2py2@units     = "..."
 copy_VarCoords(p, px2py2)
 printVarSummary(px2py2)