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)