load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
; Make a list of all files we are interested in
  DATADir = "/WRF/ECMWF_DEMETER/"

  FILES = systemfunc (" ls -1 " + DATADir + "demeter_daily_1965_* ")
  numFILES = dimsizes(FILES)
  print("numFILES = " + numFILES)
  print(FILES)
  print (" ")

a = addfiles(FILES+".nc","r")
;use join function as each file refers to a different year and/or different initial run time i.e. run starts 08(Aug) 11(Nov)
ListSetType(a, "join")
printVarSummary(a)

;decide on file type for output

type = "pdf"

wks = gsn_open_wks(type,"plt_min2m_allmodels_allensembles_all_days")

; Set some basic resources
  res = True
  res@MainTitle                   = "DEMETER_ENSEMBLES"

  pltres = True
  mpres = True

;extract minimum temperature at 2 m for each of the models used in demeter project ukmo, scwf, cnrm
;get data from file (temperature at each time (180 time steps each is 1 day) for each model (ukmo, scwf, cnrm) for each ensemble (9 for each model)time for each location
;data ukmo, scwf and cnrm are short not float so format needs to be changed
ukmo = short2flt(a[:]->ukmo)
scwf = short2flt(a[:]->scwf)
cnrm = short2flt(a[:]->cnrm)

times = a[0]->step
ensemble = a[0]->number
lat = a[0]->latitude
lon = a[0]->longitude

printVarSummary(ukmo)
printVarSummary(scwf)
printVarSummary(cnrm)

print(times)
print(ensemble)

printVarSummary(lat)
printVarSummary(lon)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; What times and how many time steps are in the data set?
  
  ntimes = dimsizes(times)         ; number of times in the file

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;First look at data from model ukmo

   do iensemble = 0,8,1	; ENSEMBLE LOOP MODELS 0,1...7,8
   do itime = 24,4320,24	; TIME LOOP DAYS 1,2....178,180 FOR EACH RUN THIS IS CHANGED TO UNITS IN FILE IE INTO HOURS 24,4320,24

    

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need        
       t2_uk = ukmo     ; T2 in Kelvin
 
       tc2_uk = t2_uk-273.16                  ; T2 in C

  

;Counts the number of consecutive days that the temperature falls to or below the stipulated temperature.
; specify cut off temperature 'tcrit'
;I am trying to tally the number of cold spells per set of data (i.e. for each ensemble for each model for each year/start of ensemble run) but I need to be able to count how many 1,2,3....35 day duration events occur

tcrit = 0

;create 0/1 array

cold_uk = where(tc2_uk.le.tcrit, 1, 0) ;create 0/1 array
;hot_uk = where(tc2_uk.gt.tcrit, 1, 0)

cold_uk0 = dim_numrun_n(cold_uk, 0, 0) ; opt = 0 means total
cold_uk1 = dim_numrun_n(cold_uk, 1, 0) ; opt = 1 means unique

;optional ... add meta data

cold_uk0@long_name ="total cold runs"
cold_uk0@cold_crit_value = tcrit
copy_VarCoords(tc2_uk, cold_uk0)

cold_uk1@long_name ="unique cold runs"
cold_uk1@cold_crit_values = tcrit
copy_VarCoords(tc2_uk, cold_uk1)


 end do	;end of time loop
 end do	;end of ensemble loop

;Next look at data from model scwf

   do iensemble = 0,8,1	; ENSEMBLE LOOP MODELS 0,1...7,8
   do itime = 24,4320,24	; TIME LOOP DAYS 1,2....178,180 FOR EACH RUN THIS IS CHANGED TO UNITS IN FILE IE INTO HOURS 24,4320,24

    

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need        
                
	t2_sc = scwf  			     ; T2 in Kelvin
 
       tc2_sc = t2_sc-273.16                  ; T2 in C

         

;Counts the number of consecutive days that the temperature falls to or below the stipulated temperature.
; specify cut off temperature 'tcrit'

tcrit = 0

;create 0/1 array

cold_sc = where(tc2_sc.le.tcrit, 1, 0)

cold_sc0 = dim_numrun_n(cold, 0, 0)
cold_sc1 = dim_numrun_n(cold, 1, 0)

;optional ... add meta data

cold_sc0@long_name ="total cold runs"
cold_sc0@cold_crit_value = tcrit
copy_VarCoords(tc2_sc, coldrun)

cold_sc1@long_name ="unique cold runs"
cold_sc1@cold_crit_values = tcrit
copy_VarCoords(tc2_sc, hotrun)


 end do	;end of time loop
 end do	;end of ensemble loop

;First look at data from model ukmo

   do iensemble = 0,8,1	; ENSEMBLE LOOP MODELS 0,1...7,8
   do itime = 24,4320,24	; TIME LOOP DAYS 1,2....178,180 FOR EACH RUN THIS IS CHANGED TO UNITS IN FILE IE INTO HOURS 24,4320,24

    

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Lastly get the variables we will need from cnrm model       
               
	t2_cn = cnrm     			; T2 in Kelvin
 
       tc2_cn = t2_cn-273.16                  ; T2 in C
    
 

;Counts the number of consecutive days that the temperature falls to or below the stipulated temperature.
; specify cut off temperature 'tcrit'

tcrit = 0

;create 0/1 array

cold = where(tc2_cn.le.tcrit, 1, 0)

cold_cn0 = dim_numrun_n(cold, 0, 0)
cold_cn1 = dim_numrun_n(cold, 1, 0)

;optional ... add meta data

cold_cn0@long_name ="total cold runs"
cold_cn0@cold_crit_value = tcrit
copy_VarCoords(tc2_cn, coldrun)

cold_cn1@long_name ="unique cold runs"
cold_cn1@cold_crit_values = tcrit
copy_VarCoords(tc2_cn, hotrun)


 end do	;end of time loop
 end do	;end of ensemble loop

; I now need to add the data from the 3 model runs

tc2_tally = cold_uk + cold_sc +cold_cn
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;I now want to produce plots showning the average number of times in a given period there has been a 1 day duration cold event.... max day duration cold event
;so if the max duration of the cold event is 20 days, I will have 20 plots. I have 3 models each with 9 ensembles so for this period tally/(3*9)

tc2_tally@description = "Average Number of occurances of Consecutive Number of days that Surface Temperature at 2m <= 0 degC"
tc2_tally@units = "duration-days"

    ; Plotting options for T                
      opts = res                         
      opts@cnFillOn = True  
      opts@ContourParameters = (/ 0., 50., 1./)
      opts@gsnSpreadColorEnd = -3  ; End third from the last color in color map
      contour_tc2_tally = wrf_contour(a,wks,coldrun,opts)
      delete(opts)

  

    ; MAKE PLOTS                                       
      plot = wrf_map_overlays(a,wks,(/contour_tc2_tally/),pltres,mpres)
     


end