[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