[ncl-talk] how to fill an array with an other

[Buzan, Jonathan]

Hi Rashed,

One way to speed up this process is to use the ncl ndtooned (n-dimension to one dimension) function, then loop over the 1d blocks by an amount equal to the spread in lats Lons, etc. When your calculations are completed you bring your data back to the original n-dimensions. It is complicated to code, and takes some effort, but you can substantially speed up your code this way.

Cheers,
-Jonathan




On May 6, 2020, at 8:52 PM, Dennis Shea via ncl-talk <ncl-talk at ucar.edu<mailto:ncl-talk at ucar.edu>> wrote:

FYI:

In general, looping in an interpreted language [NCL, Matlab, Python, R, ...] is 'slow' relative to a compiled language like fortran, etc. Sometimes, quite slow. To my knowledge, one 'key' to fortran's speed is that fortran does not allow unrestricted pointers.  This allows the compiler to do 'data look-ahead' [my words, I am not a computer science person]. The compiler can issue 'pre-fetch' instructions tp load (say) cache so when the data are needed it is in fast memory. Please correct me if I am wrong.

I would say that looping in NCL is reasonably 'quick.'  As example, see the code for the Monthly Water Balance Model described at:  http://www.ncl.ucar.edu/Applications/crop.shtml
Code:   http://www.ncl.ucar.edu/Applications/Scripts/MWBM_DRIVER.ncl

In the 3-dimension (time,lat,lon) case, the following loops aver all times and grid points

The following occurred surprisingly quickly

;========================================================
; Brute Force: Loop over each time step and all 'land' grid points
;========================================================

  ntim = dim_data(0)
[**SNIP**]
  elseif (rank_prc.eq.3) then
    nlat = dim_data(1)
    mlon = dim_data(2)

    remain  := new( (/nlat,mlon/), type_data, type_fill)
    remain   =  0.0                                   ; 25.4
    remain at long_name= "Remain Precip Water"
    remain at units    = "mm"

    prestor := new( (/nlat,mlon/), type_data, type_fill)
    prestor  =  0.0                                   ; 25.4
    prestor at long_name= "Remain Precip Water"
    prestor at units    = "mm"

    if (isscalar(lsmask) .and. lsmask.eq.1) then
        lsm := new( (/nlat,mlon/), "integer", "No_FillValue")
        lsm  = 1
    else
        lsm  = lsmask
    end if

    do nl=0,nlat-1
      do ml=0,mlon-1
        if (lsm(nl,ml).eq.1) then   ; land only
          do nt=0,ntim-1
             prain(nt,nl,ml) = prc(nt,nl,ml)-snow(nt,nl,ml)
             rodirect(nt,nl,ml)= prain(nt,nl,ml) *directfac           ; Direct Runoff:  Eq (3)
             prain(nt,nl,ml)   = prain(nt,nl,ml) -rodirect(nt,nl,ml)  ; remaining rain
             snstor(nt,nl,ml)  = snstor(nt,nl,ml)+snow(nt,nl,ml)      ; psnow accumulates as snstor (page 1)
             if (snstor(nt,nl,ml).gt.0 .and. tmp(nt,nl,ml).gt.TSNOW) then
                 if (snstor(nt,nl,ml).le.10) then
                     smelt(nt,nl,ml)  = snstor(nt,nl,ml)
                 else
                     smeltf(nt,nl,ml) = MELTMAX*((tmp(nt,nl,ml)-TSNOW)/((TRAIN-TSNOW)+0.0001))  ; Eq (5)
                     if (smeltf(nt,nl,ml).gt.MELTMAX) then
                         smeltf(nt,nl,ml) = MELTMAX
                     end if
                     smelt(nt,nl,ml)  = smeltf(nt,nl,ml)*snstor(nt,nl,ml)       ; Eq (6)
                     snstor(nt,nl,ml) = snstor(nt,nl,ml)-smelt(nt,nl,ml)
                     if (snstor(nt,nl,ml).lt.0.0) then
                         snstor(nt,nl,ml) = 0.0
                     end if
                end if
             end if
             prain(nt,nl,ml)   = prain(nt,nl,ml) + smelt(nt,nl,ml)
             pmpe(nt,nl,ml)    = prain(nt,nl,ml) - pet(nt,nl,ml)    ; net water

             if (pmpe(nt,nl,ml).lt.0) then
                 sms(nt,nl,ml) =  prestor(nl,ml) - abs((pmpe(nt,nl,ml)*prestor(nl,ml))/WHC)
                 if (sms(nt,nl,ml).lt.0)  then
                     sms(nt,nl,ml) = 0
                 end if
                 delstor = sms(nt,nl,ml)   - prestor(nl,ml)       ; change in soil moisture
                 aet(nt,nl,ml) = prain(nt,nl,ml) - delstor       ; JAVA: ae = prain + delstor * (-1.0);
                 prestor(nl,ml) = sms(nt,nl,ml)                   ; save for next iteration
                 surplus(nt,nl,ml) = 0.0
            else
                 aet(nt,nl,ml) = pet(nt,nl,ml)
                 sms(nt,nl,ml) = prestor(nl,ml) + pmpe(nt,nl,ml)
                 if (sms(nt,nl,ml).gt.WHC) then
                     surplus(nt,nl,ml) = sms(nt,nl,ml) - WHC
                     sms(nt,nl,ml)  = WHC
                 end if
                 prestor(nl,ml) = sms(nt,nl,ml)                  ; save for next iteration
            end if

            deficit(nt,nl,ml) = pet(nt,nl,ml) - aet(nt,nl,ml)
            runoff(nt,nl,ml) = (surplus(nt,nl,ml)+remain(nl,ml))*runoffFactor
            remain     = (surplus(nt,nl,ml)+remain(nl,ml))-runoff(nt,nl,ml)
            if (remain(nl,ml).lt.0) then
                remain(nl,ml) = 0
            end if

            runoff(nt,nl,ml)  = runoff(nt,nl,ml) + rodirect(nt,nl,ml)
          end do  ; nt
        end if    ; lsm
      end do      ; ml
     end do       ; nl

  end if          ; rank_data

_______________________________________________
ncl-talk mailing list
ncl-talk at ucar.edu<mailto:ncl-talk at ucar.edu>
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/20200506/bcf9003d/attachment.html>