;**********************************************************************
;2019-03-04  Cressman_Interp.ncl. 
;the wind components extracted from METAR are interpoted  using the function
; obj_anal_ic, over a domain which covers partially west Caribe, central and northern
; South America. as first guess, the wind at 10m of the GFS analisis at 2019-02-26:18z,
;is used.this wind field is shown in Fig.1.
;***********************************************************************
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
delim=" "
;***********************************************************************
; read the file containing de GFS wind 10m over ground at 20190302:18Z *
;***********************************************************************
;**************************************************
DirAnData="/usr/local/Devel/Metar_Macro/SendEmail/"
g = addfile (DirAnData+"wind10m.nc","r")
u10=g->u10
v10=g->v10
lat_g =g->lat_g
lon_g =g->lon_g

;invert the latitude order
lati=lat_g(::-1)
;convert to west latitude
loni=lon_g-360.0
lati@_FillValue    = 999.99
loni@_FillValue    = 999.99
   nlati =dimsizes(lati)
   nloni = dimsizes(loni)
;**********************************************************************
;read the file containing the wind components from METAR 20190302:18Z *
;without fill values in order to plot the barbs                       *
;**********************************************************************
DirVie="/usr/local/Devel/Metar_Macro/SendEmail/"
NamVie="VieMetar1H.txt"
Ecf1 = asciiread(DirVie+NamVie,-1,"string")
latbarb=stringtofloat(str_get_field(Ecf1, 1, delim))
lonbarb=stringtofloat(str_get_field(Ecf1, 2, delim))
latbarb@_FillValue=-99.99
lonbarb@_FillValue=-99.99
U_barb=stringtofloat(str_get_field(Ecf1, 3, delim))
V_barb=stringtofloat(str_get_field(Ecf1, 4, delim))

;**************************************************************
;read the file containing the wind extracted from METAR and   *
;and its coordinates                                          *
;************************************************************ *
DirWea="/usr/local/Devel/Metar_Macro/SendEmail/"
NameWea="VarMetar1H.txt"
Ecf = asciiread(DirWea+NameWea,-1,"string")
Wea= str_get_field(Ecf, 5, delim)
lat1d=stringtofloat(str_get_field(Ecf, 1, delim))
lon1d=stringtofloat(str_get_field(Ecf, 2, delim))
npts = dimsizes(lat1d)
lat1d@_FillValue=-999.99
lon1d@_FillValue=-999.99
Wea@_FillValue=-9999
U_m=new(npts,float)
V_m=new(npts,float)
U_m=stringtofloat(str_get_field(Ecf, 6, delim))
V_m=stringtofloat(str_get_field(Ecf, 7, delim))
U_m@_FillValue=-99.99
V_m@_FillValue=-99.99
  nloni = dimsizes(loni)
  nlati = dimsizes( lati)

;*************************************************
; declare some arrays for Cressman interpolation *
;*************************************************
   temp_u=new((/npts/),"float")
   temp_v=new((/npts/),"float")
   temp_u(:)=U_m(:)
   temp_v(:)=V_m(:)
;**************************************************
  temp_u@_FillValue    = 999.99
  temp_v@_FillValue    = 999.99
;============================================================== ;
  grid_u =new((/nlati,nloni/),typeof(temp_u))
  grid_v =new((/nlati,nloni/),typeof(temp_v))
;**********************************
;start the Cressman interpolation *
;**********************************
    zVal = temp_u(:)
;some options for rscan.
; rscan = (/10, 5, 2/)
 rscan = (/0.1,0.01/)
;    rscan = (/0.1,0.08,0.05/)
; rscan = (/0.05,0.01/)
; rscan = (/0.2,0.1,0.01/)
; rscan = (/0.25,0.2,0.1,0.01/)
    opt = True 
    opt@guess = u10(:,:)  ; first guess = GFS wind at 10m
;    opt@zonal = True       ; first guess = zonal mean
    grid_u(:,:) = obj_anal_ic(lon1d,lat1d,zVal,loni,lati, rscan, opt)  
  grid_u@_FillValue    = -999.99
  grid_u!0="lati"          ;set coordinate variables for the wind to be ploted
  grid_u!1="loni"
  grid_u&lati=lati
  grid_u&loni=loni
  grid_u&loni@units    = "degrees_east"
  grid_u&lati@units     = "degrees_north"
;*****************
;Poisson filling *
;*****************
  guess     = 1                ; use zonal means
  is_cyclic = False             ; cyclic [global]
  nscan     = 1500             ; usually much less than this
  eps       = 1.e-2            ; variable dependent
  relc      = 0.6              ; relaxation coefficient
  opt2       = 0                ; not used
   poisson_grid_fill(grid_u(:,:), is_cyclic, guess, nscan, eps, relc, opt2)

;***********************************************
;repeat the process for the V_m wind component *
;***********************************************
    zVal = temp_v(:)
; some options for rscan.
; rscan = (/10, 5, 2/)
 rscan = (/0.1,0.01/)
;    rscan = (/0.1,0.08,0.05/)
; rscan = (/0.05,0.01/)
; rscan = (/0.2,0.1,0.01/)
; rscan = (/0.25,0.2,0.1,0.01/)
    opt = True
    opt@guess = v10(:,:)  ; Viento hora anterior = "first guess"
;    opt@zonal = True ; esta opción (recomendada) permite que FG = medias zonales
    grid_v(:,:) = obj_anal_ic(lon1d,lat1d,zVal,loni,lati, rscan, opt) 
  grid_v@_FillValue    = -999.99
  grid_v!0="lati"
  grid_v!1="loni"
  grid_v&lati=lati
  grid_v&loni=loni
  grid_v&loni@units    = "degrees_east"
  grid_v&lati@units     = "degrees_north"
;*****************
;Poisson filling *
;*****************
  guess     = 1                ; use zonal means
  is_cyclic = False             ; cyclic [global]
  nscan     = 1500             ; usually much less than this
  eps       = 1.e-2            ; variable dependent
  relc      = 0.6              ; relaxation coefficient
  opt2       = 0                ; not used
   poisson_grid_fill(grid_v(:,:), is_cyclic, guess, nscan, eps, relc, opt2)

;..............................................
;**************
; create plot *
;**************
   pltType = "png"                             ; "ps", "eps", "pdf", "png"
 wks = gsn_open_wks(pltType, "Cressman_Interp") 
   gsn_define_colormap(wks,"nice_gfdl"); choose a color map 
   gsn_reverse_colormap(wks)           ; Reverse the color map.
;set up some common resources, in order to Turn off draw and 
;frame for all plots, because we are going to do overlays
  res           = True
  res@gsnDraw   = False
  res@gsnFrame  = False
  res@gsnMaximize      = True          ; maximize in frame

;resources for the base map   
          mp_res = res
   mp_res@mpDataSetName                    = "Earth..4"
   mp_res@mpDataBaseVersion           = "MediumRes"
   mp_res@mpGeophysicalLineThicknessF = 2
   mp_res@mpNationalLineThicknessF    = 2

;set domain values
   mp_res@mpMinLatF           = 0.5
   mp_res@mpMaxLatF           = 23.5
   mp_res@mpMinLonF           = -93.0
   mp_res@mpMaxLonF           = -60.0

   mp_res@mpFillOn            = False
   mp_res@mpOutlineBoundarySets = "National"    ; turn on country boundaries
   mp_res@trGridType          = "TriangularMesh"   ; Necessary b/c lat, lon
                                                ; arrays have missing values.
   mp_res@tiMainString        = "Cressman_Interp"

; resources fro wind plot
        res_vc                            = res
        res_vc@vcGlyphStyle               = "LineArrow"
        res_vc@vcLineArrowThicknessF      = 3
         res_vc@vcMinDistanceF             = 0.02
        res_vc@vcRefLengthF               = 0.03
        res_vc@vcRefMagnitudeF           = 10 

        res_vc@vcRefAnnoString1          = "10"
        res_vc@vcRefAnnoSide             = "Top"
        res_vc@vcRefAnnoString2On        = False
        res_vc@vcRefAnnoPerimOn          = False
        res_vc@vcRefAnnoOrthogonalPosF   = -0.12
        res_vc@vcRefAnnoParallelPosF     = 0.999
        res_vc@vcRefAnnoBackgroundColor  = "Purple"
  plot = gsn_csm_map(wks,mp_res)
  vector  = gsn_csm_vector(wks,grid_u,grid_v,res_vc)
 overlay(plot,vector)		; do overlay 

;********************************
;plot state boundaries
;********************************
dp = addfile("/usr/local/Devel/Metar_Macro/SendEmail/COL_adm1.shp", "r") ; Open shapefile
;
; Read data off shapefile
;
  segments = dp->segments
  geometry = dp->geometry
  segsDims = dimsizes(segments)
  geomDims = dimsizes(geometry)
;
; Read global attributes  
;
  geom_segIndex = dp@geom_segIndex
  geom_numSegs  = dp@geom_numSegs
  segs_xyzIndex = dp@segs_xyzIndex
  segs_numPnts  = dp@segs_numPnts
  numFeatures   = geomDims(0)
;*************************************************
; Section to add polylines to map.
;*************************************************
  lines = new(segsDims(0),graphic)   ; array to hold polylines

  plres             = True           ; resources for polylines
  plres@gsLineDashPattern = 2
  plres@gsLineColor = "blue"
  londp    = dp->x
  latdp    = dp->y
  segNum = 0       ; Counter for adding polylines
  do i=0, numFeatures-1  
     startSegment = geometry(i, geom_segIndex)
     numSegments  = geometry(i, geom_numSegs)
     do seg=startSegment, startSegment+numSegments-1
        startPT = segments(seg, segs_xyzIndex)
        endPT   = startPT + segments(seg, segs_numPnts) - 1
        lines(segNum) = gsn_add_polyline(wks,plot,londp(startPT:endPT),  \
                                                    latdp(startPT:endPT), plres)

        segNum = segNum + 1
     end do
  end do

draw(plot)
;****************
;plot the barbs *
;****************
   wmsetp("wdf", 1) 	; barb oriented as adopted in meteorology
   wmsetp("wbs",0.025)  ;The size, of the wind barb shaft. 
   wmsetp("blw",4.0)  ;line with of the barb:
   wmsetp("col", 2)    ; Draw in red (doesn't works!! why?)
   size = wmgetp("wdf")
   wmbarbmap(wks, latbarb, lonbarb, U_barb, V_barb)  ; Plot barbs.
frame(wks)   ; Advanced frame.
end