;***********************************************************************;
; Procedure : gsn_geop_hgt                                              ;
;                                                                       ;
; Returns geopotential height (in km) given array p (pressure in mb)    ;
; p must lie between 1013.25 mb and 2.54e-06 mb.                        ;
;                                                                       ;
; Algorithm is simply logarithmic interpolation from Standard           ;
; Atmosphere.                                                           ;
; Intended to provide an ESTIMATE of geopotential height when           ;
; temperature and geopotential data are lacking.                        ;
;                                                                       ;
; Values above 85km were obtained from the COSPAR International         ;
; Reference Atmosphere: 1986 QB495 .A38 v.10 No.12  [FL]                ;
;                                                                       ;
;***********************************************************************;
 
undef("gsn_geop_hgt")
function gsn_geop_hgt( p[*]:numeric )
local nsa,psa,zsa,ptmp,npres,found,i,j
 
begin
  if(isatt(p,"units").and.(any(p@units.eq.get_allowed_pres_units_pa())))
    ptmp = p * 0.01  ; Convert to mb
  else
    if((.not.isatt(p,"units")).or. \
       (isatt(p,"units").and. \
        .not.(any(p@units.eq.get_allowed_pres_units_hpa()).or. \
              any(p@units.eq.get_allowed_pres_units_mb()))))
      print("gsn_geop_hgt: Warning: The 'units' attribute is either not set, or it is not set")
      print("to the recognized names for 'hecto-pascals' or 'millibars', so")
      print("assuming pressure values are already converted to millibars.")
    end if
    ptmp = tofloat_wunits(p)     ; Assume already converted to mb!
  end if

  nsa = 53
  psa = new( (/nsa/), float, 1.e36 )
  zsa = new( (/nsa/), float, 1.e36 )
 
  zsa = (/ -0.3, \                      ; km
            0.0,  0.5,  1.0,  1.5,  2.0,  2.5,  3.0, \
            3.5,  4.0,  4.5,  5.0,  5.5,  6.0,  6.5, \
            7.0,  7.5,  8.0,  8.5,  9.0,  9.5, 10.0, \
           11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, \
           18.0, 19.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, \
           50.0, 60.0, 70.0, 80.0, 84.8, 87.7, 90.6, \
           93.3, 96.1, 97.5,100.4,104.9, \
          110.0,114.2,116.7,119.7/)

  psa = (/ 1050., \                     ; mb (hPa)
           1013.25, 954.61, 898.76, 845.59, 795.01, 746.91, 701.21,    \
           657.80, 616.60, 577.52, 540.48, 505.39, 472.17, 440.75,     \
           411.05, 382.99, 356.51, 331.54, 308.00, 285.84, 264.99,     \
           226.99, 193.99, 165.79, 141.70, 121.11, 103.52, 88.497,     \
           75.652, 64.674, 55.293, 25.492, 11.970, 5.746, 2.871, 1.491,\
           0.798, 0.220, 0.052, 0.010, 0.00485,0.00294,0.000178, \
           0.000108, 0.0000656, 0.0000511, 0.0000310, 0.0000146, \
           0.00000691, 0.00000419, 0.00000327, 0.00000254 /)


  if ( any(ptmp.lt.min(psa)) .or. any(ptmp.gt.max(psa))) then
    print("gsn_geop_hgt: Fatal: The pressure values do not fall between")
    print(min(psa) + " mb and " + max(psa) + " mb.")
    print("Execution halted.")
    exit
  end if
 
  npres = dimsizes(ptmp)
  gph = new(npres,float)

  do i = 0,npres-1
    found = False
    j = 0

    do while(.not.found.and.j.le.nsa-2)
     if ( ( ptmp(i) .le. psa(j) ) .and. ( ptmp(i) .ge. psa(j+1) ) ) then
       gph(i) = zsa(j) + (zsa(j+1) - zsa(j)) * \
                log( psa(j)/ptmp(i) )/log( psa(j)/psa(j+1) )
       found = True 
     end if
     j = j + 1
    end do

  end do
 
  delete(psa)
  delete(zsa)
  delete(ptmp)
 
  return(gph)
 
end