[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