**************************************************************
*Calculates and displays theta-e in a pre-defined
*MUST SET THE DIMENSION ENVIRONMENT FIRST!!!
*This script will also compute surface theta-e
*
*The variables will be:
*3-D theta-e     : thetae1
*Surface theta-e : thetae2m
*************************************************************

*This series of code is just prep for a status message that will be output after the computations are finished
'q dims'
lonline = sublin(result,2)
latline = sublin(result,3)
levline = sublin(result,4)
timeline = sublin(result,5)
lon1 = subwrd(lonline,6)
lon2 = subwrd(lonline,8)
lat1 = subwrd(latline,6)
lat2 = subwrd(latline,8)
lev1 = subwrd(levline,6)
lev2 = subwrd(levline,8)
timecheck = subwrd(timeline,3)

if (timecheck = 'fixed')
 timestep = subwrd(timeline,9)
else
 time1 = subwrd(timeline,11)
 time2 = subwrd(timeline,13)
endif

say 'Is this a WRF simulation or observed data? (type 1 for WRF simulation and anything else for observed)'
pull ans
say 'Type SFC if you want only surface theta-e, otherwise, type anything else (or nothing)'
pull surface

*************The actual computation starts here***************

Rv = 461.5
Rd = 287.05
Cp = 1005
lv = 2.5e6

if (ans != 1)
* First section is for variable names that come with observed data and the like
 if (surface != 'SFC')
* 3-D computation
  temp_K = "tmpprs"
* If qvapor is not in the binary file, then you must compute it from dewpoint.  Otherwise just use qvapor.
  'es = 6.11*exp((2.5e6/461.5)*((1/273) - (1/'temp_K')))'
  'e = rhprs*es/100'
  'mixratio = 0.622*e/(lev - e)'
  'equivt = 'temp_K' + (2.5e6*(mixratio/1005))'
  'thetae = equivt*pow(1000/lev,287.05/1005)'
 endif

* ------------------------------------------------------------

* Surface computation
 sfc_temp_K = "tmp2m"
 sfc_pressure = "(pressfc/100)"
* Most observed datasets have 2m dewpoint, so you can skip a few steps unlike above...
 'e2 = 6.11*exp((2.5e6/461.5)*((1/273) - (1/(dpt2m))))'
 'sfcmixratio = 0.622*e2/('sfc_pressure' - e2)'
 'equivt2m = 'sfc_temp_K' + (2.5e6*(sfcmixratio/1005))'
 'thetae2m = equivt2m*pow(1000/'sfc_pressure',287.05/1005)'

*-----------------------------------------------------------------
*=================================================================
*-----------------------------------------------------------------
else

* For a simulation (WRF variables)
 if (surface != 'SFC')
* 3-D
  'tempK = tc + 273'
  'es = 6.11*exp((2.5e6/461.5)*((1/273) - (1/tempK)))'
  'e = qvapor*((p+pb)/100) / (qvapor + 0.622)'
  'equivt = tempK + 'Lv'*qvapor/'Cp'
  'thetae = equivt*pow(1000/slp,'Rd'/'Cp')'
 endif

* Surface
 'equivt2m = t2 + (2.5e6*(q2(z=1)/1005))'
 'thetae2m = equivt2m*pow(1000/(psfc(z=1)/100),287.05/1005)'

endif
*******************************************************************

*Printing status message to indicate successful computation
say 'Theta-e calculations complete!'
say 'Calculations were computed between...'
say 'Longitude: 'lon1' and 'lon2
say 'Latitude: 'lat1' and 'lat2
say 'Pressure levels 'lev1' and 'lev2
say 'At time(s) 'time1' and 'time2' or 'timestep