[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