;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.

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

begin
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
  data1     = "/home/wrf/ARW/corrida_thomp/"
  data2     = "/home/wrf/ARW/corrida_lin/"
  data3     = "/home/wrf/ARW/corrida_wsm3/"
  files1    = systemfunc("ls -1 "+data1+"wrfout_d01*.nc")
  files2    = systemfunc("ls -1 "+data2+"wrfout_d01*.nc")
  files3    = systemfunc("ls -1 "+data3+"wrfout_d01*.nc")
  numfiles1 = dimsizes(files1)
  numfiles2 = dimsizes(files2)
  numfiles3 = dimsizes(files3)
  print ("Numero de Archivos Dom 1 Thompson = " + numfiles1)
  print ("Numero de Archivos Dom 1 Lin = " + numfiles2)
  print ("Numero de Archivos Dom 1 WSM3 = " + numfiles3)
  a1 = addfiles(files1,"r")
  a2 = addfiles(files2,"r")
  a3 = addfiles(files3,"r")

; We generate plots, but what kind do we prefer?
; type = "x11"
; type = "pdf"
  type = "png"
; type = "ps"
; type = "ncgm"

; Set some basic resources
  res = True
  pltres = True
  pltres@NoTitles = True
  mpres = True
  mpres@mpGeophysicalLineColor = "Black"
  mpres@mpGeophysicalLineThicknessF = "2.0"
  mpres@mpOutlineBoundarySets  = "National"
  mpres@mpNationalLineColor    = "Black"
  mpres@mpNationalLineThicknessF = 2.0
  mpres@mpUSStateLineColor     = "Black"
  mpres@mpGridLineColor        = "Black"
  mpres@mpLimbLineColor        = "Black"
  mpres@mpPerimLineColor       = "Black"
  
  
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; What times and how many time steps are in the data set?
  FirstTime = True
  times1 = wrf_user_getvar(a1[0],"times",-1)  ; get all times in the thompson files
  times2 = wrf_user_getvar(a2[0],"times",-1)  ; get all times in the Lin file
  times3 = wrf_user_getvar(a3[0],"times",-1)  ; get all times in the WSM3 file
  ntimes1 = dimsizes(times1)         ; number of times in the file
  ntimes2 = dimsizes(times2)
  ntimes3 = dimsizes(times3)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 if ((ntimes1 .and. ntimes2 .and. ntimes3) .eq. 72)
 ntimes = 72
 do it = 0,ntimes-1,1             ; TIME LOOP
   times = times1                 ; se asume que todos los tiempos son iguales y se toma el set de tiempos de thompson
   wks = gsn_open_wks(type,"ptot1hracc_ens."+times(it))
   print("Working on time: " + times(it) )
   if (FirstTime) then            ; Save some times for tracking tendencies
     times_sav = times(it)
   end if
   res@TimeLabel = times(it)   ; Set Valid time to use on plots
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need

; Get non-convective, convective and total precipitation thomp
; Calculate tendency values
   rainnc_t = wrf_user_getvar(a1[0],"RAINNC",it)
   rainc_t = wrf_user_getvar(a1[0],"RAINC",it)
   rain_tot_t = rainnc_t + rainc_t

; Get non-convective, convective and total precipitation Lin
; Calculate tendency values
   rainnc_l = wrf_user_getvar(a2[0],"RAINNC",it)
   rainc_l = wrf_user_getvar(a2[0],"RAINC",it)
   rain_tot_l = rainnc_l + rainc_l

; Get non-convective, convective and total precipitation wsm3
; Calculate tendency values
   rainnc_w = wrf_user_getvar(a3[0],"RAINNC",it)
   rainc_w = wrf_user_getvar(a3[0],"RAINC",it)
   rain_tot_w = rainnc_w + rainc_w

    if( FirstTime ) then
      if ( it .eq. 0 ) then
        rainnc_save_t = rainnc_t
        rainnc_save_l = rainnc_l
        rainnc_save_w = rainnc_w

        rainc_save_t = rainc_t
        rainc_save_l = rainc_l
        rainc_save_w = rainc_w

        rain_tot_save_t = rain_tot_t
        rain_tot_save_l = rain_tot_l
        rain_tot_save_w = rain_tot_w

      else
        rainnc_save_t = wrf_user_getvar(a1[0],"RAINNC",it-1)
        rainc_save_t = wrf_user_getvar(a1[0],"RAINC",it-1)
        rain_tot_save_t = rainnc_save_t + rainc_save_t
        
        rainnc_save_l = wrf_user_getvar(a2[0],"RAINNC",it-1)
        rainc_save_l = wrf_user_getvar(a2[0],"RAINC",it-1)
        rain_tot_save_l = rainnc_save_l + rainc_save_l
        
        rainnc_save_w = wrf_user_getvar(a3[0],"RAINNC",it-1)
        rainc_save_w = wrf_user_getvar(a3[0],"RAINC",it-1)
        rain_tot_save_w = rainnc_save_w + rainc_save_w
        
        FirstTime = False
        times_sav = times(it-1)
      end if
    end if
    
; Calculo de ensamble    
    
  ;  prep_tot     = rain_tot - rain_tot_save

    prep_tot_ens = ((rain_tot_t + rain_tot_l + rain_tot_w) / 3) - ((rain_tot_save_t + rain_tot_save_t + rain_tot_save_t)/ 3)  

    prep_tot_ens@description = "Precipitacion horaria acumulada Ensamble"
    prep_tot_ens@units       = "mm"

; Bookkeeping, just to allow the tendency at the next time step
    rainnc_save_t = rainnc_t
    rainc_save_t = rainc_t
    rain_tot_save_t = rain_tot_t

    rainnc_save_l = rainnc_l
    rainc_save_l = rainc_l
    rain_tot_save_l = rain_tot_l

    rainnc_save_w = rainnc_w
    rainc_save_w = rainc_w
    rain_tot_save_w = rain_tot_w

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
d  = systemfunc("date +%d")
y  = systemfunc("date +%Y")
m  = systemfunc("date +%b")
da = stringtointeger(d) + 2
if ( m .eq. "Jan") then
    mes = "Enero"
else if ( m .eq. "Feb") then
    mes = "Febrero"
else if ( m .eq. "Mar") then
    mes = "Marzo"
else if ( m .eq. "Apr") then
    mes = "Abril"
else if ( m .eq. "May") then
    mes = "Mayo"
else if ( m .eq. "Jun") then
    mes = "Junio"
else if ( m .eq. "Jul") then
    mes = "Julio"
else if ( m .eq. "Aug") then
    mes = "Agosto"
else if ( m .eq. "Sep") then
    mes = "Septiembre"
else if ( m .eq. "Oct") then
    mes = "Octubre"
else if ( m .eq. "Nov") then
    mes = "Noviembre"
else if ( m .eq. "Dec") then
    mes = "Diciembre"
end if
end if
end if
end if
end if
end if
end if
end if
end if
end if
end if
end if

if( .not. FirstTime ) then     ; We will skip the first time

      ; Plotting options for Precipitation
        opts_r = res
        opts_r@MainTitle = "Modelo WRF - 27km Res - 72hr Pronostico ~C~~Z80~Weather Research Forecast Model ~C~~Z80~Valido del: "+d+" - "+da+" "+mes+"~C~~Z80~Fecha y Hora: " +times(it)+ " (UTC)"
        opts_r@cnLevelSelectionMode = "ExplicitLevels"
        opts_r@cnLevels             = (/ 0.1, 0.25, 0.5, 1, 1.5, 2, 3, \
                                        4, 5, 7.5, 10, 12.5, 15, 20, 25, 30, 40/)
        opts_r@cnFillColors         = (/"White","Green","Green", \
                                        "Green","Green3","Chartreuse4", \
                                        "Green4","ForestGreen","darkgreen", \
                                        "Yellow","gold","Orange","coral", \
                                        "Red","deeppink3","deeppink","magenta","Violet"/)
        opts_r@cnInfoLabelOn        = False
        opts_r@cnConstFLabelOn      = False
        opts_r@cnFillOn             = True
        opts_r@Footer               = False
        opts_r@InitTime             = False
        opts_r@MainTitlePos         = "Right"
        opts_r@ValidTime            = False


      ; contour_tend = wrf_contour(a[0],wks,prep_tot_ens,opts_r) ; total (color)
        delete(opts_r)

       ; MAKE PLOTS

       ; Total Precipitation 1hr
       ; plot = wrf_map_overlays(a[0],wks,contour_tend,pltres,mpres)


    end if    ; END IF FOR SKIPPING FIRST TIME

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

   times_sav = times(it)
   FirstTime = False
  end do        ; END OF TIME LOOP
 end if
end