[ferret_users] error producing stereographic plots of 4-D variables

[Hella Riede <hella.riede@xxxxxxx>]

Dear ferreters,

when running mp_stereographic_north from within a repeat loop over L, I get the following error and FERRET crashes - always after the 4th round/plot, independent of the start value of L:

in FERRET v6.9: "At line 143 of file tm_world_ax.F Fortran runtime error: Array reference out of bounds for array 'line_mem', lower bound of dimension 1 exceeded (-839 < 1)"

in pyFERRET v6.9: "**ERROR: illegal limits: "(MP_Y[J=@MAX] + MP_..." does not exist at Y=25N:90N Axis extremes are Y=0 DEFINE VARIABLE/quiet mp_test = 50 - `(mp_y[j=@max] + mp_y[j=@min])/2` Command file, command group, or REPEAT execution aborted"

I call the script from the ferret shell: repeat/l=10:40 (go polar_plot_example.jnl; pause)

Background info and data structures below.


Thank you in advance,
Hella



= BACKGROUND =

My script is based on the excellent help page http://ferret.pmel.noaa.gov/Ferret/faq/polar-stereographic-plots/. Difference: I work with 4-D data (3 space, 1 time). The vertical axis is in model levels, which I transform to pressure levels using ZAXREPLACE. After this transformation, the first stereographic plot at a specific Z always works, but any subsequent plots hang with just a white canvas, no matter if inside REPEAT, called from the same ferret shell, or called for the same L that worked in the first plot. The analogous plots using the native vertical model levels work fine multiple times in a row.

I found a related email from 2004, http://www.pmel.noaa.gov/maillists/tmap/ferret_users/fu_2004/msg00687.html:
"To make stereographic plots we need to "set grid my_variable" & this will load the information about abstract axis "K" of the variable along with other axis. In mp_graticule.jnl (ferret/go dir.) one repeat loop is defined as "repeat/k=2:5". So when we use a variable with the size of abstract axis "K" as 1 (as in your case [Z=1]) Ferret will complain that the size of abstract axis "K" set by the "set grid my_variable" is 1 ..so any value of "k" greater than 1 ( set by the `k` in DEFINE VARIABLE/quiet mp_xaxdel = mp_xdel/`k`; inside the  repeat loop) is not possible. You can get rid of this problem as described below."

I simply decided to get rid of my vertical grid altogether before using mp_stereographic_north: I write the desired horizontal slab to netCDF, reload it into ferret, and then use reshape to get rid of the K definition. (Getting rid of K on the fly without an external file in-between took way too much memory.)


= DATA STRUCTURE  =

== 4-D dataset  ==

  name     title                             I         J         K         L         M         N
  LEV_BNDS                                  1:2       ...       1:47      ...       ...       ...
        (invalid axis bounds)
  A        hybrid A coefficient at layer m  ...       ...       1:47      ...       ...       ...
  PS       Surface Pressure                 1:320     1:160     ...       1:240     ...       ...
  P0       reference pressure               ...       ...       ...       ...       ...       ...
  B        hybrid B coefficient at layer m  ...       ...       1:47      ...       ...       ...
  B_BNDS                                    1:2       ...       1:47      ...       ...       ...
  VMR_CO                                    1:320     1:160     1:47      1:240     ...       ...
  A_BNDS                                    1:2       ...       1:47      ...       ...       ...

pressure axis (hPa) built by
let pressure_hpa = (($pressure_hyam) + ($pressure_hybm) * ($pressure_psurf)) / 100
define axis/from/name=pressure_z/depth/unit=hPa (($pressure_hyam) + ($pressure_hybm) * ($pressure_pref)) / 100


== Extracted slab  ==

  name     title                             I         J         K         L         M         N
  VAR_PRESS
           ZAXREPLACE(VMR_CO,PRESSURE_HPA,  1:320     1:160     1:1       1:1       ...       ...


== Grid of variable var_p (no Z axis) ==

     GRID (G013)
  name       axis              # pts   start                end
  %%        LONGITUDE          320mr   0E                   1.125W
  %%        LATITUDE           160 i   0                    0
  normal    Z
  normal    T
  normal    E
  normal    F


\can mode veri

!can var/all
!can sym/all
!can dat/all

say "region:"
show region
say

set mem/size=600

def sym bottom_lon = -60 ! for map plot, -90 -> 90degW at bottom
def sym myvar = vmr_co ! sst
def sym plevel = 200 ! hPa

use my_huge_4D_dataset.nc 
!USE coads_climatology

!!! default font
PPL dfltfnt AS !CR

!!! vertical pressure axis
go pressure a b ps
let var_press = ZAXREPLACE(($myvar),pressure_hpa,Z[gz=pressure_z])
!let var_press = ($myvar)

say; say "saving var ($myvar) at pressure ($plevel) hPa to nc file"
!let var_press_tmp = var_press[z=($plevel)]
save/clob/file=temp_polar_plot.nc var_press[z=($plevel)]
!cancel data 1
cancel var/all
use temp_polar_plot.nc

say; say "generating 'flat' variable without z axis"
let xy_stamp  = x[gx=var_press]*0 + y[gy=var_press]*0
let var_p     = RESHAPE(var_press,xy_stamp)

say; say "setting grid and region"
SET grid var_p
SET REGION/X=0:360/Y=25:90

say; say "setting up map projection and aspect ratio"
GO mp_stereographic_north ($bottom_lon) 50
GO mp_aspect
LET masked_var = var_p * mp_mask

!!! Plot data, add land mask and land outline
FILL/NOAXES/NOLAB/lev=(-inf)(0,150,5)(inf) masked_var*1e9, x_page, y_page
GO mp_land ! land outline

!!! Mark latitude and longitude with a graticule
GO mp_graticule

can reg/all
can var/all
!can dat 2


exit