I am trying to plot the stream function for my data, using an NCL
script. I am getting the error below, but I do not see what is
causing it and how to get the NCL script to run.
I have included the NCL script that produces the error and the
Fortran routine that the script calls (which had to be customized for
my model from DZPSIDRV). I also include the statement I use to
compile the Fortran code.
Does anyone know what is causing this error and how to get the script
to run ?.
Thanks,
Helen.
a) The error:
Variable: v
Type: float
Total Size: 344448 bytes
86112 values
Number of Dimensions: 4
Dimensions and sizes: [time | 1] x [lev_p | 26] x [lat | 46] x [lon
| 72]
Coordinates:
time: [13710..13710]
lev_p: [91500..0.03]
lat: [ -90.. 90]
lon: [ 0.. 355]
Number Of Attributes: 4
units : m/s
long_name : Meridional wind
cell_method : time: mean
_FillValue : -9999
warning:Argument 4 of the current function or procedure was coerced
to the appropriate type and thus will not change if the function or
procedure modifies its value
warning:Argument 6 of the current function or procedure was coerced
to the appropriate type and thus will not change if the function or
procedure modifies its value
warning:Argument 7 of the current function or procedure was coerced
to the appropriate type and thus will not change if the function or
procedure modifies its value
warning:Argument 9 of the current function or procedure was coerced
to the appropriate type and thus will not change if the function or
procedure modifies its value
fatal:VENUSDRV: dimension size of dimension (3) of V must be equal to
the value of NTIM
fatal:Execute: Error occurred at or near line 66 in file
zonstreamtest.ncl
b) The NCL script:
external HELEN "./venus_mpsitest.so"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
; input file
diri = "./"
fili = "surfha.cam2.h0.0252-08-25-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)
meridstream = new ( (/ntim,klvl,nlat/) , typeof(vavg), getFillValue
(vavg))
HELEN::VENUSDRV
(mlon,nlat,klvl,ntim,v,lat,lev_p,PS,vavg@_FillValue,meridstream)
;***********************
; Create Plot
;***********************
wks = gsn_open_wks ("pdf", "surfha2520825stream" ) ;
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
do nt=0,ntim-1
res_at_gsnCenterString = "surfha2520825stream"
plot = gsn_csm_pres_hgt(wks, meridstream(nt,:,:), res ) ;
(lev,lat)
res_at_trYReverse = True
kl = 5
res_at_mpFillOn = False
res_at_mpGridAndLimbOn = True
res_at_mpGridLineDashPattern = 2
res_at_mpOutlineBoundarySets = "NoBoundaries"
res_at_mpCenterLonF = 180.
res_at_gsnCenterString = res_at_gsnCenterString+" p="+u&lev_p(kl)
end do
end
c) The Fortran routine VENUSDRV :
C NCLFORTSTART
SUBROUTINE VENUSDRV(MLON,NLAT,KLEV,NTIM,V,LAT,P,PS,VMSG,
+ ZMPSI)
IMPLICIT NONE
c NCL: zmpsi = zon_mpsi(v,lat,p,ps)
c INPUT
INTEGER MLON,NLAT,KLEV,NTIM
DOUBLE PRECISION V(NTIM,KLEV,NLAT,MLON),LAT(NLAT),P(KLEV),
+ PS(NTIM,NLAT,MLON),VMSG
c OUTPUT
DOUBLE PRECISION ZMPSI(NTIM,KLEV,NLAT)
C NCLEND
c LOCAL
DOUBLE PRECISION PTMP(2*KLEV+1),DP(2*KLEV+1),VVPROF(2*KLEV+1)
DOUBLE PRECISION VBAR(NTIM,KLEV,NLAT),VTMP(NTIM,KLEV,NLAT,MLON)
CALL VENUSZON(MLON,NLAT,KLEV,NTIM,V,LAT,P,PS,VMSG,ZMPSI,
+ PTMP,DP,VTMP,VBAR,VVPROF)
RETURN
END
c --
SUBROUTINE VENUSZON(MLON,NLAT,KLEV,NTIM,V,LAT,P,PS,VMSG,ZMPSI,
+ PTMP,DP,VTMP,VBAR,VVPROF)
IMPLICIT NONE
INTEGER MLON,NLAT,KLEV,NTIM
c DOUBLE PRECISION V(MLON,NLAT,KLEV,NTIM),LAT(NLAT),P(KLEV),
c + PS(MLON,NLAT,NTIM),VMSG
DOUBLE PRECISION V(NTIM,KLEV,NLAT,MLON),LAT(NLAT),P(KLEV),
+ PS(NTIM,NLAT,MLON),VMSG
c output
c DOUBLE PRECISION ZMPSI(NLAT,KLEV)
c DOUBLE PRECISION ZMPSI(KLEV,NLAT)
DOUBLE PRECISION ZMPSI(NTIM,KLEV,NLAT)
c temporary / local
DOUBLE PRECISION PTMP(0:2*KLEV),DP(0:2*KLEV),VVPROF(0:2*KLEV)
DOUBLE PRECISION VTMP(NTIM,KLEV,NLAT,MLON),VBAR(NTIM,KLEV,NLAT)
c DOUBLE PRECISION VTMP(MLON,NLAT,KLEV,NTIM),VBAR(NLAT,KLEV)
DOUBLE PRECISION PI,G,A,C,RAD,CON,PTOP,PBOT
INTEGER ML,NL,KL,NT,KNT,KFLAG
G = 9.80616D0
A = 6.37122D6
PI = 4.D0*ATAN(1.D0)
RAD = PI/180.D0
CON = 2.D0*PI*A/G
PTOP = 3.D0
PBOT = 9150000.D0
c calculate press at all levels [even half levels]
KNT = 0
DO KL = 1,2*KLEV - 1,2
KNT = KNT + 1
PTMP(KL) = P(KNT)
END DO
DO KL = 2,2*KLEV - 2,2
PTMP(KL) = (PTMP(KL+1)+PTMP(KL-1))*0.5D0
END DO
PTMP(0) = PTOP
PTMP(2*KLEV) = PBOT
c dp at all levels
DP(0) = 0.D0
DO KL = 1,2*KLEV - 1
DP(KL) = PTMP(KL+1) - PTMP(KL-1)
END DO
DP(2*KLEV) = 0.D0
c make copy; set any p > ps to msg
KNT = 0
DO ML = 1,MLON
DO NL = 1,NLAT
DO KL = 1,KLEV
DO NT = 1,NTIM
IF (P(KL).GT.PS(NT,NL,ML)) THEN
KNT = KNT + 1
VTMP(NT,KL,NL,ML) = VMSG
ELSE
VTMP(NT,KL,NL,ML) = V(NT,KL,NL,ML)
END IF
END DO
END DO
END DO
END DO
c compute zonal mean v using the vtmp variable
DO NT = 1,NTIM
DO NL = 1,NLAT
DO KL = 1,KLEV
KNT = 0
VBAR(NT,KL,NL) = 0.0D0
DO ML = 1,MLON
IF (VTMP(NT,KL,NL,ML).NE.VMSG) THEN
KNT = KNT + 1
VBAR(NT,KL,NL) = VBAR(NT,KL,NL) + VTMP(NT,KL,NL,ML)
END IF
END DO
IF (KNT.GT.0) THEN
VBAR(NT,KL,NL) = VBAR(NT,KL,NL)/DBLE(KNT)
ELSE
VBAR(NT,KL,NL) = VMSG
END IF
END DO
END DO
END DO
c compute mpsi at each latitude [reuse ptmp]
DO NT = 1,NTIM
DO NL = 1,NLAT
C = CON*COS(LAT(NL)*RAD)
PTMP(0) = 0.0D0
DO KL = 1,2*KLEV
PTMP(KL) = VMSG
END DO
DO KL = 0,2*KLEV,2
VVPROF(KL) = 0.0D0
END DO
KNT = 0
DO KL = 1,2*KLEV - 1,2
KNT = KNT + 1
VVPROF(KL) = VBAR(NT,KNT,NL)
END DO
c integrate from top of atmosphere down for each level where vbar
c is not missing
DO KL = 1,2*KLEV - 1,2
KFLAG = KL
IF (VVPROF(KL).EQ.VMSG) GO TO 10
PTMP(KL+1) = PTMP(KL-1) - C*VVPROF(KL)*DP(KL)
END DO
10 CONTINUE
c impose lower boundary condition to ensure the zmpsi is 0
c at the bottom boundary
PTMP(KFLAG+1) = -PTMP(KFLAG-1)
c streamfunction is obtained as average from its' values
c at the intermediate half levels. The minus sign before
c PTMP in the last loop is to conform to a CSM convention.
DO KL = 1,KFLAG,2
PTMP(KL) = (PTMP(KL+1)+PTMP(KL-1))*0.5D0
END DO
KNT = 0
DO KL = 1,2*KLEV - 1,2
KNT = KNT + 1
IF(PTMP(KL).NE.VMSG) THEN
ZMPSI(NT,KNT,NL) = -PTMP(KL)
ELSE
ZMPSI(NT,KNT,NL) = PTMP(KL)
END IF
END DO
END DO
END DO
RETURN
END
d) Statement used to compile the Fortran code :
WRAPIT -g95 -L /Users/helenparish/g95-install/lib/gcc-lib/powerpc-
apple-darwin6.8/4.0.3 -l f95 venus_mpsitest.f
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Jul 02 2009 - 05:43:27 MDT
This archive was generated by hypermail 2.2.0 : Thu Jul 02 2009 - 11:39:38 MDT