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

case='merra2'

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

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

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

*---- open smoothed fresh water deficit in the mixed layer
'sdfopen 'case'_rain_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 sdsm
pathSL'harmonics_annual sdsm 'nh' 36.525 sdnh'

*---- 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 fresh water 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 rainocn
'set sdfwrite -rt -flt 'case'_rain_correction_1dy_smth_new.nc'
'set sdfattr defct String long_name fresh water flux deficit in mixed layer'
'set sdfattr defct String unit kg/m^2/s'
'sdfwrite defct'
'clear sdfwrite'

*---- correct surface downwards thermal heat flux
**** old way
*'dum = 'vardlw'-defct'
*'dum = maskout(dum,dum)'
*'cor'vardlw' = const(dum,0.0,-u)'
**** new way
'alpha = defct/ave('vardlw',t-15,t+15)'
say 'alpha was calculated'
*---- make  -1 < alpha < 1
pathSL'min_xy alpha  1 alpha'
say 'alpha maximum was seted up'
pathSL'max_xy alpha -9 alpha'
say 'alpha minimum was seted up'
'dum = 'vardlw'*(1-alpha)'; 'undefine alpha'
say 'corrected value was calculated'
'dum1= maskout(dum,dum)'; 'undefine dum'
say 'negative values were remuved'
'cor'vardlw' = const(dum1,0.0,-u)'; 'undefine dum1'
say 'missed values were seted up to 0.0'
*'cor'vardlw' = const(maskout(dum,dum),0.0,-u)'

*---- save corrected downward long wave radiation flux 
'set sdfwrite -rt -flt 'case'_rainocn_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 cor'vardlw' String long_name corrected daily liquid precipitation ('vardlw') [kg/m^2/S]'
'sdfwrite cor'vardlw
'clear sdfwrite'