Re: Fwd: Incompatible dimensions in Zonal_mpsi

Helen

[1]
NCL is written in C in which the slowest varying dimension
is the leftmost dimension and the fastest varying dimension
is the rightmost. Fortran is the excct opposite. On "paper"
it may seem that the arrays must be reordered. This is not true.
xncl(time,lev,lat,lon) maps directly into xfor(lon,lat,lev,time)
Further, NCL/C is zero-based while fortran is one-based.

Which language is Right/Wrong .... neither. It is arbitrary.
Matlab and IDL have different subscript bases and one is
row-major like NCL/C while the other is column major
like fortran.

Your NCL/netCDF "ZMPSI (KLEV, NLAT)" has the exact same memory
ordering as fortran's "ZMPSI (NLAT, KLEV)"

[2]
The fortran ZMPSI (NLAT, KLEV) should be expanded to
ZMPSI (NLAT, KLEV, NTIM)

The do loops should be changed appropriately.

Good luck

Mary Haley wrote:
> Is this the one you're talking about?
>
> Begin forwarded message:
>
>> *From: *Helen Parish <hparish_at_ess.ucla.edu <mailto:hparish_at_ess.ucla.edu>>
>> *Date: *February 7, 2009 12:59:39 AM MST
>> *To: *ncl-talk_at_ucar.edu <mailto:ncl-talk_at_ucar.edu>
>> *Cc: *hparish_at_ess.ucla.edu <mailto:hparish_at_ess.ucla.edu>
>> *Subject: **[ncl-talk] Incompatible dimensions in Zonal_mpsi*
>>
>> 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
>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sun Feb 08 2009 - 17:58:12 MST