Mary,
I have modified my NCL script to pass the variables as doubles to the
Fortran routine, and I have reordered the scripts in the Fortran
routine so they are the opposite way round from those in the NCL
script, as you suggested This removes some of the errors. However, in
spite of trying many different approaches, I cannot get the routine
"gsn_csm_pres_hgt" to plot the results. It keeps coming up with the
same error, as given below. My current NCL script and Fortran routine
are also given below.
Thanks,
Helen.
1). Error message:
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
(0) gsn_csm_pres_hgt: Fatal: The first dimension of the input
data must
(0) have a coordinate variable called 'lev.'
(0) Cannot create plot.
fatal:Illegal right-hand side type for assignment
fatal:Execute: Error occurred at or near line 104 in file send.ncl
2). Current 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)
vavgdouble = flt2dble(vavg)
vdouble = flt2dble(v)
levpdouble = flt2dble(lev_p)
psdouble = flt2dble(PS)
meridstream = new ( (/ntim,klvl,nlat/) , double)
meridstream!0 = "time"
meridstream!1 = "lev"
meridstream!2 = "lat"
HELEN::VENUSDRV
(mlon,nlat,klvl,ntim,vdouble,lat,levpdouble,psdouble,vavg@_FillValue,mer
idstream)
meridstreamfloat = dble2flt(meridstream)
meridstreamfloat!0 = "time"
meridstreamfloat!1 = "lev"
meridstreamfloat!2 = "lat"
meridstream_reorder = new ( (/klvl,nlat,ntim/) , float)
meridstream_reorder!0 = "lev"
meridstream_reorder!1 = "lat"
meridstream_reorder!2 = "time"
meridstream_reorder = meridstreamfloat(lev | :, lat | :, time
| :) ; (lev,lat,time)
;***********************
; 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
nt = 0
res_at_gsnCenterString = "surfha2520825stream"
plot = gsn_csm_pres_hgt(wks, meridstream_reorder(:,:,nt),
res ) ; (lev,lat,time)
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
3). Current Fortran routine:
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(MLON,NLAT,KLEV,NTIM),LAT(NLAT),P(KLEV),
+ PS(MLON,NLAT,NTIM),VMSG
DOUBLE PRECISION ZMPSI(NLAT,KLEV,NTIM)
C NCLEND
c LOCAL
DOUBLE PRECISION PTMP(2*KLEV+1),DP(2*KLEV+1),VVPROF(2*KLEV+1)
DOUBLE PRECISION VBAR(NLAT,KLEV,NTIM),VTMP(MLON,NLAT,KLEV,
+ NTIM)
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
DOUBLE PRECISION V(MLON,NLAT,KLEV,NTIM),LAT(NLAT),P(KLEV),
+ PS(MLON,NLAT,NTIM),VMSG
DOUBLE PRECISION ZMPSI(NLAT,KLEV,NTIM)
c temporary / local
DOUBLE PRECISION PTMP(0:2*KLEV),DP(0:2*KLEV),VVPROF(0:2*KLEV)
DOUBLE PRECISION VTMP(MLON,NLAT,KLEV,NTIM),
+ VBAR(NLAT,KLEV,NTIM)
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(ML,NL,NT)) THEN
KNT = KNT + 1
VTMP(ML,NL,KL,NT) = VMSG
ELSE
VTMP(ML,NL,KL,NT) = V(ML,NL,KL,NT)
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(NL,KL,NT) = 0.0D0
DO ML = 1,MLON
IF (VTMP(ML,NL,KL,NT).NE.VMSG) THEN
KNT = KNT + 1
VBAR(NL,KL,NT) = VBAR(NL,KL,NT) + VTMP(ML,NL,KL,NT)
END IF
END DO
IF (KNT.GT.0) THEN
VBAR(NL,KL,NT) = VBAR(NL,KL,NT)/DBLE(KNT)
ELSE
VBAR(NL,KL,NT) = 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(NL,KNT,NT)
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(NL,KNT,NT) = -PTMP(KL)
ELSE
ZMPSI(NL,KNT,NT) = PTMP(KL)
END IF
END DO
END DO
END DO
RETURN
END
On Jul 2, 2009, at 10:25 AM, Mary Haley wrote:
> Helen,
>
> The error messages are telling you quite a bit here. For one:
>
> 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
>
> Your variables are probably coming in as floats, and you are trying to
> pass these floats to a Fortran routine that's expecting doubles.
>
> You should convert your floats to doubles if you want to get rid
> of these warnings.
>
> Secondly, the error:
>
> fatal:VENUSDRV: dimension size of dimension (3) of V must be
> equal to the
> value of NTIM
>
> indicates that your arrays do not have the correct dimensions. I
> believe this is b/c you are passing the arrays in the wrong order.
>
> Fortran and NCL order their arrays differently (column order versus
> row order). This means if you have an array on the Fortran
> side that is ntime x nlev x nlat x nlon, then in the NCL
> script, that array must be ordered nlon x nlat x nlev x ntime.
>
> In your case, the NCL's "v" is
>
> time x lev_p lat x lon
>
> and the Fortran routine is expecting
>
> V(NTIM,KLEV,NLAT,MLON)
>
> You will either need to reorder the arrays on the Fortran side,
> or on the NCL side:
>
> vreord = v(lon|:,lat|:,lev_P|:,time|:)
>
> You have to reorder "meridstream" and probably PS as well:
>
> PSreord = PS(lon|:,lat|:,time|:)
> meridstream = new ( (/nlat,klvl,ntim/) , typeof(vavg),
> getFillValue(vavg))
> HELEN::VENUSDRV (mlon,nlat,klvl,ntim,vreord,lat,lev_p,PSreord, \
> vavg@_FillValue,meridstream)
>
>
> You could instead reorder the arrays on the Fortran side. This will
> probably run faster than reordering them on the NCL side.
>
> --Mary
>
>
>
Received on Tue Jul 07 2009 - 04:56:49 MDT
This archive was generated by hypermail 2.2.0 : Tue Jul 07 2009 - 11:13:18 MDT