Re: Error in Stream Function script

Helen,

Your variable has named dimensions, but you haven't attached any
coordinate arrays. The coordinate arrays are necessary in order for
gsn_csm_pres_hgt to work properly.

You have:

   meridstreamfloat = dble2flt(meridstream)
   meridstreamfloat!0 = "time"
   meridstreamfloat!1 = "lev"
   meridstreamfloat!2 = "lat"

This is only assigning dimension names. You need to further
assign coordinate arrays:

   meridstreamfloat&lev = levpdouble
   meridstreamfloat&time = ???
   meridstreamfloat&lat = ???

I have "???" for the last two, because I don't know what you want
these set to. Maybe copy the ones from "v":

   meridstreamfloat&lat = v&lat
   meridstreamfloat&time = v&time

To reorder "meridstreamfloat" you only need one line. The metadata
(named dimensions, attributes, coordinate arrays) will be copied
automatically:

    meridstream_reorder = meridstreamfloat(lev | :, lat | :, time | :)

--Mary

On Tue, 7 Jul 2009, Helen Parish wrote:

> 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
>>
>>
>>
>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Jul 07 2009 - 08:41:26 MDT