Hello,
I am trying to make vertical cross sections using a WRF output file located in analysis cluster Mirage. Everything works fine if I run this script on my computer. But for some reason it doesn't work in Mirage. I would for the code to work in Mirage since I will need to analyze big WRF data files. I won't be able to do it in my computer due to the low memory allocations. Is there something I can do? Any help will be greatly appreciated. Here is the ncl code.
Thanks,
Cecille
; Example script to produce plots for a WRF real-data run,
; with the ARW coordinate dynamics option.
; Plot data on a cross section
; This script will plot data from a a given point A to point B
; Vertical coordinate is height
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "./gsn_code.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "./WRFUserARW.ncl"
begin
scheme = "MOR"
city = "Jasper"
dir = "ICCP_project/"
period = "fut"
grid = "1km"
path = "/glade/scratch/cvillanu/Morrison/"
;
; The WRF ARW input file.
; This needs to have a ".nc" appended, so just do it.
a = addfile(path+"wrfout_d01_0001-01-01_00:00:00_MOR_1km.nc","r")
; We generate plots, but what kind do we prefer?
; type = "x11"
type = "pdf"
; type = "ps"
; type = "ncgm"
wks = gsn_open_wks(type,"W_vcross_1km")
; Set some basic resources
res = True
res@MainTitle = "REAL-TIME WRF"
pltres = True
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FirstTime = True
times = wrf_user_list_times(a) ; get times in the file
ntimes = dimsizes(times) ; number of times in the file
mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
nd = dimsizes(mdims)
;---------------------------------------------------------------
; do it = 0,ntimes-1,2 ; TIME LOOP
do it = 0,ntimes-1
; do it = 5, 13
print("Working on time: " + times(it) )
res@TimeLabel = times(it) ; Set Valid time to use on plots
tc = wrf_user_getvar(a,"tc",it) ; T in C
rh = wrf_user_getvar(a,"rh",it) ; relative humidity
z = wrf_user_getvar(a, "z",it) ; grid point height
w = wrf_user_getvar(a,"wa",it) ; vertical velocity
tk = wrf_user_getvar(a,"tk",it)
p = wrf_user_getvar(a,"PB",it)
rgas = 287.04 ; J Kg-1 K-1
qv = wrf_user_getvar(a, "QVAPOR",it) ;
if(isfilevar(a,"QCLOUD"))
qc = wrf_user_getvar(a, "QCLOUD",it)
end if
if(isfilevar(a,"QRAIN"))
qr = wrf_user_getvar(a, "QRAIN",it)
end if
if(isfilevar(a,"QSNOW"))
qs = wrf_user_getvar(a, "QSNOW",it)
end if
if(isfilevar(a,"QICE"))
qi = wrf_user_getvar(a, "QICE",it)
end if
if(isfilevar(a,"QGRAUP"))
qg = wrf_user_getvar(a, "QGRAUP",it)
end if
rho = p/(rgas*tk*(1+0.61*qv))
PWC = qr*rho
LWC = (qr+qc)*rho
if ( FirstTime ) then ; get height info for labels
zmin = 0.
zmax = max(z)/1000.
nz = floattoint(zmax/2 + 1)
FirstTime = False
end if
;---------------------------------------------------------------
do ip = 1, 24 ; we are doing 3 plots, specifying different start and end points
opts = True ; setting start and end times
plane = new(4,float)
if(ip .eq. 1) then
plane = (/ 0,16, 100,16 /) ; start x;y & end x;y point
end if
if(ip .eq. 2) then
plane = (/ 0,17, 100,17 /) ; start x;y & end x;y point
end if
if(ip .eq. 3) then
plane = (/ 0,18, 100,18 /) ; start x;y & end x;y point
end if
if(ip .eq. 4) then
plane = (/ 0,19, 100,19 /) ; start x;y & end x;y point
end if
if(ip .eq. 5) then
plane = (/ 0,20, 100,20 /) ; start x;y & end x;y point
end if
if(ip .eq. 6) then
plane = (/ 0,21, 100,21 /) ; start x;y & end x;y point
end if
if(ip .eq. 7) then
plane = (/ 0,22, 100,22 /) ; start x;y & end x;y point
end if
if(ip .eq. 8) then
plane = (/ 0,23, 100,23 /) ; start x;y & end x;y point
end if
if(ip .eq. 9) then
plane = (/ 0,24, 100,24 /) ; start x;y & end x;y point
end if
if(ip .eq. 10) then
plane = (/ 0,25, 100,25 /) ; start x;y & end x;y point
end if
if(ip .eq. 11) then
plane = (/ 0,26, 100,26 /) ; start x;y & end x;y point
end if
if(ip .eq. 12) then
plane = (/ 0,27, 100,27 /) ; start x;y & end x;y point
end if
if(ip .eq. 13) then
plane = (/ 0,28, 100,28 /) ; start x;y & end x;y point
end if
if(ip .eq. 14) then
plane = (/ 0,29, 100,29 /) ; start x;y & end x;y point
end if
if(ip .eq. 15) then
plane = (/ 0,30, 100,30 /) ; start x;y & end x;y point
end if
if(ip .eq. 16) then
plane = (/ 0,31, 100,31 /) ; start x;y & end x;y point
end if
if(ip .eq. 17) then
plane = (/ 0,32, 100,32 /) ; start x;y & end x;y point
end if
if(ip .eq. 18) then
plane = (/ 0,33, 100,33 /) ; start x;y & end x;y point
end if
if(ip .eq. 19) then
plane = (/ 0,34, 100,34 /) ; start x;y & end x;y point
end if
if(ip .eq. 20) then
plane = (/ 0,35, 100,35 /) ; start x;y & end x;y point
end if
if(ip .eq. 21) then
plane = (/ 0,36, 100,36 /) ; start x;y & end x;y point
end if
if(ip .eq. 22) then
plane = (/ 0,37, 100,37 /) ; start x;y & end x;y point
end if
if(ip .eq. 23) then
plane = (/ 0,38, 100,38 /) ; start x;y & end x;y point
end if
if(ip .eq. 24) then
plane = (/ 0,39, 100,39 /) ; start x;y & end x;y point
end if
gsn_define_colormap(wks, "BlWhRe") ; Define a new color map
; Vertical Velocity
w_plane = wrf_user_intrp3d(w,z,"v",plane,0.,opts)
; w_plane = 100.*w_plane
;res@TimeLabel = it
;res@vpWidthF = plot_width
;res@vpHeightF = plot_height
opts_w = res
;opts_w@FieldTitle = w@description
opts_w@UnitLabel = "m/s"
; opts_w@PlotLevelID = 0.001*height + " km"
opts_w@cnFillOn = True
opts_w@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
opts_w@cnMinLevelValF = -38 ; set min contour level
opts_w@cnMaxLevelValF = 38. ; set max contour level
opts_w@cnLevelSpacingF = 2 ; set contour spacing
; Basic Options for all cloud plots
opts_clouds = res
opts_clouds@cnFillOn = True
opts_clouds@UnitLabel = "g/kg"
if (isvar("qv"))
qv_plane = wrf_user_intrp3d(qv,z,"v",plane,0.,opts)
qv_plane = qv_plane*1000.
opts_qv = opts_clouds
end if
if (isvar("qc"))
qc_plane = wrf_user_intrp3d(qc,z,"v",plane,0.,opts)
qc_plane = qc_plane*1000.
opts_qc = opts_clouds
end if
if (isvar("qr"))
qr_plane = wrf_user_intrp3d( qr,z,"v",plane,0.,opts)
qr_plane = qr_plane*1000.
opts_qr = opts_clouds
end if
if (isvar("PWC"))
; opts_clouds@UnitLabel = "Precipitation Mass (g/m3)"
PWC_plane = wrf_user_intrp3d( PWC,z,"v",plane,0.,opts)
PWC_plane = PWC_plane*1000.
opts_PWC = opts_clouds
end if
if (isvar("LWC"))
; opts_clouds@UnitLabel = "Liquid Water Content (g/m3)"
LWC_plane = wrf_user_intrp3d( LWC,z,"v",plane,0.,opts)
LWC_plane = LWC_plane*1000.
opts_LWC = opts_clouds
end if
if (isvar("qs"))
qs_plane = wrf_user_intrp3d( qs,z,"v",plane,0.,opts)
qs_plane = qs_plane*1000.
opts_qs = opts_clouds
end if
if (isvar("qi"))
qi_plane = wrf_user_intrp3d( qi,z,"v",plane,0.,opts)
qi_plane = qi_plane*1000.
opts_qi = opts_clouds
end if
if (isvar("qg"))
qg_plane = wrf_user_intrp3d( qg,z,"v",plane,0.,opts)
qg_plane = qg_plane*1000.
opts_qg = opts_clouds
end if
tc_plane = wrf_user_intrp3d(tc,z,"v",plane,0.,opts)
dim = dimsizes(qv_plane) ; Find the data span - for use in labels
zspan = dim(0)
; Options for XY Plots
opts_xy = res
opts_xy@tiYAxisString = "Height (km)"
opts_xy@AspectRatio = 0.75
opts_xy@cnMissingValPerimOn = True
opts_xy@cnMissingValFillColor = 0
opts_xy@cnMissingValFillPattern = 11
opts_xy@tmYLMode = "Explicit"
opts_xy@tmYLValues = fspan(0,zspan,nz) ; Create tick marks
opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels
opts_xy@tiXAxisFontHeightF = 0.020
opts_xy@tiYAxisFontHeightF = 0.020
opts_xy@tmXBMajorLengthF = 0.02
opts_xy@tmYLMajorLengthF = 0.02
opts_xy@tmYLLabelFontHeightF = 0.015
opts_xy@PlotOrientation = tc_plane@Orientation
; Plotting options for RH
opts_clouds = opts_xy
opts_clouds@pmLabelBarOrthogonalPosF = -0.07
;opts_clouds@ContourParameters = (/ 10., 90., 10. /)
;opts_clouds@ContourParameters = (/ 1., 20., 1. /)
opts_w = opts_xy
opts_PWC = opts_xy
opts_LWC = opts_xy
opts_clouds@cnFillOn = True
opts_clouds@cnFillColors = (/"White","White","White", \
"White","Chartreuse","Green", \
"Green3","Green4", \
"ForestGreen","PaleGreen4"/)
; Plotting options for Temperature
opts_tc = opts_xy
opts_tc@cnInfoLabelOrthogonalPosF = 0.00
opts_tc@ContourParameters = (/ 5. /)
; Get the contour info for the rh and temp
contour_tc = wrf_contour(a,wks,tc_plane,opts_tc)
contour_w = wrf_contour(a,wks, w_plane,opts_w)
contour_qv = wrf_contour(a,wks,qv_plane,opts_clouds)
contour_qr = wrf_contour(a,wks, qr_plane,opts_clouds)
contour_qc = wrf_contour(a,wks, qc_plane,opts_clouds)
contour_qs = wrf_contour(a,wks, qs_plane,opts_clouds)
contour_qi = wrf_contour(a,wks, qi_plane,opts_clouds)
contour_qg = wrf_contour(a,wks, qg_plane,opts_clouds)
contour_PWC = wrf_contour(a,wks, PWC_plane,opts_PWC)
contour_LWC = wrf_contour(a,wks, LWC_plane,opts_LWC)
; MAKE PLOTS
plot = wrf_overlays(a,wks,(/contour_w,contour_tc/),pltres)
; Delete options and fields, so we don't have carry over
delete(opts_tc)
delete(opts_clouds)
delete(tc_plane)
delete(qv_plane)
delete(qc_plane)
delete(qr_plane)
delete(qs_plane)
delete(qi_plane)
delete(qg_plane)
delete(w_plane)
delete(PWC_plane)
delete(LWC_plane)
end do ; make next cross section
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
end do ; END OF TIME LOOP
end
> ********************************************************
> Cecille M. Villanueva Birriel
> Ph.D. Candidate
> Cloud Microphysics Research Group
> Purdue University
> Department of Earth, Atmospheric, & Planetary Sciences
> 550 Stadium Mall Drive, West Lafayette, IN 47907-2051
> email: cvillanu@purdue.edu
> ********************************************************
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Dec 17 12:16:54 2012
This archive was generated by hypermail 2.1.8 : Fri Dec 21 2012 - 10:43:23 MST