;-------------------------------
undef("eady_growth_rate")
function eady_growth_rate(th:numeric, u:numeric, z:numeric, lat:numeric,  opt[1]:integer, ndimz[1]:integer)
;
; Maximum Eady Growth Rate  
;
; Reference: 
;   R. S. Lindzen and Brian Farrell, 1980: 
;   A Simple Approximate Result for the Maximum Growth Rate of Baroclinic Instabilities. 
;   J. Atmos. Sci., 37, 1648–1654. 
;    http://dx.doi.org/10.1175/1520-0469(1980)037<1648:ASARFT>2.0.CO;2
;
;   Simmonds, I., and E.-P. Lim (2009): 
;   Biases in the calculation of Southern Hemisphere mean baroclinic eddy growth rate 
;   Geophys. Res. Lett., 36, L01707 
;   doi:10.1029/2008GL036320
;
;   Vallis, G.K. (2006) 
;   Atmospheric and Oceanic Dynamics: Fundamentals and Large-Scale Circulation
;   Cambridge Univ. Press, New York.
;    
; Nomenclature
; th   - potential temperature (K)
;        http://glossary.ametsoc.org/wiki/potential_temperature
;        https://www.ncl.ucar.edu/Document/Functions/Contributed/pot_temp.shtml
; u    - zonal wind components (m/s) ; same dimensionality as 'th'
; z    - height (m)
; lat  - latitude of each grid point ; same dimensionality as 'th'  
; opt  - =0 (egr only); 
; ndimz- dimension of 'th' for which vertical gradient is to be calculated
;        th(:), wspd(:), z(:) ..... ndimz=0
;        th(:,:,:)  , z(:) or z(:,:,:)   and the left is the vertical dim, ndimz=0
;        th(:,:,:,:), z(:) or z(:,:,:,:) and the left is the time dimension 
;                                        and the next is height,  ndimz=1 
;
;
local dimth, dimu, dimz, dimth, dimlat, rankth, ranku, rankz, rankth, ranklat \
    , u_shear, brunt, fcor, g, omega, con, rad, egr
begin
                                ; dimension checking
  dimu    = dimsizes(u)
  dimz    = dimsizes(z)
  dimth   = dimsizes(th)
  dimlat  = dimsizes(lat)

  ranku   = dimsizes(dimu)   
  rankz   = dimsizes(dimz)   
  rankth  = dimsizes(dimth)
  ranklat = dimsizes(dimlat)   

  if (.not.(rankth.eq.ranku .and. rankth.eq.rankz .and. rankth.eq.ranklat) \
      .and. all(dimth.eq.dimu) .and. all(dimth.eq.dimz) .and. all(dimth.eq.dimlat)) then
      print("eady_growth_rate: th, u, lat and must have the same rank & sizes")
      print("ranku="+ranku+"  rankth="+rankth+"  rankz="+rankz+"  ranklat="+ranklat)
      print(dimth)
      print(dimu) 
      print(dimz) 
      print(dimlat) 
      exit
  end if

  brunt   = brunt_vaisala_atm(th, z, 0, ndimz)

  dudz    = center_finite_diff_n(u , z, False, 0, ndimz) 

  copy_VarCoords(u, dudz)
  dudz@long_name  = "vertical gradient of the zonal wind (zonal wind shear): du/dz"
  dudz@units      = "1/s"
 ;print("eady_growth_rate: dudz: min="+min(dudz)+"   max="+max(dudz))

  g     = 9.80665                    ; m/s2 ; gravity at 45 deg lat used by the WMO
  con   = 0.3098
  fcor  = coriolis_param(lat)        ; (1/s)     ; coriolis parameter
                                     ; prevent 1/0
  if (any(brunt.eq.0)) then
      if (.not.isatt(brunt,"_FillValue")) then
          if (typeof(brunt).eq."double") then
              brunt@_FillValue = 1d20
          else
              brunt@_FillValue = 1e20
          end if
       end if
       brunt = where(brunt.eq.0, brunt@_FillValue, brunt)
  end if

  egr   = con*abs(fcor)*abs(dudz)/brunt
  egr@long_name = "maximum eady growth rate"
  egr@units     = ""
  copy_VarCoords(u, egr)

  if (opt.eq.0) then
      return (egr)
  else if (opt.eq.1) then
      return ( [/egr, dudz /] )
  else if (opt.eq.2) then
      return ( [/egr, dudz, brunt /] )
  end if
  end if
  end if

end