[ncl-talk] (no subject)
Dennis Shea
shea at ucar.edu
Mon Aug 4 08:44:26 MDT 2014
While not WRF, the last example at:
http://www.ncl.ucar.edu/Applications/iso.shtml
may give you an idea of how to approach the issue.
===
Also, please send *all* WRF questions to wrfhelp at ucar.edu
The know the WRF model best.
You can cc ncl-talk at ucar but, really, wrfhelp should be your primary
contact for assorted WRF related issues.
Good luck
On Sat, Aug 2, 2014 at 12:06 PM, Modise Wiston <
modise.wiston at postgrad.manchester.ac.uk> wrote:
> *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
>
> ---------------------------------------------------------------------------------------------------------------
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20140804/be554295/attachment.html
More information about the ncl-talk
mailing list