; ***********************************************
; xy_2.ncl
;
; Concepts illustrated:
;   - Drawing an XY plot with multiple curves
;   - Changing the line color for multiple curves in an XY plot
;   - Changing the line thickness for multiple curves in an XY plot
;   - Drawing XY plot curves with both lines and markers
;   - Changing the default markers in an XY plot
;   - Making all curves in an XY plot solid
;
; ***********************************************
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/csm/contributed.ncl"
;************************************************
begin
;************************************************
; read in data
;************************************************
;
   diri = "./"
   fili = systemfunc("cd "+diri+" ; ls rce-bret-2005-rad.*.nc")
   nfili= dimsizes(fili)
   print("nfili="+nfili)
;
   f    = addfiles(diri+fili,"r") 
   ListSetType (f, "join") 
;
;    vNam = getfilevarnames (f)  ; all variables on file   
;    print(vNam)
;
    vNam1 = (/"x","y","z","p","raing","prcate1","tke","radfrc","radflx","ptsflx","qvsflx"/)  ; manually specify
    print(vNam1)
;
;	  
     x   = f[0]->x
     y   = f[0]->y
     z   = f[0]->z 
     p  = f[:]->p
     PREC= f[:]->raing
     RADQR= f[:]->rnflx
     LHF = f[:]->ptsflx
     SHF = f[:]->qvsflx
;
    printVarSummary(x)
    printVarSummary(y)
    printVarSummary(z)           ; [74]
    printVarSummary(p)      
;
; create time 
      TIME3D = yyyymmddhh_time(2017, 2017, 1, "integer")
      printVarSummary(TIME3D)
;      time3D = TIME3D({2017040100:2017051000})    ; coordinate subscripting
      time3D = TIME3D({2017040100:2017040203})
      time3D!0="time"
      time3D&time=time3D
      printVarSummary(time3D)  
;
;
      print(time3D) 
;
      yrfrac = yyyymmddhh_to_yyyyfrac(time3D, 0)
      print(yrfrac)
; 
     dim3D= dimsizes(p)
     print(dim3D)
     ntim3D= dim3D(0)
     klev3D = dim3D(1)
     ny3D   = dim3D(2)
     nx3D   = dim3D(3)
     p!0="time"
     p&time=time3D
     p!1="z"
     p&z=z
     p!2="y"
     p&y=y
     p!3="x"
     p&x=x
     printVarSummary(p)
;
; create sfc pressure ps
     ps= p(z|:,y|:,x|:,time|:)
     copy_VarCoords(p,ps)
     copy_VarAtts(p,ps)
     psfc=ps(0,:,:,:)     
    printVarSummary(psfc) 
     Avepsfc= new((/dimsizes(time3D)/),"float") 
      Avepsfc= dim_avg_n_Wrap( psfc, (/0,1/))
     printVarSummary(Avepsfc)     
;
;
;    print(p)
    printVarSummary(PREC)
    printVarSummary(LHF)
    printVarSummary(SHF)
    printVarSummary(RADQR)
;
;
;
    maxt=max(time3D)
    mint=min(time3D)
    print(maxt)
    print(mint)
;
      AveRADQR= dim_avg_n_Wrap( RADQR, (/1,2/))
;
    copy_VarAtts (RADQR, AveRADQR)         ; copy variable attributes
    copy_VarCoords_2 (RADQR, AveRADQR)     ; copy left coordinate arrays if present
;     
     printVarSummary(AveRADQR) 
;    print(AveRADQR)
;
  AveQrad  = new((/dimsizes(time3D)/),"float")  
  AveQrad= (AveRADQR*1005.7*Avepsfc)/(9.81*8.64*10000.)
;
  copy_VarAtts (AveRADQR,AveQrad)         ; copy variable attributes
  copy_VarCoords (AveRADQR,AveQrad)     ; copy left coordinate arrays if present
;     
  printVarSummary(AveQrad)
;
;  print(AveQrad)
; define QSEN ans QLAT
;
     cp=1004.0
     lv=2.5*1000000
     THF=LHF*lv+SHF*cp
     copy_VarAtts (LHF, THF)         ; copy variable attributes
     copy_VarCoords (LHF, THF)     ; copy left coordinate arrays if present      
     AveTHF= dim_avg_n_Wrap( THF, (/1,2/))
     printVarSummary(AveTHF) 
     QTOT=AveTHF+AveQrad
     copy_VarAtts (AveQrad,QTOT)         ; copy variable attributes
     copy_VarCoords(AveQrad, QTOT)     ; copy left coordinate arrays if present         
     printVarSummary(QTOT) 
;*********************************************************************************
      AvePREC= dim_avg_n_Wrap( PREC, (/1,2/))
;
     copy_VarAtts (PREC, AvePREC)         ; copy variable attributes
     copy_VarCoords_2  (PREC, AvePREC)      ; copy left coordinate arrays if present         
     printVarSummary(AvePREC) 
;*******************************************************************************   
; 
 data      = new((/4,dimsizes(time3D)/),"float")
;
 data(0,:) = AvePREC(:)
 data(1,:) = AveTHF(:)
 data(2,:) = AveQrad(:)
 data(3,:) = QTOT(:)
;
   printVarSummary(data)
;********************************************************************
  wks = gsn_open_wks("png","ARPS-PROVA")   ; send graphics to PNG file
;	
  plot = new(4,graphic)	       ; Array to hold plots so we can panel them later.
;
  res           = True
  res@trXMinF   = mint       ; starting point along X axis
  res@trXMaxF   = maxt      ; ending point along X-axis
  res@vpWidthF  = 0.7          ; stretch the plot to be wider (in NDC units)
  res@vpHeightF = 0.25         ; and not as tall 
  res@vpXF      = 0.15         ; set the start point along the X-axis in NDC units
  res@gsnDraw   = False        ; do not draw the plot (plots will be paneled)
  res@gsnFrame  = False        ; do not advance the frame (plots will be paneled)
;*************************************************************************
   res@tiMainString = "Surface Precipitation          "
   res@xyLineColors         = (/"black"/) 
   res@tiXAxisString        = "days"  
   res@tiXAxisPosition      = "Center"
   res@tiYAxisString        = "Raing mm"  
   plot(0) = gsn_csm_xy(wks,time3D,data(0,:),res)	

   res@tiMainString = "AveTHF          "
   res@xyLineColors         = (/"purple"/) 
   res@tiXAxisString        = "days"  
   res@tiYAxisString        = " W/m2 "  
   res@tiXAxisPosition      = "Center"

   plot(1) = gsn_csm_xy(wks,time3D,data(1,:),res)

   res@tiMainString = " Net SW flux at sfc     "
    res@xyLineColors         = (/"red"/) 
   res@tiXAxisString        = "days"  
   res@tiXAxisPosition      = "Center"
   res@tiYAxisString        = "SWNS W/m^2"  
   plot(2) = gsn_csm_xy(wks,time3D,data(2,:),res)
;
   res@tiMainString = " LHF+SHF+AveQrad    "
    res@xyLineColors         = (/"blue"/) 
   res@tiXAxisString        = "days"  
   res@tiXAxisPosition      = "Center"
   res@tiYAxisString        = "QRAD W/m^2"  
   plot(3) = gsn_csm_xy(wks,time3D,data(3,:),res)
;
;*************************************************************************
   resP             = True
   resP@gsnMaximize = True
   resP@txString           ="bret-marat-LD_rad_cos_hdf_day_perpetual [40d]"     ; add title
   gsn_panel(wks,plot,(/3,2/),resP)
end