'reinit'

rc=gsfallow('on')
define()

say '>>> Calculando os 4 termos da frontogenese em thetae... '

** Os novos arquivos nao dependem do mes nem do ano**

xmax=_eta_nx
ymax=_eta_ny

nlevs=1
lev.1=850

lv=1
while (lv<=nlevs)
  _the.lv=openvar('the',lev.lv)
  _temp.lv=openvar('temp',lev.lv)
  _uvel.lv=openvar('uvel',lev.lv)
  _vvel.lv=openvar('vvel',lev.lv)
  _omeg.lv=openvar('omeg',lev.lv)

  lv=lv+1
endwhile
* para calcular FG4, precisa da deriva na vertical, ou seja, the 1 nivel acima e abaixo de lev.1
_thedown=openvar('the',900)
_theup=openvar('the',800)

* ajusta as dimensoes'
'set x 0 'xmax+1
'set y 0 'ymax+1
'set z 1'
'set t 1'
'q time'; rc=subwrd(result,3); say 'tempo 1='rc

* Elementos de distancia e f para calculo do vento geostrofico
'define dx=6.37e6*cdiff(lon,x)*(acos(-1)/180.)*cos(lat*acos(-1)/180.)'
'define dy=6.37e6*cdiff(lat,y)*(acos(-1)/180.)'
'define f=2*7.292e-5*sin(lat*acos(-1)/180.)'

*******************************************************************
** Os derivados serao gravados em um unico arquivo
*******************************************************************
base=_p1'_front_'_p2''_eta_begy''_p3''
fout=_dir'/'base'.bin'
'set gxout fwrite'
'set fwrite -le -st 'fout
say 'gravando no arquivo: '
say '   'fout

*******************************************************************
say 'Loop nos meses...'
*******************************************************************

tbeg=1
a=_eta_begy
m=_eta_begm
while( (a*100+m)<=_eta_end )
  
** ajusta o tempo para cobrir este mes
  tfim=tbeg+30*4-1
  'set t 'tbeg' 'tfim
  'q dims'; lin=sublin(result,5)
  timei=subwrd(lin,6)
  timef=subwrd(lin,8)
  say 'Tempo REAL: year='a' month='m
  say 'Tempo FAKE: timei= 'timei' timef= 'timef

******************************************************************
*** fazer as contas
******************************************************************
  'set x 0 'xmax+1
  'set y 0 'ymax+1

  lv=1
  while(lv<=nlevs)
    say lev.lv
    'set lev 'lev.lv
***************************************************************************
* frontogenese em temperatura potencial equivalente
    'define dthdx'lev.lv'=cdiff('_the.lv',x)/dx(z=1,t=1)'
    'define dthdy'lev.lv'=cdiff('_the.lv',y)/dy(z=1,t=1)'
    
    'define A=cdiff('_uvel.lv',x)/dx(z=1,t=1) - cdiff('_vvel.lv',y)/dy(z=1,t=1)'
    'define B=cdiff('_vvel.lv',x)/dx(z=1,t=1) + cdiff('_uvel.lv',y)/dy(z=1,t=1)'
    'define D=cdiff('_uvel.lv',x)/dx(z=1,t=1) + cdiff('_vvel.lv',y)/dy(z=1,t=1)'

***************************************************************************
* frontogenese por deformacao
    'FG3'lev.lv'=((pow(dthdx'lev.lv',2)-pow(dthdy'lev.lv',2))*A + 2*dthdx'lev.lv'*dthdy'lev.lv'*B)/(-2*mag(dthdx'lev.lv',dthdy'lev.lv'))*86400.*1.e5'
    'FG2'lev.lv'=(pow(dthdx'lev.lv',2)+pow(dthdy'lev.lv',2))*D/(-2*mag(dthdx'lev.lv',dthdy'lev.lv'))*86400.*1.e5'

    'undefine A'
    'undefine B'
    'undefine D'
    if (lev.lv!='850')
      'undefine dthdx'lev.lv
      'undefine dthdy'lev.lv
    endif
***********************************
    lv=lv+1
  endwhile

* FG1 e FG4 apenas em 850mb
  'set lev 850'
  'define dthdt=('_the.1'(t+1)-'_the.1'(t-1))/12.0'
  'define dthdtdx=cdiff(dthdt,x)/dx(z=1,t=1)'
  'define dthdtdy=cdiff(dthdt,y)/dy(z=1,t=1)'
  'undefine dthdt'
  'FG1=(dthdx850*dthdtdx + dthdy850*dthdtdy)/mag(dthdx850,dthdy850)*24.*1.e5'
  'undefine dthdtdx'
  'undefine dthdtdy'

  'define wvel=-('_omeg.1')*287.04*('_temp.1')/(lev*100)/9.8'
  'define dwdx=cdiff(wvel,x)/dx(z=1,t=1)'
  'define dwdy=cdiff(wvel,y)/dy(z=1,t=1)'
  'undefine wvel'
  'define dthdp=('_thedown'(lev=900)-'_theup'(lev=800))/(90000-80000)'
  'FG4=dthdp*(dthdx850*dwdx+dthdy850*dwdy)/(-mag(dthdx850,dthdy850))*86400.*1.e5'
  'undefine dthdx850'
  'undefine dthdy850'
  'undefine dwdx'
  'undefine dwdy'
  'undefine dthdp'

*******************************************************************
**** Grava no arquivo
*******************************************************************
  'set x 1 'xmax
  'set y 1 'ymax

* Obtem o numero de tempos no dia
  'q dims'
  lin = sublin(result,5)
  first= subwrd(lin,11)
  last = subwrd(lin,13)
  say '    first time = 'first
  say '    last time  = 'last

*  Loop para escrever no arquivo de saida
  t=first
  while(t<=last)
    'set t 't
      
    'set lev 850'
    'd const(FG1,9.999e20,-u)' 
      
    lv=1
    while(lv<=nlevs)
      'set lev 'lev.lv
      'd const(FG2'lev.lv',9.999e20,-u)'
      lv=lv+1
    endwhile

    lv=1
    while(lv<=nlevs)
      'set lev 'lev.lv
      'd const(FG3'lev.lv',9.999e20,-u)'
      lv=lv+1
    endwhile

    'set lev 850'
    'd const(FG4,9.999e20,-u)' 
    
    t=t+1
  endwhile

***************************************************************************
* Libera Memoria
***************************************************************************
  'undefine fg1'
  'undefine fg4'
  lv=1
  while(lv<=nlevs)
    'undefine FG2'lev.lv
    'undefine FG3'lev.lv
    lv=lv+1
  endwhile

  m=m+1
  if (m>12)
    a=a+1
    m=1
  endif
  tbeg=tbeg+30*4

endwhile 

'disable fwrite'
'set gxout contour'

*******************************************************************
**** GERA O CTL
*******************************************************************
nt=tbeg-1
fctl=_dir'/'base'.ctl'
list=''; lv=1
while(lv<=nlevs)
  list=list' 'lev.lv
  lv=lv+1
endwhile

nm=nomemes(_eta_begm)
rc=write(fctl,'dset ^'base'.bin')
rc=write(fctl,'undef  9.999E+20',append)
rc=write(fctl,'xdef  'xmax' linear '_eta_x0' '_eta_dx'',append)
rc=write(fctl,'ydef  'ymax' linear '_eta_y0' '_eta_dy'',append)
rc=write(fctl,'zdef 'nlevs' levels 'list,append)
rc=write(fctl,'tdef 'nt' linear 00Z01'nm''_eta_begy' 6hr',append)
rc=write(fctl,'vars 4',append)
rc=write(fctl,'fg1 0 99 ** frontogenese 1 por deformacao em thetae [K/100km/dia]',append)
rc=write(fctl,'fg2 'nlevs' 99 ** frontogenese 2 por deformacao em thetae [K/100km/dia]',append)
rc=write(fctl,'fg3 'nlevs' 99 ** frontogenese 3 por deformacao em thetae [K/100km/dia]',append)
rc=write(fctl,'fg4 0 99 ** frontogenese 4 por deformacao em thetae [K/100km/dia]',append)
rc=write(fctl,'ENDVARS',append)
rc=close(fctl)

*******************************************************************
**** GRIBA
*******************************************************************
*'!bin2grib_levs 'fctl

'reinit'
out=_dir'/'_p1'_fg1_'_p2''_eta_begy''_p3'_850hpa_grib'
'!lats4d.sh -i 'fctl' -o 'out' -levs 850 -vars fg1 -format grads_grib -ftype ctl -mxtimes 50000'

i=2
while(i<=3)
  lv=1
  while(lv<=nlevs)
    'reinit'
    out=_dir'/'_p1'_fg'i'_'_p2''_eta_begy''_p3'_'lev.lv'hpa_grib'
    '!lats4d.sh -i 'fctl' -o 'out' -levs 'lev.lv' -vars fg'i' -format grads_grib -ftype ctl -mxtimes 50000'
    lv=lv+1
  endwhile
  i=i+1
endwhile

'reinit'
out=_dir'/'_p1'_fg4_'_p2''_eta_begy''_p3'_850hpa_grib'
'!lats4d.sh -i 'fctl' -o 'out' -levs 850 -vars fg4 -format grads_grib -ftype ctl -mxtimes 50000'

'!rm 'fout
'!rm 'fctl

* fim

* funcao que calcula o numero de dias no mes.
function fimdomes(m,a)

if (m=1|m=3|m=5|m=7|m=8|m=10|m=12)
  fimdomes=31
endif

if (m=4|m=6|m=9|m=11)
  fimdomes=30
endif

if (m=2)
  if (math_mod(a,4)=0)
    fimdomes=29
  else 
    fimdomes=28
  endif
endif  

return(fimdomes)

*fim     

* funcao que atribui o nome do mes com tres letras.
function nomemes (m)

if (m=1);  nomemes="jan"; endif
if (m=2);  nomemes="feb"; endif
if (m=3);  nomemes="mar"; endif
if (m=4);  nomemes="apr"; endif
if (m=5);  nomemes="may"; endif
if (m=6);  nomemes="jun"; endif
if (m=7);  nomemes="jul"; endif
if (m=8);  nomemes="aug"; endif
if (m=9);  nomemes="sep"; endif
if (m=10);  nomemes="oct"; endif
if (m=11);  nomemes="nov"; endif
if (m=12);  nomemes="dec"; endif

return(nomemes)



*fim