Re: speed up maxind loop?

From: Seifert Axel <Axel.Seifert_at_nyahnyahspammersnyahnyah>
Date: Sat Dec 18 2010 - 01:54:22 MST

Will,

Thanks for your comment. But before writing my own dim_maxind_n I thought I could ask whether somebody else already stumbled across this issue. Measuring layer thicknesses seems not such an uncommon problem.

The reversal of the coordinate in minind(cwc(kwmax(j,i):0,j,i) is not a typo, because minind returns the index of the FIRST occurrence of a minimum value, i.e.

ncl 0> a = (/3,2,1,0,0,0/)
ncl 1> print(minind(a(0:)))
(0) 3
ncl 2> print(minind(a(:0)))
(0) 0

Nevertheless, thanks for trying to help.

I should have mentioned that I do mask the w array before the loop

; mask vertical velocity as conditional sampling of deep convection
; (1) w larger than 1 m/s
; (2) cloud water larger than 0.5 g/kg
; (3) CAPE larger than 10 J/kg
      
      w = mask(w,cwc.lt.0.5e-3.or.w.lt.1,False)
      w = mask(w,cape.lt.10,False)

Maybe ind_resolve could be used to collapse the double-loop over the whole array into a single-loop over the relevant columns? Thats at least what I would try on a vector machine. Would this help to speed up ncl, too?

Axel

-----Ursprüngliche Nachricht-----
Von: ncl-talk-bounces@ucar.edu im Auftrag von Hobbs, Will R (3244-CalTech)
Gesendet: Fr 17.12.2010 19:36
An: Axel Seifert; ncl-talk
Betreff: Re: [ncl-talk] speed up maxind loop?
 
Axel

You have a couple of options; the first being to write your own dim_maxind_n function!

The second would be to remove as many operations as possible from the loop; for example the statement ' zcdepth(j,i) = (/ zctop(j,i) - zcbase(j,i) /)' would be better after the loop has closed. There are a few places where you could do similar things (another example is subtracting 1e-6 from cwc before the loop starts; that's one operation as opposed to xdim*ydim operations).

I also notice an error in one of your lines that might be really slowing you down:

     ktop(j,i) = kwmax(j,i) - minind(cwc(kwmax(j,i):0,j,i) - 1e-6)

SHOULD read:

     ktop(j,i) = kwmax(j,i) - minind(cwc(0:kwmax(j,i),j,i) - 1e-6) ;order of the zero in indexing cwc has changed.

I 'speculate' that because the indexing has been inverted, NCL might be doing an unnecessary reversal at every time step

Will

On 12/17/10 6:57 AM, "Axel Seifert" <Axel.Seifert@dwd.de> wrote:

Hi NCL folks:

I am using maxind and minind to measure cloud depth in CRM data. This
works fine, but it is kind of slow. Unfortunately, maxind and minind do
not support multidimensional arrays. There is no dim_maxind_n available,
is there?

This is the loop that I want to speed up (I have some TB of data to analyze):

      print(" calculate convective cloud depth (please be patient, this can be slow)")
      do i=0,xdim-1
        do j=0,ydim-1

; find maximum of vertical velocity in the updraft
          kwmax(j,i) = maxind(w(:,j,i))

          if (.not.ismissing(kwmax(j,i))) then
; find cloud top and cloud base starting from max updraft position
            ktop(j,i) = kwmax(j,i) - minind(cwc(kwmax(j,i):0,j,i) - 1e-6)
            kbase(j,i) = kwmax(j,i) + minind(cwc(kwmax(j,i):,j,i) - 1e-6)

; store geometric height of cloud base, top and depth
            zctop(j,i) = (/ z(ktop(j,i),j,i) /)
            zcbase(j,i) = (/ z(kbase(j,i),j,i) /)
            zcdepth(j,i) = (/ zctop(j,i) - zcbase(j,i) /)
          end if

        end do
      end do

Cheers, Axel

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk

**********************************************************
Will Hobbs, Ph.D. William.R.Hobbs@jpl.nasa.gov
Jet Propulsion Laboratory
4800 Oak Grove Dr. office: 300-324a
M/S 300-323 phone: (818) 354-0466
Pasadena, CA 91109 fax: (818) 354-0966
**********************************************************

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sat Dec 18 01:55:20 2010

This archive was generated by hypermail 2.1.8 : Wed Dec 22 2010 - 16:10:23 MST