Tom Micevski wrote:

here is my promised test program. i've just stripped down an old
program (but you could strip even more away, eg. data statements and
write statements). hope it's helpful.

in case the code gets wrapped by my news program, i've also uploaded it
here:
http://www.mytempdir.com/1176091

tom

program forall_test_flike
implicit none
integer, parameter :: rk = selected_real_kind(p=10,r=100) ! double precision
! test data (0.0 = missing data)
real(rk), parameter :: q1(30) = (/
34.1,45.1,38.5,60.9,75.4,71.6,365.7,93.1,49.0,98.2,516.9,&

51.6,361.1,75.5,0.0,26.7,30.4,29.7,49.0,45.9,50.0,18.5,73.4,26.1,92.4,71.0,48.0,0.0,1216.3,0.0
/)
real(rk), parameter :: q2(30) = (/
169.2,469.9,305.7,566.1,317.0,492.5,150.0,795.4,764.3,1372.9,&

3366.0,62.2,418.9,261.8,158.5,345.3,8516.7,206.6,71.3,42.5,210.9,0.0,0.0,0.0,0.0,0.0,670.3,712.6,178.6,748.3
/)
real(rk), parameter :: q3(30) = (/
353.8,0.0,2135.5,0.0,3677.0,2123.0,1631.0,3693.7,1826.6,1852.4,&

146.6,433.4,614.3,393.5,54.8,1073.4,1368.6,1492.5,155.6,361.5,1616.8,182.8,259.3,570.6,1302.4,1301.8,0.0,0.0,0.0,0.0
/)
real(rk), parameter :: q4(30) = (/
38.4,12.7,12.4,39.1,175.3,0.0,0.0,0.0,0.0,0.0,&

0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,101.1,16.2,59.5,366.4,0.0
/)
real(rk), parameter :: q5(30) = (/ 0.0,0.0,0.0,0.0,0.0,0.3,0.6,1.2,0.9,4.5,&

4.3,5.6,1.5,1.6,2.4,3.0,1.6,4.4,10.1,2.1,3.9,4.4,2.4,0.7,3.1,3.2,1.6,1.8,5.0,3.2
/)
real(rk), parameter :: q6(30) = (/
1372.9,279.9,202.4,40416.9,7522.6,2533.5,3304.0,1231.3,1390.4,12536.6,&
3639.8,448.7,478.4,180.2,164.2,229.3,2123.0,965.3,27416.5,49.0,&
76.4,571.5,925.6,468.8,348.2,5412.3,3632.3,942.6,1006.5,335.2 /)
real(rk), parameter :: q7(30) = (/ 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,&

0.0,22.9,53.4,44.5,14.4,36.6,216.5,34.0,25.0,33.0,8.5,10.7,25.1,0.0,0.0,4.5,44.0,0.0,0.0,0.0
/)
real(rk) :: q(30,7) = reshape( (/q1,q2,q3,q4,q5,q6,q7/), (/30,7/) )
! vars
integer :: npos_min = 15
real(rk) :: zero = 0.0_rk, miss = -999.0_rk
integer :: nsite, nyear, i, j
integer, allocatable :: npos(:)
real(rk), allocatable :: mn(:), sd(:)
character*20, allocatable :: text(:)
! start program
nyear = size(q, dim=1)
nsite = size(q, dim=2)
! mark/flag missing data (2 versions [1 using forall])
where (q <= zero .and. q > miss) q = miss - 10.0
forall (i=1:nsite, j=1:nyear, q(j,i) <= zero .and. q(j,i) > miss)
q(j,i) = miss - 10.0
end forall
do i=1,nyear
write(*,'(999f8.1)') (real(q(i,j)), j=1,nsite)
enddo
! count no. of positive values (2 versions [1 using forall])
allocate (npos(nsite))
npos = count(q > zero, dim=1)
write(*,*) npos
forall (i=1:nsite) npos(i) = count(q(:,i) > miss)
write(*,*) npos
! calc mead/sd (outer forall [with mask] calling pure subroutines [1
with inner forall])
allocate (mn(nsite), sd(nsite))
mn = zero; sd = zero
forall (i=1:nsite, npos(i) >= npos_min)
mn(i) = getmean (x=pack(q(:,i), q(:,i) > miss))
sd(i) = getsd (x=pack(q(:,i), q(:,i) > miss), m=mn(i))
end forall
write(*,*) real(mn)
write(*,*) real(sd)
! now just add a string-based forall
allocate (text(nsite))
forall (i=1:nsite) text(i) = 'site='//number_string(i)
write(*,*) text
continue
stop
contains
! ======================================================
pure function getmean (x)
IMPLICIT NONE
REAL(rk), INTENT(IN) :: x(:)
REAL(rk) :: getmean
real(rk) :: nr
nr = real(size(x),rk)
getmean = sum(x)/nr
end function getmean
! ======================================================
pure function getsd (x, m)
IMPLICIT NONE
REAL(rk), INTENT(IN) :: x(:), m
REAL(rk) :: getsd
INTEGER :: i
real(rk) :: nr1, resid(size(x))
nr1 = real(size(x)+1,rk)
forall (i=1:size(x)) resid(i) = x(i) - m
getsd = dot_product(resid,resid)/nr1
getsd = sqrt(getsd)
end function getsd
! ======================================================
pure function number_string (int)
implicit none
integer, intent(in) :: int
character(len=20) :: number_string
write(number_string,'(i0)') int
END function number_string

end program forall_test_flike