[ncl-talk] (no subject)
Modise Wiston
modise.wiston at postgrad.manchester.ac.uk
Sat Aug 2 12:06:36 MDT 2014
To whom it may concern
Hi all,
I am trying to determine the cloud droplet number concentration (QNDROP) from my [WRF-Chem] model simulation. The script [see the code below] is such that the code loops through the x-y-z (i-j-k) planes and finds the max (& min) QNDROP along the vertical (model levels). The code is such that it exits/stops just after the first cloud band only -without continuing through all the layers above the first cloud band. This also calculates CCN number concentrations below the cloud base.
How can I make it (or write a function/code) to calculate the height at which the max QNDROP occurs..? I have been trying to go through several NCL examples (functions) but couldn't figure out out how to do this.. (I would like to plot 'height vs max QNDROP' on the horizontal lat-lon domain)
Any help is appreciated please..
Thank you,
M. Wiston
The University of Manchester
-----------------------------------------------------------------------------------------------------------
dt = 6
do ifil = 0, 0; numFILES-48, dt
a = addfile(FILES(ifil)+".nc","r")
;list time step to output data
times = wrf_user_list_times(a)
ntimes = dimsizes(times)
do it = 0,ntimes-1,2
do ispec = 0, nspec-1
species = speclist(ispec)
;Obtain the variables from the file
if(isfilevar(a,"QNDROP"))
N = wrf_user_getvar(a,"QNDROP",it)
Nd = (N*1.2923)/1e6 ; convert from /kg to /sm^3 then to /scm^3
Nd at units = "/cm^3" ; converted units
end if
;define an array to hold data
dimsNd = dimsizes(Nd)
print(dimsNd)
;----------------Determine the domain size (grid points and vert. levels)----------------;
; (321 x 251) ; grid points along "i" (hor) by "j" (vert) direction ;
; (e_we:lat, e_sn:lon) ;
; 40 vert levels:(k) ; model levels ;
; Print output gives: k = 40 ;NB: this is read starting from zero (i.e. 0-39) ;
; j = 250 ;
; i = 320 ;
;------------------------------------------------------------------------------------------------;
min_Nd = new( (/dimsNd(1),dimsNd(2)/),float )
max_Nd = new( (/dimsNd(1),dimsNd(2)/),float )
c_btm = new( (/dimsNd(1),dimsNd(2)/),integer )
c_top = new( (/dimsNd(1),dimsNd(2)/),integer )
do j = 0, dimsNd(1)-1 ; loop along j
do i = 0, dimsNd(2)-1 ; loop along i
; cldfrc = (i*j)
c_btm(j,i) = -1 ;define cloud base and cloud top limits to confine the execution
c_top(j,i) = -1
do k = 0, dimsNd(0)-2 ; loop [vertically] through levels
; print("k= "+k)
if(Nd(k,j,i) .le. cut_pnt .and. Nd((k+1),j,i) .gt. cut_pnt)
c_btm(j,i) = (k+1) ;[cloud base]
; print("c_btm = "+c_btm(j,i))
else if (Nd(k,j,i) .gt. cut_pnt .and. Nd((k+1),j,i) .le. cut_pnt) then
c_top(j,i) = k ;[cloud top]
break ; exit the the process after the first cloud band
end if
end if
end do
; min_Nd(i) = min(Nd(:,i))
; max_Nd(i) = max(Nd(:,i))
; min_Nd(j) = min(Nd(:,j))
; max_Nd(j) = max(Nd(:,j))
; min_Nd(j,i) = min(Nd(:,j,i))
; max_Nd(j,i) = max(Nd(0:3,j,i))
if(c_btm(j,i) .eq. -1 .and. c_top(j,i) .eq. -1)
max_Nd(j,i) = 0.0
else if (c_btm(j,i) .gt. -1 .and. c_top(j,i) .gt. -1) then
max_Nd(j,i) = max( Nd(c_btm(j,i):c_top(j,i),j,i) )
;print(min_Nd(j,i))
;print("max_Nd = " +max_Nd(j,i))
print("| "+sprintf("%5.3f",max_Nd(j,i))+ " | "+sprintf("%2.1f", c_btm(j,i))+" | "+sprintf("%2.1f", c_btm(j,i))+" |" )
end if
end if
end do
end do
;------------------------------------------------------------------------------------------
;;get the CCN number concentration from the file
Cn1 = wrf_user_getvar(a,"CCN1",it)
Cn2 = wrf_user_getvar(a,"CCN2",it)
Cn3 = wrf_user_getvar(a,"CCN3",it)
Cn4 = wrf_user_getvar(a,"CCN4",it)
Cn5 = wrf_user_getvar(a,"CCN5",it)
Cn6 = wrf_user_getvar(a,"CCN6",it)
z = wrf_user_getvar(a,"z",it)
w = wrf_user_getvar(a,"wa",it) ; vertical [updraft..??] velocity
w at units = "m/s"
p = wrf_user_getvar(a, "p",it) ;extract p and tk to convert to STP
tk = wrf_user_getvar(a,"tk",it)
C1 = Cn1*((tk/p)*(P_std/T_std))
C2 = Cn2*((tk/p)*(P_std/T_std))
C3 = Cn3*((tk/p)*(P_std/T_std))
C4 = Cn4*((tk/p)*(P_std/T_std))
C5 = Cn5*((tk/p)*(P_std/T_std))
C6 = Cn6*((tk/p)*(P_std/T_std))
C1 at units = "/cm^3"
C1_btm = new( (/dimsNd(1),dimsNd(2)/),float )
C2_btm = new( (/dimsNd(1),dimsNd(2)/),float )
C3_btm = new( (/dimsNd(1),dimsNd(2)/),float )
C4_btm = new( (/dimsNd(1),dimsNd(2)/),float )
C5_btm = new( (/dimsNd(1),dimsNd(2)/),float )
C6_btm = new( (/dimsNd(1),dimsNd(2)/),float )
SS = new( (/dimsNd(1),dimsNd(2)/),float )
do j = 0, dimsNd(1)-1
do i = 0, dimsNd(2)-1
if(c_btm(j,i) .eq. -1)
C1_btm(j,i) = 0.0
C2_btm(j,i) = 0.0
C3_btm(j,i) = 0.0
C4_btm(j,i) = 0.0
C5_btm(j,i) = 0.0
C6_btm(j,i) = 0.0
else if(c_btm(j,i) .ne. -1) then
C1_btm(j,i) = ( C1(c_btm(j,i),j,i) )
C2_btm(j,i) = ( C2(c_btm(j,i),j,i) )
C3_btm(j,i) = ( C3(c_btm(j,i),j,i) )
C4_btm(j,i) = ( C4(c_btm(j,i),j,i) )
C5_btm(j,i) = ( C5(c_btm(j,i),j,i) )
C6_btm(j,i) = ( C6(c_btm(j,i),j,i) )
; print("Cn1 = " +C1_btm(j,i))
; print("Cn2 = " +C2_btm(j,i))
; print("Cn3 = " +C3_btm(j,i))
; print("Cn4 = " +C4_btm(j,i))
; print("Cn5 = " +C5_btm(j,i))
; print("Cn6 = " +C6_btm(j,i))
print("| "+sprintf("%4.3f", C1_btm(j,i))+" | "+sprintf("%4.3f", C2_btm(j,i))+" | "+sprintf("%4.3f", C3_btm(j,i))+" | "+\
sprintf("%4.3f", C4_btm(j,i))+" | "+sprintf("%4.3f", C5_btm(j,i))+" | "+sprintf("%4.3f", C6_btm(j,i))+" |" )
end if
end if
end do
end do
;`````````````````````````````````````````````````````````````````````````````````
; if(max_Nd(j,i) .le. C1_btm(j,i))
; SS(j,i) = 0.01
; else if(max_Nd(j,i) .gt. C1_btm(j,i) .and. max_Nd(j,i) .le. C2_btm(j,i))
; SS(j,i) = "0.02%"
; else if(max_Nd(j,i) .gt. C2_btm(j,i) .and. max_Nd(j,i) .le. C3_btm(j,i))
; SS(j,i) = "0.1%"
; else if (max_Nd(j,i) .gt. C3_btm(j,i) .and. max_Nd(j,i) .le. C4_btm(j,i))
; SS(j,i) = "0.2%"
; else if (max_Nd(j,i) .gt. C4_btm(j,i) .and. max_Nd(j,i) .le. C5_btm(j,i))
; SS(j,i) = "0.5%"
; else if (max_Nd(j,i) .gt. C5_btm(j,i) .and. max_Nd(j,i) .le. C6_btm(j,i))
; SS(j,i) = "1.0%"
; end if
; end if
; end if
; end if
; end if
; end if
;-------------------------------------------------------------------------------------------
;Set some Basic Plot options
res = True
res at tiMainString = "CCN nmbr @: " +times(it) ; main field title
mpres = True
pltres = True
res at gsnMaximize = False
res at gsnFrame = False
res at gsnDraw = False
res at cnFillOn = True
res at cnLinesOn = False ; turn off contour lines
res at cnLineLabelsOn = False ; lables the data on contours/map grids
res at gsnSpreadColors = True
res at lbOrientation = "Horizontal" ; orientates the linear scale
res at cnLevelSelectionMode = "ExplicitLevels"
;-------------------------------------------------------------------------------------------
res at tiYAxisString = ""
res at tiXAxisString = "Nmbr cont'n (/cm^3)"
res at tiXAxisSide = "Bottom"
;res at tiMainString = "Droplet nmbr at: " +times(it)
contour_Nd = wrf_contour(a,wks,max_Nd,res)
plot_Nd = wrf_map_overlays(a,wks,(/contour_Nd/),pltres,res) ; plots data on a map backgrnd
; contour_Cn1 = wrf_contour(a,wks,C1_btm,res)
; plot_Cn1 = wrf_map_overlays(a,wks,(/contour_Cn1/),pltres,res)
; contour_Cn2 = wrf_contour(a,wks,C2_btm,res)
; plot_Cn2 = wrf_map_overlays(a,wks,(/contour_Cn2/),pltres,res)
; contour_Cn3 = wrf_contour(a,wks,C3_btm,res)
; plot_Cn3 = wrf_map_overlays(a,wks,(/contour_Cn3/),pltres,res)
; contour_Cn4 = wrf_contour(a,wks,C4_btm,res)
; plot_Cn4 = wrf_map_overlays(a,wks,(/contour_Cn4/),pltres,res)
; contour_Cn5 = wrf_contour(a,wks,C5_btm,res)
; plot_Cn5 = wrf_map_overlays(a,wks,(/contour_Cn5/),pltres,res)
; contour_Cn6 = wrf_contour(a,wks,C6_btm,res)
; plot_Cn6 = wrf_map_overlays(a,wks,(/contour_Cn6/),pltres,res)
; delete(contour_Nd)
; delete(contour_Cn1)
; map = wrf_map(wks,a,res)
; overlay(map,contour_Nd)
; overlay(map,contour_Cn1)
; overlay(plot0,plot3)
; gsn_panel(wks,plots,(/3,2/),True)
; draw(plot)
; frame(wks)
end do
end do
end
---------------------------------------------------------------------------------------------------------------
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20140802/a61c7651/attachment.html
More information about the ncl-talk
mailing list