Re: Error in Stream Function script - How to Overlay ?

Mary,

After printing out values from the Fortran stream function routine I
realised that the pressures are defined in the opposite sense in the
NCL script and the Fortran routine. Now I have reversed them for the
Fortran routine the plots appear more realistic. I had to (I believe
correctly) reverse the direction of the pressure levels which are
output from the Fortran routine before plotting the results (so that
the surface pressure is at the bottom of the plot). What I want to do
now is to:

1. Turn the stream function plots from colour contours to black and
white contour lines with numbers on them, and

2. Overlay these black and white stream function contour lines on top
of colour contour plots of, for example, "U" the zonal wind.

I still want the overlaid plots to be a function of pressure/height
and latitude, so I assume I will need to continue to use
"gsn_csm_pres_hgt" ?. Do you know how I can do this ?. I enclose my
current NCL script below.

Thanks,
Helen.

Here is my current NCL script:

external HELEN "./venus_mpsitest.so"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
                                        ; input file
   diri = "./"
   fili = "zeusrand09125.cam2.h0.0172-01-05-00000.nc"
   fi = addfile (diri+fili, "r")
                                        ; output file
   diro = "./"
   filo = "helen.nc"
   system ("/bin/rm -f "+diro+filo) ; remove any pre-exist file
   fo = addfile (diro+filo, "c")
                                        ; add any file attributes
   fo_at_title = fi_at_title + ": Selected Variables at pressure levels"
   fo_at_history = systemfunc ("date")
   fo_at_source = fi_at_source
   fo_at_case = fi_at_case
   fo_at_Conventions = fi_at_Conventions

   Var = (/ "T" , "Q", "U", "V"/) ; select variable to be
interpolated.
   nVar = dimsizes (Var)
                                        ; desired output levels
lev_p = (/ 91500., 90000., 87500., 80000., 67500., 50000., 30000.,
20000., 10000., 5000., 3000., 1500., 750., 390., 200., 100., 50.,25.,
14., 7., 3.5, 1.9, 0.95, 0.4, 0.1, 0.03 /)
   lev_p!0 = "lev_p" ; variable and dimension name
the same
   lev_p&lev_p = lev_p ; create coordinate variable
   lev_p_at_long_name = "pressure" ; attach some attributes
   lev_p_at_units = "hPa"
   lev_p_at_positive = "down"

   hyam = fi->hyam ; read hybrid info
   hybm = fi->hybm
   PS = fi->PS
   P0mb = 0.01*fi->P0

   do n=0,nVar-1 ; loop over the variables
      X = fi->$Var(n)$
      Xp = vinth2p (X, hyam, hybm, lev_p ,PS, 1, P0mb, 2, True)
      copy_VarAtts(X, Xp)
      fo->$Var(n)$ = Xp ; write to netCDF file
      print (Var(n)+": interpolated and written to netCDF")
   end do

end
;***********************
; zonal.ncl
;***********************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;***********************
begin
   f = addfile ("helen.nc","r")
   print(f)
   v = f->V
   lat = f->lat
   printVarSummary( v )

   dimu = dimsizes( v )
   ntim = dimu(0)
   klvl = dimu(1)
   nlat = dimu(2)
   mlon = dimu(3)

  vavg = dim_avg_Wrap(v)
  vavgdouble = flt2dble(vavg)

  levpreverse = lev_p(25:0)

  vdouble = flt2dble(v)
  levpdouble = 100.0*flt2dble(levpreverse)
  psdouble = 100.0*flt2dble(PS)

meridstream = new ( (/ntim,klvl,nlat/) , double)
meridstream!0 = "time"
meridstream!1 = "lev"
meridstream!2 = "lat"
meridstream&time = v&time
meridstream&lev = levpreverse
meridstream&lat = v&lat

  HELEN::VENUSDRV
(mlon,nlat,klvl,ntim,vdouble,lat,levpdouble,psdouble,vavg@_FillValue,mer
idstream)

  meridstream_reorder = meridstream(lev | :, lat | :, time
| :) ; (lev,lat,time)

;***********************
; Create Plot
;***********************
   wks = gsn_open_wks ("pdf", "zeusrand.172.stream" ) ;
open workstation
   gsn_define_colormap(wks,"rainbow") ; choose colormap

   res = True
   res_at_cnFillOn = True
   res_at_lbLabelAutoStride = True
   res_at_gsnMaximize = True ; if [ps, eps, pdf] make large
   res_at_gsnSpreadColors = True ; span color map

    nt = 0

      res_at_gsnCenterString = "zeusrand.172.stream"
      plot = gsn_csm_pres_hgt(wks, meridstream_reorder(:,:,nt),
res ) ; (lev,lat,time)

      res_at_trYReverse = True

end

On Jul 8, 2009, at 8:38 AM, Mary Haley wrote:

>
> Helen,
>
> I think your new script looks okay. It doesn't make sense that
> you have to attach the coordinate arrays before calling the Fortran
> routine, but it also shouldn't hurt anything.
>
> If the plot is not making sense, maybe you have a missing value
> set incorrectly?
>
> Make sure that the _FillValue attribute is set if there are
> potential missing values in your data.
>
> --Mary
>
Received on Thu Jul 09 2009 - 05:29:43 MDT