NOTES:

Tasking region causes each instance of each procedure call to spawn a
task.  Because U and V are inputs, there are no loop-carried
dependences, first ilayer loop can be unrolled
statically.

Note that istep loop cannot be unrolled until runtime, because of
dynamic loop bound. However, since nsteps is an input value, it
can be unrolled at the beginning of execution.

Also, nhrs is input value, so nhrs instances of unrolled ilayer loops
can be further unrolled before execution.

Manual transpose: use IDXPERM for fast transpose.


-----------------------------------------------------------------------------
parameter nlayers=5, nspecies=50, ncols=3500

real A(nlayers,nspecies,ncols)      ! concenteration vector 
real B(nspecies,ncols)              ! input emissions vector
real U(nlayers,ncols),V(nlayers,ncols)    ! wind velocities  
real K(nlayers,ncols,ncols)         ! global stiffness matrix for layer l 

integer nlayers                     ! Constant parameter

distribute A(BLOCK,BLOCK,:)
distribute U(BLOCK,:), V(BLOCK,:) , K(BLOCK,:,:)

Init(A,nlayers,nspecies, ncols)

begin tasking
do hour = 1, nhrs
   getemission(B)
   OUT(B)

   update(A,B)
   INOUT(A),IN(B)

   read(nsteps)
   OUT(nsteps)

   do ilayer = 1,nlayers            ! loop over layers to allow pipelining
     call getvelocities(U(ilayer),V(ilayer))
     OUT(U(ilayer),V(ilayer))

     call assemble(K(ilayer),U(ilayer),V(ilayer))
     IN(U(ilayer),V(ilayer)), OUT(K(ilayer))

     call factor(K(ilayer))
     INOUT(K(ilayer))

   enddo

   do istep = 1, nsteps

      call alltrans(A,K,U,V)
      INOUT(A), IN(K,U,V)

      call allchem(A)
      INOUT(A)

      call alltrans(A,K,U,V)
      INOUT(A), IN(K,U,V)

   enddo
enddo
end tasking


subroutine alltrans(A,K,U,V)
distribute A(BLOCK,BLOCK,:)
distribute U(BLOCK,:), V(BLOCK,:), K(BLOCK,:,:)

pdo (ilayer = 1:nlayers, ispec = 1:nspecies)
  INOUT(A(ilayer,ispec,:), IN(U(ilayer,:), V(ilayer,:), K(ilayer,:,:))
  call transport(A(ilayer,ispec), U(ilayer), V(ilayer), K(ilayer,:,:))
enddo

subroutine getemission(B)
<input B>
end


subroutine update(A,B)
distribute A(BLOCK,BLOCK,:), B(BLOCK,:)
   A(1,:,:) += B(:,:)
end

subroutine chemistry(A)
distribute A(:,:,BLOCK)


pdo icol = 1, ncols
   INOUT(A(:,:,icol)
   call chemistry(A(:,:,icol)  ! local computation 
end


subroutine factor(K)
distribute K(BLOCK,:,:)

   pdo ilayer = 1,nlayers

      factor(K(ilayer))       ! local computation 
      INOUT(K(ilayer,:,:)

   enddo

end
