Incompatible dimensions in Zonal_mpsi

From: Helen Parish <hparish_at_nyahnyahspammersnyahnyah>
Date: Fri, 6 Feb 2009 23:59:39 -0800

I have been trying to convert zonal_mpsi for a different number of
vertical levels (with different top and bottom pressure level
limits). The version of routine zonal_mpsi which I have, and which I
have been trying to modify, assumes that the input meridional wind
has only 3 dimensions (latitude, longitude and level). However, the
meridional winds in my netcdf output files have 4 dimensions, and
include time.

I am wondering if there is an alternative version of the zonal_mpsi
code which includes the time variable?. If I have an example it would
hopefully be much easier to convert this into a usable form.

Particularly confusing is the fact that the version of zonal_mpsi
which I am trying to convert, has variable ZMPSI defined as ZMPSI
(NLAT, KLEV). However, in the online NCL documentation, ZMPSI is
defined as ZMPSI (KLEV,NLAT). I am not sure where exactly to put the
time into the dimensions of each variable, which variables to put
time into, and which do loops I need to include time within. I keep
making adjustments to the routine, and each time I do I come up with
a slightly different error in which the dimensions are inconsistent
between the NCL routine which reads in the data and the routine
zonal_mpsi. It seems to be this inconsistency of dimensions which is
causing problems. Perhaps there is another version available in which
the variable "time" is included ?.

I include below the partially converted version of the routine
zonal_mpsi which I have been attempting to modify, and following
that, the ncl script which I am using to read in the data.

Thanks,
Helen.

1) Here is the partially converted version of zonal_mpsi:

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
c OUTPUT
       DOUBLE PRECISION ZMPSI(NLAT,KLEV)
C NCLEND
c LOCAL
       DOUBLE PRECISION PTMP(2*KLEV+1),DP(2*KLEV+1),VVPROF(2*KLEV+1)
       DOUBLE PRECISION VBAR(NLAT,KLEV),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
c output
       DOUBLE PRECISION ZMPSI(NLAT,KLEV)
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)

       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 NL = 1,NLAT
           DO KL = 1,KLEV
                 DO NT = 1,NTIM
               KNT = 0
               VBAR(NL,KL) = 0.0D0
               DO ML = 1,MLON
                   IF (VTMP(ML,NL,KL,NT).NE.VMSG) THEN
                       KNT = KNT + 1
                       VBAR(NL,KL) = VBAR(NL,KL) + VTMP(ML,NL,KL,NT)
                   END IF
               END DO
               IF (KNT.GT.0) THEN
                   VBAR(NL,KL) = VBAR(NL,KL)/DBLE(KNT)
               ELSE
                   VBAR(NL,KL) = VMSG
               END IF
             END DO
           END DO
       END DO

c compute mpsi at each latitude [reuse ptmp]

       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)
           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) = -PTMP(KL)
               ELSE
                  ZMPSI(NL,KNT) = PTMP(KL)
               END IF
           END DO
       END DO

       RETURN
       END

2) NCL scrip that reads in the data :

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.,
20., 15.,10., 8., 5.,3.,1. /)
;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.01 /)
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 = (/ 91500., 90000., 87500., 80000., 67500., 50000., 30000.,
20000., 10000., 5000., 3000., 1500., 750., 390., 200., 100., 50.,25.,
14., 7., 3.5, 2.0 /)
   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)
; nlon = dimu(2)
   nlat = dimu(2)
   mlon = dimu(3)

;; meridstream = zonal_mpsi_Wrap(v,lat,lev_p,PS)
; meridstream = zonal_mpsi_Wrap(v(:,::-1,:,:),lat,lev_p(::-1),PS)

; meridstream = new ( (/klvl,nlat,mlon/) , typeof(v), getFillValue(v))
;; meridstream = new ( (/ntim,klvl,nlat,mlon/) , typeof(v),
getFillValue(v))
;;; meridstream = new ( (/nlat,klvl/) , typeof(v), getFillValue(v))
  meridstream = new ( (/nlat,klvl/) , float)
; HELEN::DZONPMPSI(mlon,nlat,klvl,v,lat,P,PS,v@_FillValue,meridstream)
; HELEN::VENUSZON(mlon,nlat,klvl,v,lat,P,PS,v@_FillValue,meridstream)
;;; HELEN::VENUSDRV
(mlon,nlat,klvl,ntim,v,lat,P,PS,v@_FillValue,meridstream)
  HELEN::VENUSDRV
(mlon,nlat,klvl,ntim,v,lat,lev_p,PS,v@_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, u(nt,:,:,40), res ) ; (lev,lat)
       plot = gsn_csm_pres_hgt(wks, meridstream(nt,:,:), res ) ;
(lev,lat)
; plot = gsn_csm_pres_hgt_streamline(wks, uavg(nt,:,:),
res ) ; (lev,lat)
; plot = gsn_csm_contour (wks, uavg(t,:,:), res ) ; (lev,lat)

      res_at_trYReverse = True
; plot = gsn_csm_contour (wks, u(nt,:,:,0), res ) ; (lev,lat)
; plot = gsn_csm_contour (wks, u(nt,:,0,:), res ) ; (lev,lon)

      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)
; plot = gsn_csm_contour_map_ce (wks, u(nt,kl,:,:), res )

   end do

end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sat Feb 07 2009 - 00:59:39 MST

This archive was generated by hypermail 2.2.0 : Thu Feb 12 2009 - 10:16:44 MST