Subversion Repositories lagranto.ocean

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"
;-------------------------------------------------------------------------------------;
;-------------------------------------------------------------------------------------;


;-------------------------------------------------------------------------------------;
; FUNCTION TO READ IN TRAJECTORIES
;-------------------------------------------------------------------------------------;
undef("read_traj")
function read_traj (lslfile:string, flag:string, nshift:integer)
;-------------------------------------------------------------------------------------;
;       Input : Lagranto ECMWF lsl output, old or new Lagranto ?
;       Output: 4d Trajectorie array with Traj(number of traj, timestep, fields)
;-------------------------------------------------------------------------------------;
local nhead, head, delim, data, n, ifirst, ifinal, timestp, ntra, traj, vars, flag, nshift
;-------------------------------------------------------------------------------------;
begin
;---------- Read file ----------------------------------------------------------------;
print(" ")
if (flag .eq. "new") then
print("Assuming Lagranto 2.0 header")
end if
if (flag .eq. "old") then
print("Assuming old Lagranto header")
end if
if (flag .eq. "oldlm") then
print("Assuming old Lagranto LM header")
end if

; Delimiting character
delim   = " "

; Number of header lines
nhead   = 4

; Read header lines
head    = systemfunc("head -n"+nhead+" "+ lslfile)

; Number of columns
ncol    = str_fields_count(head(2),delim)

; Read data in as float
data    = asciiread(lslfile,-1,"float")

; Number of entries in data
n       = dimsizes(data)

; First data element
if (flag .eq. "new") then
ifirst  = 3 + nshift
end if
if (flag .eq. "oldlm") then
ifirst  = 4 + nshift
end if
if (flag .eq. "old") then
ifirst  = 5 + nshift
end if

; Last data element
ifinal  = n-1

; Period
period  = data(ifinal-ncol+1)  - data(ifirst)

; Time step of trajectories
timestp = data(ifirst+ncol) - data(ifirst)

; Number of times
ntim    = floattoint( period/timestp + 1 )

; Number of trajectories
ntra    = floattoint( (n-ifirst) / (ntim * ncol ) )

; 3d trajectory array
traj    = onedtond(data(ifirst:ifinal),(/ntra,ntim,ncol/))

; Field names (columns)
vars    = str_split(head(2),delim)

;-------------- Add meta information 
if (flag .eq. "new") then
; Add some attributes to the variable
traj@HOUR        = stringtoint( str_get_cols(head(0),24,25) )
traj@DAY         = stringtoint( str_get_cols(head(0),21,22) )
traj@MONTH       = stringtoint( str_get_cols(head(0),19,20) )
traj@YEAR        = stringtoint( str_get_cols(head(0),15,18) )
traj@FIELD_NAMES = vars
traj@TIMEPERIOD  = period
traj@TIMESTEP    = timestp
traj@MODEL       = "LAGRANTO 2.0"
end if

if (flag .eq. "old") then
traj@HOUR        = stringtoint( str_get_cols(head(0),22,23) )
traj@DAY         = stringtoint( str_get_cols(head(0),19,20) )
traj@MONTH       = stringtoint( str_get_cols(head(0),17,18) )
traj@YEAR        = stringtoint( str_get_cols(head(0),13,16) )
traj@FIELD_NAMES = vars
traj@FIELD_NAMES = vars
traj@TIMEPERIOD  = period
traj@TIMESTEP    = timestp
traj@MODEL       = "Lagranto"
end if

; Set names of dimensions and assign values
traj!0           = "TRA"
traj!1           = "TIME"
traj!2           = "FIELD"
traj&TIME        = fspan(0,period,ntim)

; print info
printVarSummary(traj)

return(traj)
print(" ")

;---- Clean
delete(traj)
delete(vars)
delete(data)

;---- end
end
;-------------------------------------------------------------------------------------;
;-------------------------------------------------------------------------------------;


;-------------------------------------------------------------------------------------;
; START OF NCL SCRIPT
;-------------------------------------------------------------------------------------;

system("\rm ._*")


; --------- Name of trajectory file --------------------------------------------------;
filename="LSL20040108_00"

; --- read trajectories
traj = read_traj(filename,"new",0)

; -- get infos
dims = dimsizes(traj)
ntra = dims(0)
ntim = dims(1)
fld  = dims(2)

; set field to plot
cnIndex = 3

; set the contour range
cnLev   = fspan(900,1210,31)

; --- create the plot
wks  = gsn_open_wks("eps","test.") 
gsn_define_colormap(wks,"MPL_YlGnBu")

; --- the map
res                        = True  
res@gsnDraw                        = False                         
res@gsnFrame                       = False                        
    
res@mpLandFillColor        = "gray50" 
res@mpDataBaseVersion      = "MediumRes"  
res@mpGreatCircleLinesOn   = True 
 
res@mpMinLonF              = -120
res@mpMaxLonF              =  -50
res@mpMinLatF                      =   10
res@mpMaxLatF              =   80
res@mpCenterLonF           =    0     
   
map             = gsn_csm_map_ce(wks,res)  

; --- plot starting region

gres                 = True
gres@gsLineColor     = "black"
gres@gsLineDashPattern = 0
gres@gsLineThicknessF  = 1.
gres@tfPolyDrawOrder = "PreDraw"
gres@gsMarkerColor = "red"
gres@gsMarkerIndex = 15

la = (/ 30,  35,  35,  30,  30/)
lo = (/-75, -75, -70, -70, -75/)

;dum0 = gsn_add_polyline(wks,map,lo,la,gres)
 dum0 = gsn_add_polymarker(wks,map,-72.5,35,gres)    

; --- the trajectories
pres                                       = True                         
pres@gsLineThicknessF      = 0.25 
                       

do i = 0,ntra-1
if (traj(i,0,cnIndex).gt.1000 .and. traj(i,0,cnIndex).lt.1100  ) then

  x = traj(i,0:(ntim-1),1)    ; Longitudes for trajectory i
  y = traj(i,0:(ntim-1),2)    ; Latitudes for trajectory i
 
 do t = 0,ntim-2

        pres@gsLineColor = get_color_index("MPL_YlGnBu",cnLev,traj(i,t,cnIndex))  
    gsn_polyline  (wks,map,(/x(t),x(t+1)/),(/y(t),y(t+1)/),pres)  
    
end do
end if     

end do  

; --- Draw labelbar
cnLabels                          = flt2string(cnLev)        
cmap                                              = read_colormap_file("MPL_YlGnBu")

csizes = dimsizes(cmap)
nboxes = dimsizes(cnLev)+1                ; # of labelbar boxes
clen   = csizes(0)                    ; # of colors in color map
stride = ((clen-1) - 2) / nboxes      ; Start at color index 2 and end
                                      ; near color index clen-1.
                                      
fill_colors                               = ispan(2,clen-1,stride)

lres                          = True
lres@lbTitleString            = "[m]"           
lres@lbTitleFontHeightF       = 0.015
lres@vpWidthF                 = 0.60                    
lres@vpHeightF                = 0.10                    
lres@lbPerimOn                = False                   
lres@lbOrientation            = "Horizontal"
                
lres@lbFillColors             = fill_colors             
lres@lbMonoFillPattern        = True    
lres@lbLabelFontHeightF       = 0.015                   
lres@lbLabelAutoStride        = True                   

lbid = gsn_create_labelbar(wks,nboxes,cnLabels,lres)
                                      
amres = True
amres@amZone           = 2
amres@amParallelPosF   = 0.5            ; Center labelbar.
amres@amOrthogonalPosF = 0.1            ; Move down, away from plot

annoid = gsn_add_annotation(map,lbid,amres)



; ---- create the plot
draw(map)

frame(wks)