PRO Hydrostatic

; ===============================================================================
; Program to print vertical pressure gradients and density *g profiles
; ===============================================================================

; ----------------------------------------------------------------------
; Set-up the colours and the postscript output
; ----------------------------------------------------------------------
PRINT, 'Setting-up the colours and the postscript device'
r = [0,255,0,0,0,0,66,184,245,255,255,255,230,191,148,105,161,184]
g = [0,255,84,199,255,255,255,255,255,209,135,28,0,0,0,0,0,0]
b = [0,255,255,255,255,178,0,0,0,0,0,0,0,0,0,20,135,186]
TVLCT, r,g,b


SET_PLOT, 'ps'
PostscriptFile = 'Hydrostatic.ps'
DEVICE, FILE=PostscriptFile, /PORTRAIT, /COLOR, XSIZE=20.0, YSIZE=30.0, $
        XOFFSET=2.0, YOFFSET=2.0, FONT_SIZE=8, /ENCAPSULATED

!P.POSITION = [0.15, 0.05, 0.95, 0.95]

x_pos = 0       ; Greenwich
y_pos = 171     ; UK
;y_pos = 109     ; Near equator

; ----------------------------------------------------------------------
; Read-in the wind, pressure and orography data
; ----------------------------------------------------------------------
PRINT, 'Reading in the data'
fileid     = NCDF_OPEN('Fc.nc')
varidx     = NCDF_VARID(fileid, 'longitude_1')
varidy     = NCDF_VARID(fileid, 'latitude')
varidz     = NCDF_VARID(fileid, 'hybrid_ht')
varidrho   = NCDF_VARID(fileid, 'unspecified') ; density
varidp     = NCDF_VARID(fileid, 'field7')      ; pressure
varidorog  = NCDF_VARID(fileid, 'ht')          ; orography

NCDF_VARGET, fileid, varidx, x_degrees
NCDF_VARGET, fileid, varidy, y_degrees
NCDF_VARGET, fileid, varidz, vertlevs
NCDF_VARGET, fileid, varidrho, density
NCDF_VARGET, fileid, varidp, press
NCDF_VARGET, fileid, varidorog, orog
NCDF_CLOSE, fileid


; ----------------------------------------------------------------------
; Calculate the height field (fn of long and lat) from sea level
; ----------------------------------------------------------------------
PRINT, 'Calculating the height field at position ', x_degrees(x_pos), y_degrees(y_pos)
eta                 = FLTARR(38)
ztop                = vertlevs(37)
eta(0:37)           = vertlevs(0:37) / ztop
interface           = 30-1
height              = FLTARR(38)
height(0:37)        = eta(0:37) * ztop
height(0:interface) = height(0:interface) +                       $
                      (1.0 - eta(0:interface) / eta(interface)) * $
                      (1.0 - eta(0:interface) / eta(interface)) * $
                      orog(x_pos,y_pos)

; ----------------------------------------------------------------------
; Calculate -dp/dz, and pure density * g
; ----------------------------------------------------------------------
PRINT, 'Calculating dp/dz'

; -dp/dz
dpdz    = FLTARR(37)
height1 = FLTARR(37)
FOR z = 0, 36 DO BEGIN
  dpdz(z)    = -1.0 * ( press(x_pos,y_pos,z+1) - press(x_pos,y_pos,z) ) / $
                      ( height(z+1) - height(z) )
  height1(z) = ( height(z+1) + height(z) ) / 2.0
ENDFOR

; rho * g
rhog     = FLTARR(38)
a        = 6371000.0
g        = 9.8
err      = FLTARR(38)
FOR z = 0, 37 DO BEGIN
  rhog(z) = g * density(x_pos,y_pos,z) / ( (height(z) + a) * (height(z) + a) )
ENDFOR

; Relative error between these quantities
FOR z = 0, 36 DO BEGIN
  err(z)  = dpdz(z) - (rhog(z) + rhog(z+1)) / 2.0
ENDFOR


; ----------------------------------------------------------------------
; Plot fields that should be in hydrostatic balance
; ----------------------------------------------------------------------
PRINT, 'Plot fields that should be in hydrostatic balance'

;PLOT, dpdz, height1, /NOERASE


PLOT, rhog, height, XSTYLE=5, YTITLE='height (m)', /NOERASE
AXIS, XAXIS=0, XTITLE='density * g'

PLOT, err, height1, XSTYLE=5, YSTYLE=5, COLOR=13, /NOERASE
AXIS, XAXIS=1, XTITLE='hydrostatic residual', COLOR=13

DEVICE, /CLOSE
END