<html>
  <head>

    <meta http-equiv="content-type" content="text/html; charset=ISO-8859-1">
  </head>
  <body text="#000000" bgcolor="#FFFFFF">
    Dear folks,
    <br>
    <br>
    Although I have carefully read around the corresponding pages on the
    topic (.../2Dvertcoord.shtml), I am getting into troubles for
    plotting a vertical contour section of the 3D variable w from ROMS
    model and overlay the model topography.
    <br>
    <br>
    1) It looks that the topography doesn't match the last sigma layer
    and I would bet on an error overlaying the topography. (Figures
    with_topo_overlaid_mesh.ps and with_topo_overlaid_nomesh.ps with
    respect to without_topo_overlaid.ps)
    <br>
    <br>
    <br>
    2) Another error is that I "loose" my xticks labels when switching
    trGridType from triangularMesh (to turn the edges on as in
    2dvertcoords_2) to curvilinear grid (to get a nicer plot) with that
    error message:
    <br>
    warning:tmXBLabelFontHeightF is not a valid resource in
    w.ps_contour.PlotManager at this time
    <br>
    warning:tmXBMode is not a valid resource in w.ps_contour.PlotManager
    at this time
    <br>
    warning:tmXBValues is not a valid resource in
    w.ps_contour.PlotManager at this time
    <br>
    warning:tmXBLabels is not a valid resource in
    w.ps_contour.PlotManager at this time
    <br>
    <br>
    (Figures with_topo_overlaid_mesh.ps and with_topo_overlaid_nomesh.ps
    again)
    <br>
    <br>
    <br>
    3) Last error of mine, there is no frame to my contour plot when no
    grid nor tickmarks when no topography overlaid. Do see figure
    without_topo_overlaid_nomesh.ps
    <br>
    <br>
    <br>
    Hereafter you will find the full script. Many thanks in advance for
    any help
    <br>
    <br>
    Best regards,
    <br>
    Laurent
    <br>
    <br>
    <br>
    <br>
    <br>
    <br>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    <br>
    ; ;
    <br>
    ;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Copyright (C)&nbsp;
    2016&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ;
    <br>
    ;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Centre National de la Recherche Scientifique
    CNRS&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ;
    <br>
    ; ;
    <br>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    <br>
    ;
    <br>
    ;&nbsp;&nbsp; File:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; croco_vslice.ncl
    <br>
    ;&nbsp;&nbsp; Author:&nbsp;&nbsp;&nbsp;&nbsp; Laurent ROBLOU
    <br>
    ;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Centre National de la Recherche Scientifique CNRS
    <br>
    ;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Laboratoire d Aerologie, Toulouse
    <br>
    <br>
    ; this program interpolates the vertical section of a 3D variable
    along <br>
    ; the transect line approximated by a great circle section and
    overlay<br>
    ; the underlying sea floor topography<br>
    ; interpolation is made at each time step for the specified 3D var<br>
    <br>
    <br>
    ; INIT LIBS AND PLOT
    <br>
    ; ========================
    <br>
    load "/opt/ncl/ncl-6.1.2-LA/lib/ncarg/nclscripts/csm/gsn_code.ncl"
    <br>
    load "/opt/ncl/ncl-6.1.2-LA/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
    <br>
    load
    "/opt/ncl/ncl-6.1.2-LA/lib/ncarg/nclscripts/csm/contributed.ncl"
    <br>
    load
"/opt/ncl/ncl-6.1.2-LA/lib/ncarg/nclscripts/contrib/calendar_decode2.ncl"<br>
    <br>
    <br>
    <br>
    begin
    <br>
    <br>
    <br>
    ;************************************
    <br>
    ; get model variables
    <br>
    ;************************************
    <br>
    <br>
    f = addfile("test_noxios.nc","r")
    <br>
    <br>
    var&nbsp; = f-&gt;$"w"$
    <br>
    if (var@units .eq. 1) then
    <br>
    &nbsp;&nbsp;&nbsp; var@units="-"
    <br>
    end if
    <br>
    <br>
    if(isatt(var,"long_name"))
    <br>
    &nbsp;&nbsp;&nbsp; varname=var@long_name
    <br>
    else
    <br>
    &nbsp;&nbsp;&nbsp; if(isatt(var,"standard_name"))
    <br>
    &nbsp;&nbsp;&nbsp; varname=var@standard_name
    <br>
    &nbsp;&nbsp;&nbsp; else
    <br>
    &nbsp;&nbsp;&nbsp; varname=var@short_name ; last tentative
    <br>
    &nbsp;&nbsp;&nbsp; end if
    <br>
    end if
    <br>
    <br>
    lat = f-&gt;lat_rho&nbsp; ; get lat
    <br>
    lon = f-&gt;lon_rho&nbsp; ; get lon
    <br>
    time = f-&gt;time(0:2)&nbsp; ; get time frame
    <br>
    <br>
    vlev=f-&gt;s_rho
    <br>
    f = addfile ("test_noxios.nc", "r")
    <br>
    <br>
    &nbsp;print("process s coordinates")
    <br>
    &nbsp;h=f-&gt;h;
    <br>
    &nbsp;h=where(h.eq.0, 1.e-14, h)
    <br>
    &nbsp;h2=h+f-&gt;hc;
    <br>
    &nbsp;h2inv=1./h2;
    <br>
    &nbsp;zeta=f-&gt;zeta
    <br>
    <br>
    &nbsp;; grid depth at rho-points
    <br>
    &nbsp;sc_r=f-&gt;sc_r
    <br>
    &nbsp;Cs_r=f-&gt;Cs_r
    <br>
    &nbsp;z_r = new( dimsizes(f-&gt;rho), float, default_fillvalue("float") )
    <br>
    &nbsp;nlevels=dimsizes(f-&gt;s_rho)
    <br>
    &nbsp;do k=0,nlevels-1
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; z0=1.e+16*sc_r(k)+Cs_r(k)*h;
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; dum1=(1 + z0 * h2inv)
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; dum2=z0 * h * h2inv
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; z_r(:,k,:,:)=tofloat(&nbsp;&nbsp;&nbsp; conform_dims(dimsizes(zeta),
    dum2, (/1,2/) )&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; + zeta * conform_dims(dimsizes(zeta), dum1,
    (/1,2/) )&nbsp;&nbsp;&nbsp; )
    <br>
    &nbsp;&nbsp;&nbsp; end do
    <br>
    &nbsp;z_r@standard_name="grid_point_depth_at_RHO-points";
    <br>
    &nbsp;z_r@long_name="grid point depth at RHO-points";
    <br>
    &nbsp;z_r@units=zeta@units;
    <br>
    &nbsp;z_r@axis="XYZT";
    <br>
    end if
    <br>
    <br>
    <br>
    topo=f-&gt;h * (-1) ; move to negative values
    <br>
    interp_topo = rcm2points (lat, lon, topo, cross_section_gc@gclat,
    cross_section_gc@gclon, 2)
    <br>
    <br>
    <br>
    ; PLOT
    <br>
    ; ========================
    <br>
    <br>
    <br>
    ; define cross-section
    <br>
    ;********************************
    <br>
    <br>
    if ( "xx" .eq. "xTruex" ) ; Turn on edges on
    <br>
    &nbsp;&nbsp;&nbsp; npts = 50
    <br>
    else
    <br>
    &nbsp;&nbsp;&nbsp; dum=dimsizes(lat)
    <br>
    &nbsp;&nbsp;&nbsp; d1=dum(0)
    <br>
    &nbsp;&nbsp;&nbsp; d2=dum(1)
    <br>
    &nbsp;&nbsp;&nbsp; d=sqrt(d1*d1+d2*d2)
    <br>
    &nbsp;&nbsp;&nbsp; npts&nbsp;&nbsp; =&nbsp; min ( (/ toint(d), 1000 /)) ; number of points in
    resulting cross-section
    <br>
    end if
    <br>
    cross_section_gc&nbsp;&nbsp;&nbsp;&nbsp; = gc_latlon(35.9,-6,35.89,-5.6,npts,-2)
    <br>
    points&nbsp;&nbsp; = ispan(0,npts-1,1)*1.0
    <br>
    <br>
    <br>
    ; create contour plot
    <br>
    ;********************************
    <br>
    print("create contour plot")
    <br>
    wks= gsn_open_wks("x11", "figure")
    <br>
    gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
    <br>
    <br>
    res&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = True&nbsp;&nbsp;&nbsp;&nbsp; ; Use plot options
    <br>
    <br>
    ; high-level control
    <br>
    res@gsnDraw&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = False&nbsp;&nbsp;&nbsp; ; dont draw
    <br>
    res@gsnFrame&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = False&nbsp;&nbsp;&nbsp; ; dont advance frame
    <br>
    res@gsnMaximize&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = True&nbsp;&nbsp;&nbsp;&nbsp; ; maximize plot in frame
    <br>
    res@gsnBoxMargin&nbsp;&nbsp;&nbsp;&nbsp; = 0.0
    <br>
    res@gsnLeftString = " "
    <br>
    res@gsnRightString = " "
    <br>
    res@gsnScale&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = False
    <br>
    res@trGridType&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = "curvilinear"; "TriangularMesh" ;
    <br>
    <br>
    ; xticks labels
    <br>
    nLabels=10
    <br>
    XBValues&nbsp;&nbsp;&nbsp; = toint( fspan(0,npts-1,nLabels) )
    <br>
    XBLabels&nbsp;&nbsp;&nbsp; = new(nLabels,"string")
    <br>
    do i=0,nLabels-1
    <br>
    &nbsp;&nbsp;&nbsp; x = cross_section_gc@gclon(XBValues(i))
    <br>
    &nbsp;&nbsp;&nbsp; y = cross_section_gc@gclat(XBValues(i))
    <br>
    &nbsp;&nbsp;&nbsp; XBLabels(i) = sprintf("%5.1f",
    y)+"~S~o~N~N"+"~C~"+sprintf("%5.1f", x)+"~S~o~N~E"
    <br>
    end do
    <br>
    <br>
    print(XBValues)
    <br>
    print(XBLabels)
    <br>
    <br>
    <br>
    res@tmXBMode&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = "Explicit"&nbsp;&nbsp;&nbsp; ; explicitly label x-axis
    <br>
    res@tmXBValues = XBValues
    <br>
    res@tmXBLabels = XBLabels
    <br>
    res@tmXBLabelFontHeightF = 0.01
    <br>
    <br>
    ; contours
    <br>
    res@cnFillOn&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = True
    <br>
    <br>
    if ( "xx" .eq. "xTruex" ) ; Turn on edges on
    <br>
    &nbsp;&nbsp; res@cnFillMode&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = "CellFill"
    <br>
    &nbsp;&nbsp; res@cnCellFillEdgeColor= "red" ; "grey22"&nbsp;&nbsp;&nbsp;&nbsp; ; "black
    <br>
    &nbsp;&nbsp; res@cnLinesOn&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = False&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; turn off countour lines
    <br>
    &nbsp;&nbsp; res@trGridType&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = "TriangularMesh" ;
    <br>
    end if
    <br>
    res@cnFillDrawOrder = "PreDraw"
    <br>
    res@cnLineDrawOrder = "PreDraw"
    <br>
    <br>
    <br>
    res@cnInfoLabelOn&nbsp;&nbsp; = False&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; turn off contour
    label
    <br>
    res@cnMaxLevelCount&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 21
    <br>
    varmin=min(var)
    <br>
    varmax=max(var)
    <br>
    maxlev=50
    <br>
    ;---Make sure contour levels are fixed for all time steps
    <br>
    mini=3.40282347e+38
    <br>
    maxi=-1.17549435e-38
    <br>
    i= 0
    <br>
    do while (i .le. 2 )
    <br>
    &nbsp; varmin=min( var(i+0,:,:,:) )
    <br>
    &nbsp; varmax=max( var(i+0,:,:,:) )
    <br>
    &nbsp; mini=min((/tofloat(mini), tofloat(varmin)/)) ; minimum value of
    mini,varmin array
    <br>
    &nbsp; maxi=min((/tofloat(maxi), tofloat(varmax)/)) ; maximum value
    <br>
    &nbsp; i=i+1
    <br>
    end do
    <br>
    mnmxint = nice_mnmxintvl(varmin, varmax, maxlev, False)
    <br>
    <br>
    res@cnLevelSelectionMode = "ManualLevels"
    <br>
    res@cnMinLevelValF&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = mnmxint(0) ; min( (<i
      class="moz-txt-slash"><span class="moz-txt-tag">/</span>mnmxint(0),
      mini<span class="moz-txt-tag">/</span></i>) )
    <br>
    res@cnMaxLevelValF&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = mnmxint(1)
    <br>
    res@cnLevelSpacingF&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = mnmxint(2)
    <br>
    <br>
    <br>
    ; label bar
    <br>
    res@lbLabelAutoStride&nbsp;&nbsp; = True&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; nice label bar label stride
    <br>
    res@gsnSpreadColors&nbsp;&nbsp;&nbsp;&nbsp; = True&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; use full range of colormap
    <br>
    <br>
    res@lbLabelBarOn = True
    <br>
    res@lbTitlePosition = "Bottom"
    <br>
    res@lbTitleDirection = "Across"
    <br>
    res@lbOrientation = "Vertical"
    <br>
    res@lbTopMarginF = 0.0
    <br>
    <br>
    ; titles
    <br>
    res@tiYAxisString&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = ("z_r" + " ["+z_r@units+"]")
    <br>
    res@gsnYAxisIrregular2Linear = True ; Y axis will be linearized
    <br>
    <br>
    <br>
    <br>
    ; another ressource for bathymetry
    <br>
    tres&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = True&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; plot mods desired
    <br>
    tres@gsnDraw&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = False&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; don t draw yet
    <br>
    tres@gsnFrame = False&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; don t advance frame yet
    <br>
    tres@gsnMaximize&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = True&nbsp;&nbsp;&nbsp;&nbsp; ; maximize plot in frame
    <br>
    tres@xyCurveDrawOrder = "PostDraw"
    <br>
    tres@gsnYRefLine&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = min(interp_topo)&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; create a
    reference line
    <br>
    tres@gsnAboveYRefLineColor = "tan";"gray75"&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; above ref
    line fill color
    <br>
    tres@trYReverse&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = True&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; reverse y axis
    <br>
    <br>
    i=0
    <br>
    do while (i .le. 2 )
    <br>
    <br>
    &nbsp;&nbsp;&nbsp; ; interpolate data along a great circle
    <br>
    &nbsp;&nbsp;&nbsp; ;************************************
    <br>
    &nbsp;&nbsp;&nbsp; print("interpolate data along the cross-section")
    <br>
    &nbsp;&nbsp;&nbsp; interp = new ((<i class="moz-txt-slash"><span
        class="moz-txt-tag">/</span>nlevels,npts<span
        class="moz-txt-tag">/</span></i>), double,
    default_fillvalue("double"))
    <br>
    &nbsp;&nbsp;&nbsp; z = new ((<i class="moz-txt-slash"><span class="moz-txt-tag">/</span>nlevels,npts<span
        class="moz-txt-tag">/</span></i>), double,
    default_fillvalue("double"))
    <br>
    &nbsp;&nbsp;&nbsp; do k = 0,nlevels-1
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; interp(k,:)&nbsp; = rcm2points_Wrap (lat, lon, var(i+0,k,:,:),
    cross_section_gc@gclat, cross_section_gc@gclon, 2)
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; z(k,:) = rcm2points_Wrap (lat, lon, z_r(i+0,k,:,:),
    cross_section_gc@gclat, cross_section_gc@gclon, 2)
    <br>
    &nbsp;&nbsp;&nbsp; end do
    <br>
    <br>
    &nbsp;&nbsp;&nbsp; interp!0 = "levels"
    <br>
    &nbsp;&nbsp;&nbsp; interp!1 = "points"
    <br>
    &nbsp;&nbsp;&nbsp; interp&amp;levels=vlev
    <br>
    <br>
    &nbsp;&nbsp;&nbsp; z!0 = "levels"
    <br>
    &nbsp;&nbsp;&nbsp; z!1 = "points"
    <br>
    &nbsp;&nbsp;&nbsp; z&amp;levels=vlev
    <br>
    <br>
    <br>
    &nbsp;&nbsp;&nbsp; copy_VarAtts(var,interp)&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; copy attributes
    <br>
    &nbsp;&nbsp;&nbsp; copy_VarAtts(z_r,z)&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; copy attributes
    <br>
    <br>
    &nbsp;&nbsp;&nbsp; ; plot section
    <br>
    &nbsp;&nbsp;&nbsp; ;************************************
    <br>
    <br>
    &nbsp;&nbsp;&nbsp; ; update title
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; t=time(i)
    <br>
    &nbsp;&nbsp; utc_date = cd_calendar(t, 0)
    <br>
    &nbsp;&nbsp; year&nbsp;&nbsp; = tointeger(utc_date(:,0))
    <br>
    &nbsp;&nbsp; month&nbsp; = tointeger(utc_date(:,1))
    <br>
    &nbsp;&nbsp; day&nbsp;&nbsp;&nbsp; = tointeger(utc_date(:,2))
    <br>
    &nbsp;&nbsp; hour&nbsp;&nbsp; = tointeger(utc_date(:,3))
    <br>
    &nbsp;&nbsp; minute = tointeger(utc_date(:,4))
    <br>
    &nbsp;&nbsp; sdate = sprinti("%0.2i/", day) + sprinti("%0.2i/",month)&nbsp; +
    sprinti("%0.4i", year)
    <br>
    &nbsp;&nbsp; stime = sprinti("%0.2ih", hour)&nbsp; + sprinti("%0.2i", minute)
    <br>
    &nbsp;&nbsp; titledate = ("time = " + sdate + " - " + stime)
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; res@tiMainString&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = (varname +" ~C~" + titledate)
    <br>
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; print out stats as subtiltle
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; statsline=("min="+sprintf("%5.2f", min(interp)) + "&nbsp;&nbsp;&nbsp; mean="
    + sprintf("%5.2f", avg(interp)) + "&nbsp;&nbsp;&nbsp; max=" + sprintf("%5.2f",
    max(interp)) +"&nbsp;&nbsp;&nbsp; ["+ interp@units +"]")
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; res@gsnCenterString = statsline
    <br>
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; plot sliced variable
    <br>
    <br>
    &nbsp;&nbsp;&nbsp; res@sfXArray&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = tofloat(conform_dims(dimsizes(z),
    points, 1) ); data for x-axis tms
    <br>
    &nbsp;&nbsp;&nbsp; res@sfYArray&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = z&nbsp; ; data for y-axis tms
    <br>
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; plot sliced variable
    <br>
    &nbsp;&nbsp;&nbsp; profile = gsn_csm_contour(wks,interp,res)
    <br>
    <br>
    <br>
    ; add terrain topography
    <br>
    ;************************************
    <br>
    topo_profile = gsn_csm_y(wks,interp_topo,tres)
    <br>
    overlay(profile, topo_profile)
    <br>
    <br>
    &nbsp;&nbsp;&nbsp; maximize_output(wks,True)
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; i=i+1
    <br>
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ; clean up variables
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; delete(profile)
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; delete(interp)
    <br>
    &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; delete(z)
    <br>
    end do
    <br>
    <br>
    end
  </body>
</html>