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

Dennis Shea shea at ucar.edu
Wed May 6 12:52:51 MDT 2020


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*
<http://www.ncl.ucar.edu/Applications/crop.shtml>
Code:   *http://www.ncl.ucar.edu/Applications/Scripts/MWBM_DRIVER.ncl*
<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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20200506/fca1599f/attachment.html>


More information about the ncl-talk mailing list