'reinit'; 'set display color white'; 'c'

case='merra2'

if case='merra2'
ver='3.3.2a'
filedlw='merra2_lwga_1dy_1-1980-1-2019.nc'
vardlw='lwgab'
endif

***************************************************************
nh = 2
pathSL = '/homes/metofac/chepurin/GRADS/GS/'
pathFF = '/aosc/greenland/forcing/'case'/'

* nh - number harmonics to rectore seasonal cycle
***************************************************************

*---- open smoothed heat deficit in the mixed layer
'sdfopen 'case'_heat_dfct_10dy_smth.nc'

*---- set environment 
pathSL'set_file 1'; 'set z 1'
if ( case = 'merra2' )
  'set lon -181 181'
endif

*---- calculate annual harmonics from hdsm
pathSL'harmonics_annual hdsm 'nh' 36.525 hdnh'

*---- open forcing file for correction
'sdfopen 'pathFF''filedlw
pathSL'set_file 2'

*---- get time and time dependant functions
'set x 1'; 'set y 1'; 'set t 1'
'pi=atan2(1,1)*4'
'set t 1 last'; 'tn=tloop(gint(1,t=1,t+0)/(24*60))'
i=1
while(i<=nh)
  's'i'=sin(tn*2*'i'*pi/365.25)'
  'c'i'=cos(tn*2*'i'*pi/365.25)'
  i=i+1
endwhile

*---- interpolate harmonics to the forcing grid
pathSL'set_file 2'; 'set t 1'
'ai0 = lterp(a0,'vardlw')'
i=1
while(i<=nh)
  'ai'i' = lterp(a'i','vardlw')' 
  'bi'i' = lterp(b'i','vardlw')' 
  i=i+1
endwhile


*---- calculate correction for surface downwards thermal heat flux
'set t 1 last'
'defct = ai0'
i=1
while(i<=nh)
   'defct = defct + (ai'i'*c'i'+bi'i'*s'i')'
   i=i+1
endwhile

*---- reduce maximum correction between -60 and 60 W/m^2
pathSL'min_xy defct  60 defct'
pathSL'max_xy defct -60 defct'

*---- cut correction by 65 north
'defct = maskout(defct,65-lat)'

*---- fill missing values with 0.0
'defct = const(defct,0.0,-u)'

*---- save correction for surface downwards thermal heat flux
'set sdfwrite -rt -flt 'case'_'vardlw'_correction_1dy_smth_new.nc'
'set sdfattr defct String long_name heat flux deficit in mixed layer'
'set sdfattr defct String float actual_range -60 60 W/m^2'
'set sdfattr defct String unit W/m^2'
'sdfwrite defct'
'clear sdfwrite'

*---- correct surface downwards thermal heat flux
' 'vardlw' = 'vardlw'+defct'

*---- save corrected downward long wave radiation flux 
'set sdfwrite -rt -flt 'case'_'vardlw'_corrected_1980-2019_new.nc'
*'set sdfattr lon String standard_name longitude'
*'set sdfattr lon String units degree_north'
*'set sdfattr lon String axis X'
*'set sdfattr lat String standard_name latitude'
*'set sdfattr lat String units degree_east'
*'set sdfattr lat String axis Y'
*'set sdfattr time String calendar julian'
*'set sdfattr time String axis T'
'set sdfattr 'vardlw' String long_name 'case' corrected downward long wave radiation flux (W/m^2)'
'sdfwrite 'vardlw
'clear sdfwrite'