Re: Error in Stream Function script

From: Helen Parish <hparish_at_nyahnyahspammersnyahnyah>
Date: Tue, 7 Jul 2009 03:56:49 -0700

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