;This scripts calculates from 03 hourly data the daily minimum values of temperature and detect heatwaves
;Indice:Indice 1
;Data: ERA-I
;Resolution:1*1
;Parameter:Tn
;Minimum duration: 04 days

T=addfile("/research/geog/data2/DATA/kg312/ERA-I/Tmax_Tmin/ERA-I_1_1_1979-2017_Tmax_Tmin_NTA.nc","r")
mn2t=T->mn2t
Tmin=mn2t*mn2t@scale_factor+mn2t@add_offset
Tmin=Tmin-273.15
copy_VarCoords(mn2t,Tmin)
mx2t=T->mx2t

;================================================;Computation of the daily minima of the dry bulb temperature================================================================================================
Tmin_min=calculate_daily_values(Tmin,"min",0,False)
Tn:=Tmin_min
n:=dimsizes(Tn)
Tn:=Tn(0:n(0)-2,:,:)

;;;;;Spatiotemporal attributes

time:=Tn&time
calendar=cd_calendar(time,0)
yyyy:=floattointeger(calendar(:,0))
mm:=floattointeger(calendar(:,1))
day:=floattointeger(calendar(:,2))
yyyymm:=tointeger(cd_calendar(time,1))
yyyymmdd:=tointeger(cd_calendar(time,2))
yyyyfrac=cd_calendar(time,4)     ;fraction of year
ddd:=floattointeger(day_of_year(yyyy,mm,day)) ; julian days
yyyyddd=yyyy*1000+ddd
latitude:=Tn&latitude
nlat=dimsizes(latitude)
longitude:=Tn&longitude
nlon=dimsizes(longitude)

;================================================================================================================================================================================================================================= 
;=======================================================================================Heatwaves detection ========================================================================================= 
;================================================================================================================================================================================================================================= 

;==========================================================================================Index 1 =================================================================================== 

;==================Fixed percentiles===========
;;;;;daily anomalies
Tn_daily_clm:=clmDayTLL(tofloat(Tn),yyyyddd)
Tn_daily_smooth:=smthClmDayTLL(Tn_daily_clm,2) ;2 harmonics used meaning to use the annual and semi-annual harmonics only. 
Tn_daily_ano=calcDayAnomTLL(tofloat(Tn),yyyyddd,Tn_daily_smooth)
;percentiles
Tn_mon_90_per=new((/12,nlat,nlon/),"double")
Tn_mon_95_per=new((/12,nlat,nlon/),"double")
Tn_mon_99_per=new((/12,nlat,nlon/),"double")
do i=1,12
do j= 0,nlat-1
do k= 0,nlon-1
mon_ind:=ind(mm .eq. i)
stat=stat_dispersion(Tn(mon_ind,j,k),False)
Tn_mon_90_per(i-1,j,k)=stat(13)
Tn_mon_95_per(i-1,j,k)=stat(25)
Tn_mon_99_per(i-1,j,k)=stat(26)
end do
end do
end do
Tn_mon_per=new((/3,12,nlat,nlon/),"double")
Tn_mon_per(0,:,:,:)=Tn_mon_90_per
Tn_mon_per(1,:,:,:)=Tn_mon_95_per
Tn_mon_per(2,:,:,:)=Tn_mon_99_per

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Heatwave days
;;;;;;;;;;;;;;;;;;;;;;;;;;;;04 day duration
;;;;;;;Tn

Tn_hw_90_per_4d=new((/dimsizes(Tn)/),typeof(Tn))
do j=0,dimsizes(Tn(0,:,0))-1
do k=0,dimsizes(Tn(0,0,:))-1
wakat:=new((/1/),"integer")
do i=0,dimsizes(Tn(:,0,0))-1
if (Tn(i,j,k) .gt. Tn_mon_per(0,mm(i)-1,j,k)) then
wakat:=array_append_record (wakat,i,0)
end if
end do
wakat:=wakat(1:dimsizes(wakat)-1)
wakato:=new((/dimsizes(wakat)/),"integer")
do m=1,dimsizes(wakat)-3
if(wakat(m+2)-wakat(m-1) .eq. 2)then
wakato(m-1:m+2)=wakat(m-1:m+2)
end if
end do
mmm:=wakato(ind(.not.ismissing(wakato)))
Tn_hw_90_per_4d(mmm,j,k)=Tn(mmm,j,k)
end do
end do