#########################################################
#
# grads_sample.gs
#
# sample script for GrADS
#
# draw 4-panel contour map of
# climatological mean air temperature (NCEP reanalysis)
# at 500 hPa for JFM, AMJ, JAS, and OND
#
#########################################################

######################################################
FUNCTION MAIN()
######################################################

'reinit'
#'set imprun imprun.sc'
'c'

##################################
# definitions
##################################

MDIR    = '/home/reichler/NNR_monthly/'
SDIR    = 'pressure/'

VAR     = 'air'
CLEVS   = '-40 -35 -30 -25 -20 -15 -10 -5 0'
UNIT    = 'C'

STAT    = 'ltm'
LEVEL   = 500
LON1    = 0
LON2    = 360
LAT1    = -90
LAT2    = 90

PRINT   = 0      
PFILE   = VAR '_' STAT '_' LEVEL

#################################
# page layout, color, map proj., etc.
#################################

# number of panels in x and y
NX=2
NY=2

# all units are inches
# landscape
'set vpage 0.0 11.0 0.0 8.5'
# portrait
#'set vpage 0.0 8.5 0.0 11.0'

# size of individual panels: landscape
X1=0.5
DX=4.4
Y1=4.5
DY=3.3

# thickness of text lines
TTHICK=3

# main title
XPOS=5.5
YPOS=8.2
'set string 1 c 'TTHICK' 0'
'set strsiz 0.2'
'draw string ' XPOS ' ' YPOS '    ' VAR ' ' LEVEL ' hPa ' 

#'set mproj latlon'
'set mproj scaled'
#'set mproj mollweide'
'set mpdset lowres'
'set mpdraw off'
#'set grid off'
'set grid on'

# open file
FILE    = MDIR SDIR VAR '.mon.' STAT '.nc'
say     FILE
'sdfopen ' FILE

#################################
# plot loop
#################################

PY=0
WHILE(PY<NY)
#############

PX=0
WHILE(PX<NX)
#############

T=PY*NX+PX+1
say T

T1=(T-1)*3+1
T2=T1+2

X=X1+PX*DX*1.28
Y=Y1-PY*DY*1.2
'set parea ' X ' ' X+DX ' ' Y ' ' Y+DY

# select plot dimensions
'set lev ' LEVEL
'set lon ' LON1 ' ' LON2
'set lat ' LAT1 ' ' LAT2
'set t ' 1

# shading
'set gxout shaded'
'set clevs ' CLEVS
'd ave(' VAR ',t='T1',t='T2')'

# contours
'set gxout contour'
'set clevs ' CLEVS
'set ccols 1 1 1 1 1 1 1 1 1 1 1 1 1'
#'set clab off'
'set cstyle 1'
'd ave(' VAR ',t='T1',t='T2')'

# title
XPOS=X+DX/2
YPOS=Y+DY*1.04
'set string 1 c 'TTHICK' 0'
'set strsiz 0.15'
'draw string ' XPOS ' ' YPOS '    ' T1'-'T2

# colorbar
CX=X+DX*1.05
CY=Y+DY/2*0.9
'run cbarn 0.4 1 ' CX ' ' CY
XPOS = CX
YPOS = Y+DY*0.85
'set string 1 c 'TTHICK' 0'
'set strsiz 0.1'
'draw string ' XPOS ' ' YPOS ' [' UNIT ']'

# y label
XPOS=X-0.11*DX
YPOS=Y+DY/2
'set string 1 c 'TTHICK' 90'
'set strsiz 0.08'
S='    lat. [deg.]'
'draw string ' XPOS ' ' YPOS S

# x label
XPOS=X+DX/2
YPOS=Y-DY*0.09
'set string 1 c 'TTHICK' 0'
'set strsiz 0.08'
S='    lon. [deg.]'
'draw string ' XPOS ' ' YPOS S

'set map 16 1' 
'set mpdraw on'
'draw map'

PX=PX+1
ENDWHILE
#############

PY=PY+1
ENDWHILE
#############

#################################
# make a hardcopy
#################################

if PRINT=1 
  line = 'PRINTING to ' PFILE
  say line
  '!rm ' PFILE '.mta'
  'enable print ' PFILE '.mta'
  'print'
  'disable print'
  '!rm ' PFILE '.ps'
  '!gxps -i ' PFILE '.mta -o ' PFILE '.ps'
  '!lpr -P cirp-bw ' PFILE '.ps'
endif

#################################
# end of script
#################################

RETURN