#############################################################################
#Script to plot wind chill.                                                 #
#Note-uses modified version of a colorbar cbarm.gs                          #
#Martynas K. 2006-11-09                                                     #
#############################################################################
cat >> ${takas3}./windchill.gs <<EOF
'open ${file1}.ctl'
'q dims'
initial=subwrd(result,46)
_times=getinfo()
valandos=0
timestep=1
while (timestep <= _times)
'set datawarn off'
'set grads off'
'set grid off'
'run ${takas3}map.gs 1 1 6'
'set t ' timestep
'q dims'
name=subwrd(result,46)
'run ${takas3}colorgen.gs'

*Here you calculate the wind chill

'define temp=0.6215*(t2-273.16)'
'define wind=mag(u10,v10)*3.6'
'define wp=pow(wind,0.16)'
'define chill=13.12+temp-(11.37*wp)+(0.3965*temp*wp)'

* Plotting shaded wind chill, levels -40 to +10, interval 1.

'set gxout shaded'
'run ${takas3}colorgen.gs'
'set clevs -40 -39 -38 -37 -36 -35 -34 -33 -32 -31 -30 -29 -28 -27 -26 -25 -24 -23 -22 -21 -20 -19  -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10'
'set ccols 74 73 72 71 70 69 68 67 66 65 64 63 62 61 60 59 58 57 56 36 35 34 33 32 31 30 29 28 27 26 25 24 23 22 21 20 19 18 55 54 53 52 51 50 49 48 47 46 45 43 42 41 40 39 38 37'
'd chill'
'run ${takas4}cbarm.gs'


*Plotting contoured wind chill

'set gxout contour'
'set clevs -40 -38 -36 -34 -32 -30 -28 -26 -24 -22 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10'
'set rgb 17 191 191 191'
'set ccols  17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17'
'set cthick 5'
'set cstyle 1'
'set clab on'
'd chill'


*Masking temperatures that are higher than +10 degrees Celcius -
*Assume that the grey filled areas equal the actual 2 meter temperature
*
'define tp=temp'
'define tp=temp+40'
'define tmp=maskout(temp,tp-50)'
'set rbcols 16'
'd tmp'

*Masking wind speeds that are lower than 5 km/hr.
*Assume that the grey filled areas equal the actual 2 meter temperature
*

'set gxout shaded'
'set rbcols 16'
'define lowspeed=wind*(-1)'
'define lowspeed=lowspeed+150'
'define wnd=maskout(wind,lowspeed-145)'
'd wnd'


* Other stuff, printing, etc...
*Not related with wind chill directly.



'set clopts 1 1 0.1'
'run ${takas1}title.gs Windchill 'initial' 'valandos' 'name''
'exec ${takas4}regulars. reglon reglat 0 -30'
'run ${takas4}platlon reglon reglat 1 1 3 10 10 off'
'run ${takas4}ispejimas.gs'
'run ${takas1}putmark1.gs 3 0.05'
'printim ${takas2}windchill+'valandos'.png white'
'clear'

timestep=timestep+1
valandos=valandos+3
endwhile
*funkcija laiko zingsniu skaiciui
function getinfo()
"set gxout stat"
"set stnprint on"
"q file"
res=result
time=sublin(res,5)
times=subwrd(time,12)
*say times
return(times)
EOF