;----------------------------------------------------------
; This function calculates a contamination rose from two
; time series concentration and wind direction it also
; The limits are 16 groups between 0 and 360
;----------------------------------------------------------

undef("contamination_rose")
function contamination_rose(serie0[*]:float, serie1[*]:float, serie2[*]:float)
;serie0: concentration of compound
;serie1: wind direction
;serie2: wind intensity
local data, ndata0, ndata1, limits
limits      = (/11.25, 33.75, 56.25, 78.75,101.25,123.75,146.25,168.75,\
191.25,213.75,236.25,258.75,281.25,303.75,326.25,348.75/)

ndata0      = dimsizes(serie0)
ndata1      = dimsizes(serie1)
ndata2      = dimsizes(serie2)
stmt0       = .not.(ndata0 .eq. ndata1)
stmt1       = .not.(ndata0 .eq. ndata2)
stmt2       = .not.(ndata1 .eq. ndata2)
stmt        = stmt0 .or. stmt1 .or. stmt2
if (stmt) then
   print("ERROR, your data must have the sime lenght")
   print("Returning an empty data")
   data = new((/1/),"float")
else
   data        = new((/3,17/),"float")
   ; 0: Accumulated Amount 1: N values 2: Accumated wind intensity 
   data(0:,0:) = 0.0
   do it = 0, ndata0-1
      if  (.not. ismissing(serie0(it)) )  then
         if ((serie1(it) .gt. limites(15)) .or. (serie1(it) .le. limites(0))) then
            if (serie2(it) .ne. 0.) .xor. (serie2(it) .gt. 360.) then
               data(0,0) = data(0,0) + serie0(it)
               data(1,0) = data(1,0) + serie2(it)
               data(2,0) = data(2,0) + 1
            end if
         end if
         if ((serie1(it) .gt. limites(0)) .and. (serie1(it) .le. limites(1))) then
            data(0,1) = data(0,1) + serie0(it)
            data(1,1) = data(1,1) + serie2(it)
            data(2,1) = data(2,1) + 1
         end if
         if ((serie1(it) .gt. limites(1)) .and. (serie1(it) .le. limites(2))) then
            data(0,2) = data(0,2) + serie0(it)
            data(1,2) = data(1,2) + serie2(it)
            data(2,2) = data(2,2) + 1
         end if
         if ((serie1(it) .gt. limites(2)) .and. (serie1(it) .le. limites(3))) then
            data(0,3) = data(0,3) + serie0(it)
            data(1,3) = data(1,3) + serie1(it)
            data(2,3) = data(2,3) + 1
         end if
         if ((serie1(it) .gt. limites(3)) .and. (serie1(it) .le. limites(4))) then
            data(0,4) = data(0,4) + serie0(it)
            data(1,4) = data(1,4) + serie2(it)
            data(2,4) = data(2,4) + 1
         end if
         if ((serie1(it) .gt. limites(4)) .and. (serie1(it) .le. limites(5))) then
            data(0,5) = data(0,5) + serie0(it)
            data(1,5) = data(1,5) + serie2(it)
            data(2,5) = data(2,5) + 1
         end if
         if ((serie1(it) .gt. limites(5)) .and. (serie1(it) .le. limites(6))) then
            data(0,6) = data(0,6) + serie0(it)
            data(1,6) = data(1,6) + serie2(it)
            data(2,6) = data(2,6) + 1
         end if
         if ((serie1(it) .gt. limites(6)) .and. (serie1(it) .le. limites(7))) then
            data(0,7) = data(0,7) + serie0(it)
            data(1,7) = data(1,7) + serie2(it)
            data(2,7) = data(2,7) + 1
         end if
         if ((serie1(it) .gt. limites(7)) .and. (serie1(it) .le. limites(8))) then
            data(0,8) = data(0,8) + serie0(it)
            data(1,8) = data(1,8) + serie2(it)
            data(2,8) = data(2,8) + 1
         end if
         if ((serie1(it) .gt. limites(8)) .and. (serie1(it) .le. limites(9))) then
            data(0,9) = data(0,9) + serie0(it)
            data(1,9) = data(1,9) + serie2(it)
            data(2,9) = data(2,9) + 1
         end if
         if ((serie1(it) .gt. limites(9)) .and. (serie1(it) .le. limites(10))) then
            data(0,10) = data(0,10) + serie0(it)
            data(1,10) = data(1,10) + serie2(it)
            data(2,10) = data(2,10) + 1
         end if
         if ((serie1(it) .gt. limites(10)) .and. (serie1(it) .le. limites(11))) then
            data(0,11) = data(0,11) + serie0(it)
            data(1,11) = data(1,11) + serie2(it)
            data(2,11) = data(2,11) + 1
         end if
         if ((serie1(it) .gt. limites(11)) .and. (serie1(it) .le. limites(12))) then
            data(0,12) = data(0,12) + serie0(it)
            data(1,12) = data(1,12) + serie2(it)
            data(2,12) = data(2,12) + 1
         end if
         if ((serie1(it) .gt. limites(12)) .and. (serie1(it) .le. limites(13))) then
            data(0,13) = data(0,13) + serie0(it)
            data(1,13) = data(1,13) + serie2(it)
            data(2,13) = data(2,13) + 1
         end if
         if ((serie1(it) .gt. limites(13)) .and. (serie1(it) .le. limites(14))) then
            data(0,14) = data(0,14) + serie0(it)
            data(1,14) = data(1,14) + serie2(it)
            data(2,14) = data(2,14) + 1
         end if
         if ((serie1(it) .gt. limites(14)) .and. (serie1(it) .le. limites(15))) then
            data(0,15) = data(0,15) + serie0(it)
            data(1,15) = data(1,15) + serie2(it)
            data(2,15) = data(2,15) + 1
         end if
      end if
   end do
   data(0,16) = data(0,0)
   data(1,16) = data(1,0)
   data(2,16) = data(2,0)
end if
return (data)
end

;------------------------------------------------------------------------------
; Summary of statistics to compare modeled timeseries and observed
; Timeseries. Both arrays need to be 1D and have the same length
;------------------------------------------------------------------------------
undef("calc_statistics")
function calc_statistics ( xm[*]:float, xo[*]:float  )
; xm Serie modelada 
; xo Serie Observada
local nobs, nmod, x_stats, fraction, condition, n_data

begin
nobs  = dimsizes(xo)
nmod  = dimsizes(xm)
if (.not.(nobs .eq. nmod)) then
   print("ERROR, your data must have the sime lenght")
   print("Returning an empty data")
   x_stats = new((/1/),"float")
else
   x_stats = new((/10/),"float")
   n_data       = nobs

   ; Correlacion de pearson (R)
   x_stats(0)   = escorc(xm,xo)

   ; Mean Error (ME)
   x_stats(1)   = dim_avg_Wrap(xm - xo)

   ; Mean Absolute Error (MAE)
   x_stats(2)   = dim_avg_Wrap(abs( xm - xo ))

   ; Root Mean Square Error (RMSE)
   x_stats(3)   = dim_rmsd(xm,xo)

   ; Normalized Mean Square Error (NMSE)
   x_stats(4)   = dim_rmsd(xm,xo)/( dim_avg(xm) * dim_avg(xo) )

   ; Index of Agreement (IOA)  Wilmott et al., 1985 with absolute value in numerador
   x_stats(5)   = 1. - dim_sum( abs(xm - xo))/( dim_sum( abs(xm - dim_avg_Wrap(xo) ))+\
   dim_sum( abs( xm - dim_avg_Wrap(xo) )) )

   ; Fractional Bias (FB)
   x_stats(6)   = (dim_avg_Wrap(xo) - dim_avg_Wrap(xm))/\
   0.5*(dim_avg_Wrap(xo) + dim_avg_Wrap(xm))

   ; Geometric Mean Bias (MG)
   x_stats(7)   = exp( log(dim_avg_Wrap(xo)) - log(dim_avg_Wrap(xm)) )

   ; Geometric Variance (VG)
   x_stats(8)   = exp( dim_rmsd( log(xm),log(xo) ) )
   ; Factor of Two Observations (FAC2)

   if any(xo .eq. 0.) then
      fraction  = xm/where(xo .ne. 0., xo,xo@_FillValue)
   else
      fraction  = xm/xo
   end if
   condition    = num((fraction .ge. 0.5) .and. (fraction .le. 2.0))
   x_stats(9)   = condition/tofloat(n_data)
end if

return(x_stats)

end