;***********************************************
; roms_4.ncl
;
; Concepts illustrated:
;   - Plotting ROMS data
;   - Drawing curly vectors
;   - Loading NCL functions from another script
;***********************************************
; Example of using a donated library: ROMS_utils.ncl
;
; Specifically: roms depth slice using roms_3d_interp
;***********************************************
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/ROMS_utils.ncl"


begin
                     
   ; mes = (/ "Jan", "Fev", "Mar", "Abr", "Mai", "Jun", "Jul", "Ago", "Set", "Out", "Nov", "Dez" /)
   ; ano = (/ "1980", "1981", "1982", "1983", "1984", "1985" /)

   mes   = (/ "Dez" /)
   ; mes_n = (/"01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12" /)
   mes_n = (/"12" /)
   ano   = (/ "1981" /)

   do y=0,dimsizes(ano)-1
    do m=0,dimsizes(mes)-1

      ;***********************************************
      ; User settings
      ;***********************************************
         date     = ano(y) + mes_n(m) + "16"
         path1    ="/home/leilane/RESULTADOS/006nested/"
         
         fhis1    = path1 + "nested_avg_" + date + "_LEILANE.nc"
         
         outfile = "roms"
         variable= "temp"
         depth   = -1
         rec     = 0   
      ;***********************************************
      ; Read file date and use ROMS_utils.ncl
      ;***********************************************
         his1   =  addfile (fhis1,"r")
         
         grid  = addfile ("/home/leilane/ROMS/GRD_files/mesh_donor_grd_004nested_LEILANE.nc","r")
         lon2d = grid->lon_rho
         lat2d = grid->lat_rho
         
         out1   = roms_3d_interp(his1,grid,variable,rec,depth)
         

         print("Done with roms_3d_interp")
         out1@lat2d = lat2d
         out1@lon2d = lon2d

         ; .....out1 is calculated as before....

         out1Missing = ismissing(out1)
         out1MissingInt = where(out1Missing.eq.True, 1, 0)
         print(dimsizes(out1MissingInt))
         ; wks=gsn_open_wks("x11","foo")
         ; plot= gsn_csm_contour(wks,out1MissingInt(:,:),False)
         ; exit
         ;print(dimsizes(out1))

         ; sst = out1(58,106)
         ; asciiwrite("temp.txt", sst)

         
      ;************************************************
      ; create plot
      ;************************************************
        wks_type = "newpng"        ; or "ps"
        wks_type@wkOrientation = "Portrait"
        wks = gsn_open_wks("png","TSM_" + ano(y) + mes_n(m) + "verifica")           ; send graphics to PNG file
                                                        
        vres1 = True               ; plot mods desired
        vres1@gsnDraw              = False
        vres1@gsnFrame             = False
        vres1@gsnMaximize          = True    ; Maximize plot in frame
        vres1@gsnPaperOrientation  = "Portrait"
        vres1@cnFillDrawOrder      = "Draw"
        vres1@cnFillOn             = True               ; turn on color for contours
        ; vres1@cnConstFEnableFill   = True
        vres1@cnFillPalette        = "BlGrYeOrReVi200"
        vres1@cnLinesOn            = False              ; turn off contour lines
        vres1@cnLineLabelsOn       = False              ; turn off contour line labels
        vres1@cnFillMode           = "AreaFill"
        vres1@cnRasterSmoothingOn  = True
        ; vres1@gsnScalarContour     = True               ; contours desired
        vres1@cnInfoLabelOn        = False            ; turn off cn info label
        vres1@lbLabelBarOn         = True            ; turn off individual cb's
       
        lon1   = min(lon2d)
        lat1   = min(lat2d)
        lon2   = max(lon2d)
        lat2   = max(lat2d)


      ; MAP
        vres1@mpMinLonF                        = lon1             ; longitude oeste
        vres1@mpMaxLonF                        = lon2             ; longitude leste
        vres1@mpMinLatF                        = lat1             ; latitude sul
        vres1@mpMaxLatF                        = lat2             ; latitude norte
        vres1@mpDataBaseVersion      = "LowRes"          ; use high resolution coast
        vres1@pmTickMarkDisplayMode  = "Always"           ; turn on tickmarks
        vres1@mpLandFillColor      = "gray"
         
        vres1@cnLevelSelectionMode =  "AutomaticLevels"   
        ; vres1@cnMinLevelValF       = 18.
        ; vres1@cnMaxLevelValF       = 28.
        ; vres1@cnLevelSpacingF      =   .5 

        vres1@gsnLeftString  = "Simulacao"
        plot = gsn_csm_contour_map(wks,out1,vres1)     

        draw(plot)
        frame(wks)

      end do
  end do 

end