\documentstyle [12pt, sane]{article}
\setlength{\textheight}{8.7in}
\setlength{\topmargin}{-1cm}
\setlength{\oddsidemargin}{+7mm}
\setlength{\textwidth}{6.2in}
\pagestyle{fancy}
\input{psfig}

\begin{document}
\title{\bf{SCCS 694} \\
       Case Studies on Parallelization of NASA's Weather Code:\\
       Results on Interpolation Techniques}
\author{{\it Director} : Geoffrey Fox \\
{\it Project Leader} : Miloje Makivic \\
{\it Research Assistant} : Manu Natarajan \\ \\
Northeast Parallel Architectures Center \\
	Syracuse University, 111, College Place \\
	Syracuse, NY, 13244-4100 }


\date{}
\maketitle
\clearpage
\tableofcontents

\newpage
\noindent
\section{Abstract}

The purpose of this project is to investigate suitable domain decompositions
that lead to efficient parallelizations of the interpolation techniques
used in the Semi-Lagrangian advection scheme.  In this report we present
the results of our investigation using the lasagna decomposition of the
input domain of the bilinear and bicubic interpolations used in NASA's 
2-D Shallow Water Model.  

Our approach to the usage of the lasagna decomposition reduces the amount
of communication needed for either the bilinear or bicubic interpolations.
The criteria for goodness of the algorithm is based on the ratio of the
useful work done to the amount of communication required.

\section{Introduction}

Weather forecasts are essential for the planning of many activities. They need
to be accurate and available in time.  The basis of weather forecasting are
Numerical Weather Prediction (NWP) models which assimilate data and 
predict future conditions of atmosphere. Given the requirement that forecasts
be available in real-time, the numerical models must be efficient.  There is
also considerable interest in building quality data-sets or model predictions
of weather conditions for research purposes and as case studies.  These involve
processing of voluminous data using the numerical models in an efficient
manner.

One chronic problem with NWP models is that the longest time step that can be
chosen is decided by stability rather than accuracy.   This is alleviated to
a considerable degree in the Semi-Lagrangian Transport Model using a 
semi-implicit time differencing scheme \cite{rob-81,rob-82}. This was 
demonstrated for a limited area shallow water model (one layer, homogeneous
 fluid).  An
excellent review of Semi-Lagrangian integration schemes for atmospheric models 
is treated by \cite{stan-91}.  

Later work addressed the deficiencies associated with polar singularity 
in the latitude-longitude co-ordinate system \cite{mc-89}.  NASA
implemented a modified model based on this paper 

Another notable issue is that for accuracy to be improved, among other
things higher resolutions are required. This calls for efficient
implementations on powerful machines.  Effort has been underway to implement 
these models on parallel computers \footnote{Weather Modelling project 
undertaken at NASA Goddard Laboratory, Maryland}. 


In order to establish the most suitable paradigm, certain modules of NASA's
serial implementation were identified as test-cases. 
Initial studies at Syracuse University 
and NASA have indicated that the interpolation techniques deployed occupied a 
total running time of 25\% of running time in the serial code.  It was also 
observed that the usage of the construct "forall" did not have
a good speed-up in the data parallel prototype on the CM-5.  This can be 
attributed to the irregular, non-local pattern of
data accesses.  In addition the decomposition was carried out in the 
North-South axis as shown in Fig. 1; as one expects winds to be mostly in the 
East-West direction.
This lead to a very poor utilization of the nodes when there was 
a small grid size and  one of the grid axes 
declared as serial or local to a
node.  In such a case the data is distributed among only a few nodes along 
the news 
or parallel axis leaving remaining processors with no work \cite{cm5-ref}
The default news:news (square or ravioli) decomposition 
showed better results.  This is because the run time library is able to 
ensure a better utilzation of nodes with the default decomposition. But 
this decompostion does not take advantage the fact that the winds are of a
East-West nature.


Hence, the message passing model was proposed as an alternate
parallelization technique
as it offers a greater flexibility in choice of domain decomposition.  One 
does not need to depend on the vagaries of run-time libraries or systems
to decide the utilization of nodes.  Scalable decompositions may be 
developed.  Additionally, the decomposition can be adapted to ensure that 
the pole regions are treated differently than regions near the equator.  
This is because, at the poles the winds move around the pole (entire latitude).
This is another disadvantage using data parallel 
techniques which permit only simple decompositions.  Hence, nearest neighbor
communication is not possible at the pole for the default square decomposition
on CM-5 when using a latitude-longitude grid.

A lasagna decomposition along the N-S axis was adopted as the 
first test case as seen in Fig. 2.  This was postulated to have mostly 
nearest neighbor 
communication in the N-S direction.  The polar region problem is
avoided as the lasagna strips are in E-W direction. The algorithm 
uses a message-passing version of the data-parallel construct "cshift".
The maximum number of "cshifts" is determined at run-time at every time
-step for each component.  The necessary data from the circulated 
wind "strips" is collected and a local interpolation is performed at
every time step.  Two strategies and their merits are discussed in the
next section.

\subsection{Metrics}

We propose some measurements for assessing the goodness of the algorithm
as well as giving direction to future work :

\begin{itemize}
 \item Overall running time
 \item Ratio of calculation to communication.
 \item The maximum displacement in terms of nodal distance in the N-S
       direction.
 \item Speed-up of interpolation
 \item Accuracy of calculations 
\end{itemize}

The first criterion shows how well the algorithm performs with respect to 
the serial code as the problem size grows.  The numerical value is of not
much interest.  But, a very low figure would indicate that the algorithm
was ill-conceived.  The second may give a clue as to what could
be a suitable ratio of the partition dimensions in more complex decompositions.
It would also verify if the assumption that the winds are mostly E-W is true.
The third criterion indicates how well the algorithm will scale when a global
 3-D Semi-lagrangian model is used.  As far as the fourth criterion is 
concerned, speedups will be obtained by comparing the performance of the
code on a 32 node CM-5 with the performance on a single node on the CM-5.
The last criterion is important as it indicates the presence of errors other
than roundoff.  This done by finding out if the difference between the 
interpolated values obtained from both the serial and parallel implementations
are homogenous for a given resolution.

\section{Design}

The problem of parallelizing the bilinear and bicubic interpolation is first
considered by using a simple approach.  In this lasagna decomposition each node receives an input partition of the 
latitude, longitude grid and each of the two wind components [Fig 2].  Each
node then computes a set of points (ordered pairs) {(ii,jj)} from the
latitude, longitude grid, based on which an interpolation is carried out.
In case of the bilinear interpolation, associated with the set of points 
{(ii,jj)} there exist four sets of points {(ii,jj+2)}, {(ii,jj+3)}, 
{(ii+1),(jj+2)}, {(ii+1,jj+3)} which are required for the interpolation 
[Fig. 3].
In case of the bicubic interpolation there are sixteen such sets of points
[Fig. 4]
\footnote{ Hereafter,
the discussion will consider only the bilinear interpolation, but the same
reasoning holds good for the bicubic interpolation. Appendix A covers
the bicubic interpolation algorithm. }.

In the first approach to interpolation, the algorithm is outlined below.

\begin{itemize}
\item {Do for all (ii,jj)}
\begin{itemize}	
\item	Find the 4 points (ii,jj+2), (ii,jj+3), (ii+1,jj+2), (ii+1,jj+3)
\item   Obtain the values of the wind component at above points.
\begin{itemize}
	\item If the value is not available locally, {\em fetch} from
         the appropriate neighboring node
\end{itemize}
\item  Find the interpolation point.
\item  Compute the value at interpolation point.
\end{itemize}
\item End of do loop
\end{itemize}

It is obvious that though the above approach is simple, it is impractical
to implement because of the huge communication overhead involved.  For, every
point that is not available locally, a communication call will have to be
generated to fetch the value at that point.  There are at least two approaches
which lower the communication overhead. In both the
approaches strips of the wind components are shifted from node to node.  The
shifting occurs in both the southern and northern directions.  The amount of
shifting is determined by the approach discussed in the next section.

\subsection{Pre-Determination of Degree of Circulation}

If the assumption that the winds in the North-South direction mostly result in
nearest neighbor communicationis valid , then the following strategy may be 
used.
For any given global indice (ii, jj), we determine the node it belongs to.
The inter-nodal distance is then simply the difference between the owner node
and user node.  But, this distance has a direction, i.e. the ``owner'' node 
may be
reached in the North or South Direction [Fig. 5].  Hence, this is converted to
 the
inter-nodal distance between the two nodes when going in the North direction.
Moreover, if the inter-nodal distance exceeds half the number of nodes, the
shorter path is in the Southern direction.  For, either interpolation there are
two extreme points in the North-South direction [Fig. 3,4]. The inter-nodal 
distance is computed for both the points separately, for all the indices on a
node.  The maximum inter-nodal distance in a direction (both North and South
are evaluated) is then found by a global reduction. The two maxima in the 
North-South direction govern the circulation.

The algorithm for {\em the pre-determination of degree of circulation} is below
\begin{itemize}
  \item { ny = Number of nodes (processors) used }
  \item { Let 'myaddr' denote the address of the node executing the code}
  \item {Let 'ownp0' denote nodal address where value 'jj(i,j)+0' resides }
  \item {Let 'ownp3' denote nodal address where value 'jj(i,j)+3' resides }
  \item { Number of elements in E-W direction for a partition is imp3 }
  \item { Number of elements in N-S direction for a partition is range }
  \item {for j = 1, range}
  \begin {itemize}
    \item {for i = 1, imp3}
      \begin {itemize}
        \item {ownp0 = jj(i,j) / partition width }
        \item {ownp3 = (jj(i,j) + 3) / partition width }
        \item { Find the inter-nodal distance in {\em northern direction}
                from 'myaddr' to both ownp0 and ownp3}
	\begin{itemize}
          \item {distp0 = ownp0 - myaddr}
          \item {distp3 = ownp3 - myaddr}
          \item {if (distp0 .lt. 0),distp0 = ny + distp0}
          \item {if (distp3 .lt. 0),distp3 = ny + distp3}
        \end{itemize}

	\item {Find the shortest inter-nodal distance for both ownp0, ownp3, 
               and update the local maximum inter-nodal distances if required}
      \end{itemize}
    \end{itemize}

    \item {Find the global max. northern dist. (reduction operation)}
    \item {Find the global max. southern dist. (reduction operation)}
\end {itemize}


Expansion of {Updation of local maximum inter-nodal distances}
        \begin {itemize}
          \item {Let distp0 be the inter-nodal distance in 
                 {\em northern direction} from 'myaddr' to ownp0}
          \item {Let distp3 be the inter-nodal distance in 
                 {\em northern direction} from 'myaddr' to ownp3}
	  \item { If distp0 .gt. (ny / 2)}
            \begin{itemize}
              \item {distp0 = ny - distp0}
	      \item {if (max. southern dist .lt. distp0) max. southern dist =
                     distp0}
            \end{itemize}
          \item {Else}
          \begin{itemize}
            \item{if (max. northern dist .lt. distp0) max. northern dist =
                     distp0}
          \end{itemize}

	  \item { If distp3 .gt. (ny / 2)}
            \begin{itemize}
              \item {distp3 = ny - distp3}
	      \item {if (max. southern dist .lt. distp3) max. southern dist =
                     distp3}
            \end{itemize}
          \item {Else}
          \begin{itemize}
            \item{if (max. northern dist .lt. distp3) max. northern dist =
                     distp3}
          \end{itemize}
        \end{itemize}

\subsection { Reconstruction of the Sets of Points Approach }

In the second approach we wish to fully reconstruct the four associated sets 
of points {(ii,jj+2)}, {(ii,jj+3)}, {(ii+1),(jj+2)}, {(ii+1,jj+3)} in four
arrays  for every node.  Then the four arrays are indexed to obtain the four 
points required for interpolation on every node.  As a wind component has
been partitioned as shown in Fig. 2, there is an upper and a lower bound
on the index jj in the North-South direction that decides which values of
a wind component may be found locally.  Thus, at every shift of the wind
strips, the four arrays maybe partially reconstructed using the new bounds 
obtained with every shift.  At the end of the shifting the arrays should be
 fully reconstructed.  However, there will be points such that the points
 (ii,jj) will be local but (x,jj+2) or (x,jj+3) may not be local as seen in
 figure 3 \footnote
{ x maybe ii or ii + 1.  This does not matter as in the East-West direction
 the elements of the wind component are local to a node. }. Thus it will not 
 be possible to reconstruct {\em entirely locally } the four associated sets of 
 points for all jj that that fall within the bounds of the node's partition.
 {\em If the reconstruction is not done entirely locally, there would be additional
 communication overhead to fetch the missing values one at a time, as in the
 first approach. }  This hitch may be
 rectified by initially creating a 3-element wide strip in the Northern 
 direction by fetching those elements from the neighboring nodes [Fig. 6]. 
 The 
 creation of the Northern buffer is created before performing any shifts.
 The original partition and this northern buffer is now termed a segment.
 These segments are shifted instead of the original partition. The algorithm
 for the second approach is outlined as follows, which is executed once for
 each of the two wind components. \newline

{\em First level refinement} for the algorithm.

 \begin{itemize}
   \item { Number of elements in E-W direction for a partition is imp3 }
   \item { Number of elements in N-S direction for a partition is range }
   \item { Pre-determine the degree of circulation }
   \item { Construct the northern buffer. }
   \item {Select elements from segment}
   \item {Fetch segments and select elements from southern direction}
   \item {Fetch segments and select elements from northern direction}
   \item {Perform local interpolation}
  \end{itemize}

In this approach, one is guaranteed that if the value at jj(i,j) is resident
on a node, then the value at jj(i,j)+3 is resident \footnote{As there is the
3-element northern buffer appended to the original partition.}.  Therefore,
only the node where jj(i,j) is resident, need be considered for determination
of the local maximum inter-nodal distance.  This the only deviation from the
strategy outlined in the previous subsection. \newline

Algorithm for {\em construction of northern buffer}

 \begin{itemize}
   \item {Copy original wind partition to temporary buffer (depart)}
   \item {If partition width less than 3 then}
   \begin {itemize}
   \item {While northern buffer incomplete}
   \begin{itemize}
     \item { Send depart to southern node }
     \item { Receive partition into northern buffer from northern node}
     \item { Copy received partition into depart }
   \end{itemize}
   \item {End of while}
   \end{itemize}
   \item {Else (partition width is .ge. 3) }
   \begin{itemize}
     \item { Send original partition to southern node }
     \item { Receive new partition in northern buffer from northern node }
   \end{itemize}
   \item {Append northern buffer to original wind partition. This will
            be called segment. }
 \end{itemize}

\pagebreak

Algorithm for {\em fetching of segments from southern direction}

     \begin {itemize}
       \item {Current partition bounds = Original partition bounds }
       \item {Copy original segment into temporary buffer (depart)}
       \item {For number of hops in southern direction }
       \begin {itemize} 
	 \item {Send upper \& lower bounds of current partition to northern
                node.}     
         \item {Receive upper \& lower bounds of new partition from southern 
                node.}
         \item {Send 'depart' to northern node.}
         \item {Receive new segment from southern node in buffer 'landing'.}
	 \item {Current partition bounds = new partition bounds}
         \item {Copy 'landing' into 'depart'}
         \item {Select elements from landing (new segment)}
       \end {itemize}
      \end {itemize}


Algorithm for {\em fetching of segments from northern direction}

     \begin {itemize}
       \item {Current partition bounds = Original partition bounds }
       \item {Copy original segment into temporary buffer (depart)}
       \item {For number of hops in northern direction }
       \begin {itemize} 
	 \item {Send upper \& lower bounds of current partition to southern
                node.}     
         \item {Receive upper \& lower bounds of new partition from northern 
                node.}
         \item {Send 'depart' to southern node.}
         \item {Receive new segment from northern node in buffer 'landing'.}
	 \item {Current partition bounds = new partition bounds}
         \item {Copy 'landing' into 'depart'}
         \item {Select elements from landing (new segment)}
       \end {itemize}
      \end {itemize}

\pagebreak

Algorithm for {\em selection of elements from segment }
  \begin{itemize}
   \item { Number of elements in E-W direction for a partition is imp3 }
   \item { Number of elements in N-S direction for a partition is range }
   \item { Input : }
   \begin{itemize}
     \item {lower \& upper partition  bounds [Fig. 6]}
     \item {segment (corresponding to partition as in Fig 6)}
     \item {set of points \{(ii,jj)\}}
     \item {arrays to hold the four sets of points - w1, w2, w3, w4}
   \end{itemize}	
   \item {For i = 1, imp3}
     \begin{itemize}
       \item {For j = 1, range}
       \begin{itemize}
         \item {If (jj(i,j) .le. upper) .and. (jj(i,j) .ge. lower)}
         \begin{itemize}
	   \item {w1 = segment(ii(i,j), jj(i,j) +2)}
           \item {w2 = segment(ii(i,j), jj(i,j) +3)}
	   \item {w3 = segment(ii(i,j)+1, jj(i,j) +2)}
           \item {w4 = segment(ii(i,j)+1, jj(i,j) +3)}
         \end{itemize}
       \end{itemize}
%       \item {end of for i = ...}
      \end{itemize}
%      \item {end of for j = ...}
    \end{itemize}
   

Algorithm for {\em bilinear interpolation}.  Fig. 3 explains the 
terms delx, dely.

\begin{itemize}
   \item { Number of elements in E-W direction for a partition is imp3 }
   \item { Number of elements in N-S direction for a partition is range }
   \item {for j = 1, range}
   \begin{itemize}
     \item {for i = 1, imp3}
	\begin{itemize}
	  \item {value at interpolated point = 
             (1.0 - delx(i,j)) * 
                ((1.0 - dely(i,j)) * w1(i,j) + dely(i,j) * w2(i,j)) 
             + delx(i,j) * ((1.0 - dely(i,j)) * w3(i,j) + dely(i,j) * w4(i,j))}
        \end{itemize}
   \end{itemize}
\end{itemize}

The drawback with this approach is that the communication overhead for 
creating the northern buffer increases when the width of an individual strip 
is smaller than the overlap region [Fig. 7]. 

\subsection{Reconstruction of Segment Approach}

The third approach recognizes the fact that given the set of computed points
{(ii,jj)} on each node, the node only needs to see a segment of the entire
wind component to perform a local interpolation; there is no necessity for 
reconstructing the associated four sets of points.  This segment is not the
same as in the second approach, in that its size is not a static quantity. The
size of the segment depends on how far back one needs to go along
the fluid trajectory to perform an interpolation, and hence is a dynamic 
quantity.  Thus, no extra communication is needed to establish a northern
buffer. Additionally, only the {\em original } partitions of the wind 
component need be shifted to 
create the segment.  The circulated wind strips are appended as shown in 
Fig. 8. Some, book-keeping needs to be done in order to the relative position
of the partition in the segment and for generation of local indices.  
The interpolation is then carried out locally. The four points are obtained by
 the use of the local indices.  The algorithm is outlined below. \newline

The {\em first-level refinement} for the third approach. The individual 
steps are elaborated subsequently.

\begin{itemize}
   \item { Number of elements in E-W direction for a partition is imp3 }
   \item { Number of elements in N-S direction for a partition is range }
   \item { Pre-determine the degree of circulation }
   \item { Copy the partition into a large buffer called segment}
   \item { Fetch partitions from the southern direction and append to
	   segment. }
   \item { Fetch partitions from the northern direction and append to
	   segment. }
   \item { Convert global indices to local indices. }
   \item { Perform the interpolation locally. }
\end{itemize}

The pre-determination of the degree of circulation is exactly as outlined in 
section 3.1. 

\pagebreak

The algorithm for {\em fetching of partitions from the southern direction and
appending to the segment} is as below.

\begin {itemize}
       \item {Current partition bounds = Original partition bounds }
       \item {Copy partition into temporary buffer (depart)}
       \item {For number of hops in southern direction }
       \begin {itemize} 
	 \item {Send upper \& lower bounds of current partition to northern
                node.}     
         \item {Receive upper \& lower bounds of new partition from southern 
                node.}
         \item {Send 'depart' to northern node.}
         \item {Receive new partition from southern node in buffer 'landing'.}
         \item {Update parameters needed for generation of local indices.}
	 \item {Move up the previously filled columns in the segment by width
                of the new partition (southern partitions occupy lower 
                indices in the segments [Fig. 8]). }
         \item {Append the new partition to the 'southern end' of the segment.}
         \item {Update width of segment.}
         \item {Copy 'landing' into 'depart'}
       \end {itemize}
      \end {itemize}


The algorithm for {\em fetching of partitions from the northern direction and
appending to the segment} is as below.

\begin {itemize}
       \item {Current partition bounds = Original partition bounds }
       \item {Copy partition into temporary buffer (depart)}
       \item {For number of hops in northern direction }
       \begin {itemize} 
	 \item {Send upper \& lower bounds of current partition to southern
                node.}     
         \item {Receive upper \& lower bounds of new partition from northern 
                node.}
         \item {Send 'depart' to southern node.}
         \item {Receive new partition from northern node in buffer 'landing'.}
         \item {Update parameters needed for generation of local indices.}
         \item {Append 'landing' to the 'northern end' of the segment; 
                northern partitions occupy higher indices in the segment
                 [Fig 8].}
         \item {Update width of segment.}
         \item {Copy 'landing' into 'depart'}
       \end {itemize}
      \end {itemize}

\pagebreak

The local indices are generated with the following considerations :

\begin {itemize}
\item { The first column in any segment does not necessarily correspond
        to column '1' in the wind component. }
\item { There maybe a wrap around in the southern direction; i.e. the southern
most node communicated may cross the south pole. }

\item { Similarly there may be a wrap around in the northern direction. }

\end {itemize}

The algorithm {\em generation of local indices} is as below.  Fig. 9
 illustrates the algorithm.

\begin {itemize}
  \item {co\_lower = lower bound of non-wrapped part of segment }
  \item {wr\_lower = lower bound of wrapped part of segment }
  \item {wrapd = width of wrapped part of segment.}
  \item {Let 'njj' denote the set of new (local) indices}
  \item {If there is no wrap around in either direction}
    \begin{itemize}      
      \item {for i = 1, imp3}
	\begin {itemize}
	  \item {for j = 1, range}
	    \begin{itemize}      
		\item {njj(i,j) = jj(i,j) - co\_lower + 1}
            \end{itemize}
         \end{itemize}
     \end{itemize}

  \item {If there is wrap around in southern direction}
    \begin{itemize}      
      \item {for i = 1, imp3}
	\begin {itemize}
	  \item {for j = 1, range}
	    \begin{itemize}      
                \item {if jj(i,j) .ge. wr\_lower
                       then njj(i,j) = jj(i,j) - wr\_lower + 1}
		\item {Else njj(i,j) = jj(i,j) + wrapd}
            \end{itemize}
         \end{itemize}
     \end{itemize}

  \item {If there is wrap around in northern direction}
    \begin{itemize}      
      \item {for i = 1, imp3}
	\begin {itemize}
	  \item {for j = 1, range}
	    \begin{itemize}      
                \item {if jj(i,j) refers to the wrapped part
                       then njj(i,j) = jj(i,j) + jmp3 - co\_lower + 1}
		\item {Else njj(i,j) = jj(i,j) - co\_lower + 1}
            \end{itemize}
         \end{itemize}
     \end{itemize}

\end{itemize}

\pagebreak

The algorithm for {\em the bilinear interpolation} is as below.  This is
similar to the algorithm in the serial code, excepting that local indices
are used to index a segment of the wind component on every processor.

\begin{itemize}
  \item {Let 'seg' denote segment}
  \item {Let 'njj' denote new (local) jj}
  \item { Number of elements in E-W direction for a partition is imp3 }
  \item { Number of elements in N-S direction for a partition is range }
  \item {for j = 1, range}
  \begin {itemize}
    \item {for i = 1, imp3}
    \begin{itemize}
            value at interpolated point =  
             (1.0 - delx(i,j)) * 
                ((1.0 - dely(i,j)) * seg(ii(i,j), njj(i,j)+2)) + 
                  dely(i,j) * seg(ii(i,j), njj(i,j)+3) +
                  delx(i,j) * ((1.0 - dely(i,j)) * seg(ii(i,j)+1, njj(i,j)+2) +
                  dely(i,j) * seg(ii(i,j)+1, njj(i,j)+3))
   \end{itemize}
  \end{itemize}
\end{itemize}

\section{Implementation}

The project was implemented on CM-5 (CMOST 7.3) using the native message 
passing library \cite{cmmd-ref}. This was developed as a proof of concept
that the message passing paradigm was a suitable way of implementing a 
parallel version of the Semi-Lagrangian Transport code.  It is hence
highly modular and has about five levels of subroutine calls.  Though 
this allows for the easy testing of the code, it conflicts with the
requirements of high performance computing in that the overall running 
time drastically increases.  This prompted us in measuring computations
without the function calling overhead.  Therefore the running time of
a production code would be between the timings with the function call
overhead and without the overhead. Ideally we would like it to be 
closer to the timings without the overhead.

The third approach was implemented and 
instrumented for measuring performance as discussed in the next
section.

\pagebreak

\section {Performance}
	The timings were obtained on CM-5 on CMOST 7.3 on a 32 node
partition.  The same code was run on a single node of the CM-5 in order 
to measure the speedup \footnote{ The flow of program execution within a
subroutine is however, not necessarily the same. This results in less 
than ideal speedup for some parts of the code. }. 
The performance is discussed for both bicubic and bilinear
interpolations with respect to the five criteria 
discussed in the Introduction.  They are as follows :
%
%  Change introduction

\begin{itemize}

 \item The overall running time. 
 \item The maximum displacement in terms of nodal distance in the 
       N-S direction.
 \item Ratio of calculation to communication.
 \item Speed-up of interpolation 
 \item Accuracy of calculations 

\end{itemize}


Before discussing the performance issues, we will first define the
following terms :

\begin{itemize}

%Resolution 

%    \(T_{pkt} \) = Amount of time required to communicate a packet of 
%                   specified size.
%
%    \(T_{bnd} \) = Amount of time required to communicate a bound of
%                   the partition

   \item{ Averaged Maximum Number of hops to northern nodes \( NI \). }
   \item{ Averaged Maximum Number of hops to southern nodes \( SI \). }
   \item{ Total Number of hops in both directions \(TI = (NI + SI)\). }
   \item{ The total communication time per iteration of fetching
          or gathering partitions from southern (or northern nodes) \(
          T_{exp} \). This is the sum of the time for
          sending and receiving a bound of the partition plus
          sending and receiving a partition.  This time is derived
          by measuring the communication time for various message sizes. }

   \item{ Time required for bicubic (or bilinear) interpolation
	  on the 32 node partition (without function overhead) \(T_{int} \). }

   \item{ Time required for bicubic (or bilinear) interpolation
          on the single node (without function overhead) \(T_{int_s} \). }

   \item{ Speedup of bicubic (or bilinear) interpolation
           \( S_{int} \) = \(T_{int_s} / T_{int}\). }

   \item{ Time required for reduction of address calculations
	  on the 32 node partition (without function overhead)
          \(T_{ra} \).  This overhead is a result of the implementation
          of the code for the Bicubic Interpolation alone
          and does not show up at the design stage. }

   \item{ Time required for reduction of address calculations
          on the single node (without function overhead) \(T_{ra_s} \). }

   \item{ Speedup of reduction of address calculations \( S_{ra} \) 
           = \(T_{ra_s} / T_{ra}\). }

   \item{ Time required for conversion of the global indices to local 
             indices on the 32 node partition (without function overhead)
          \( T_{li} \). }

   \item{ Time required for conversion of the global indices to local 
          indices on the single node (without function overhead) 
          \( T_{li_s} \). }

   \item{ Speedup of conversion to local indices \( S_{li} \) 
          = \(T_{li_s} / T_{li}\). }

   \item{ Time taken for node-to-node communication (circulation of
          partitions of winds and the bounds of the partition) on the
          32 node partition \(T_{co_n} \). }

   \item{ The product of \(TI\) and \(T_{exp}\) is \(T_{co_{ne}}\).  This
          is approximately equal to \(T_{co_n}\) }.

   \item{ Time required for book-keeping and buffer copying in creation
          of segment (no function calling overhead or node-to-node 
          comunnication is accounted for) on the 32 node partition 
          \( T_{bc} \). }

   \item{ Time required for book-keeping and buffer copying in creation
          of segment (no function calling overhead is accounted for) on 
          the single node \( T_{bc_s} \).  There is no communication required 
          on the single node for creation of segment. }

   \item{ Time required for creation of segment on the 32 node partition 
          \( T_{cs} \) = \( T_{bc} + T_{co_n} \). }

   \item{ Speedup of buffer copying in creation of segment \( S_{bc} \) 
         = \(T_{bc_s} / T_{bc}\).  This indicates how well the buffer 
         copying 
         alone scales when comparing the timings on the 32 node and the 
         single node in CM-5.  The next speedup includes the cost of the
         node-to-node communication which is a part of the ``creation of 
         segment'' function { It is a ``gathering'' process.} . }

   \item{ Speedup of creating a segment \( S_{cs} \) = \(T_{bc_s} / T_{bc}\).
          This includes the node-to-node communication overhead in the 
          computation of the speedup. }

   \item { Time taken for global maximum operation on the 32 node partition
           \( T_{co_g} \). }

   \item{ Time taken for finding local maxima in determining degree of
          circulation (no global reduction involved) on 32 node partition
          \( T_{dcl} \).}

   \item{ Time taken for determining degree of circulation \( T_{dc} \) 
            = \( T_{dcl} + T_{co_g} \). }

   \item{ Time taken for finding local maxima in determining degree of
          circulation on single node partition.  Note the single node
          does not need a global maximum operation \( T_{dcl_s} \).}

   \item{ Speedup of finding local maxima (calculations) \( S_{dcl} \) 
         = \(T_{dcl_s} / T_{dcl}\).  This indicates how well the computation
         alone scales when comparing the timings on the 32 node and the 
         single node in CM-5.  The next speedup includes the cost of the
         global maximum operation. }

   \item{ Speedup of finding degree of circulation \( S_{dc} \) 
         = \(T_{dcl_s} / T_{dc}\).  This speedup includes the cost of 
         the global maximum operation. }

   \item{ Time taken for latitude, longitude index calculations (no 
          function calling overhead) on the 32 node partition \( T_{ll} \). }

   \item{ Time taken for latitude, longitude index calculations (no 
          function calling overhead) on the single node \( T_{ll_s} \). }

   \item{ Speedup of latitude, longitude index calculations \( S_{ll} \) 
         = \(T_{ll_s} / T_{ll}\). }

   \item{ Total calculation (overall) time 
          (with function calling overhead and timing
          overhead) on the 32 node partition \(T_{ca_o}\). 
          This is the sum of \(T_{ca_a}\) and function calling overhead. }

   \item{ Total communication time (includes both global communication 
	  and node to node communication) on the 32 node partition
          \(T_{co}\) = \( T_{co_g} + T_{co_n} \) }

   \item{ Ratio of calculation time (overall) to 
        communication time on the 32-node CM-5
        \(R_{a}\) = \(T_{ca_o} / T_{co}\) }

   \item{ Total calculation time (adjusted, i.e. no function calling 
          overhead) on the 32 node partition 
          \(T_{ca_a}\) = \( T_{ll} + T_{dcl} + T_{bc} + T_{li} + [T_{ra}] + 
          T_{int} \) \footnote{ \(T_{ra}\) is shown in square brackets
         as there is no reduction of address calculations
         for Bilinear Interpolation. Similarly \(T_{ra_s}\) is shown in square
         brackets for the definition of \(T_{ca_{as}}\).}}

   \item{ Ratio of calculation time (adjusted) to 
        communication time on the 32-node CM-5
        \(R_{a}\) = \(T_{ca_a} / T_{co}\). }

   \item{ Total calculation time (adjusted, i.e. no function calling 
          overhead) on the single node \(T_{ca_{as}}\) = 
         \( T_{ll_s} + T_{dcl_s} + T_{bc_s} + T_{li_s} + [T_{ra_s}] 
         + T_{int_s} \). } 

   \item{ Speedup (adjusted) \(S_{a}\) = \(T_{ca_{as}}/ (T_{ca_a} + T_{co})\).}

   \item{ Total (Overall) running time \(T_{ot}\) (includes function calling
 overhead, communication overhead and timing overhead) on a 32 node partition
 on CM5. \( T_{gt} = T_{ca_o} + T_{co} \) .  This is also the sum of 
 \(T_{ca_a}, T_{co}\) and function calling overhead. }

   \item{ Total (Overall) running time \(T_{ot_s}\) (includes function calling
 overhead and timing overhead) on a single node on CM5. This is the sum
 of \(T_{ca_{as}}\) and function calling overhead.}

   \item{ Speedup (overall) \(S_{o}\) = \(T_{gt_s} / T_{gt}\). }

\end{itemize}


The first table shows the non-linear nature of the communication overhead
versus message size.  The overhead decreases as message size increases from 
0 to 8 bytes.  Somewhere between a size of 8 to 16 bytes, the overhead
starts increasing again and again becomes lower at 32 bytes.  The machine 
also exhibits a high latency; for example messages of size 64 bytes takes 
nearly the same amount of time as a message of 256 bytes for a synchronous
send and receive. 

	\paragraph\noindent{
	{\bf Table 1 : Communciation Time }\newline
        \indent {\bf vs Message Size} \newline
	\begin {tabular}{||l|l|l|l||} \hline
	Message Size & Communication \\ 
	(bytes)      & Time (ms) \\ \hline
	0 & 0.226 \\ 
        1 & 0.200 \\
	2 & 0.176 \\
        4 & 0.172 \\
        8 & 0.155 \\ 
        16 & 0.221 \\
        32 & 0.218 \\
        64 & 0.232 \\
        256 & 0.250 \\
        512 & 0.300 \\
        1280 & 0.428 \\
        2048 & 0.799 \\
        5120 & 1.584 \\
        7168 & 1.751 \\  
        10000 & 2.591 \\ 
        22528 & 5.218 \\ \hline
	\end{tabular}} \\\\

Table 2 shows the message size used at a given resolution.  From Tables 1
and 2, we derive the values for \(T_{exp}\).  \(T_{exp}\) is the sum of
the times for communicating the bound of a partition (4 bytes long) and
the partition itself (message size at given resolution).

	\paragraph\noindent{
	{\bf Table 2 : Message Size Used At}\newline
        \indent {\bf Given Resolution} \newline
	\begin {tabular}{||l|l|l|l||} \hline
	 Resolution & Message Size \\ 
	(bytes)      & (bytes) \\ \hline
	(8x64) & 256 \\ 
        (64x64) & 2048 \\
	(128x128) & 5120 \\
        (128x256) & 7168 \\
        (256x512) & 22528 \\  \hline
	\end{tabular}} \\\\


\subsection {Bicubic Interpolation}


%   The third approach to selecting data points is adopted for the bicubic
   The reconstruction of segment approach was adopted for the bicubic
   interpolation.  Table 3 shows the overall results with the function 
   calling overheads.  Table 4 shows the overall results without 
   function calling overhead \footnote{ The function calling overhead 
   does not apply to the timing of the communication routines.} . Tables
   5 to 9 show the functional breakup of the timings in Table 4.\\


{\bf Table 3 : Overall Results } \newline
\begin{tabular}{||l|l|l|l|l|l|l||} \hline
 Resolution&\(T_{ot}\)&\(T_{ot_s}\)&\(S_{o}\)&\(T_{ca_o}\)&\(T_{co}\)
&\(R_{o}\)\\ 
 (I,J)     &  ms      &   ms     &         &   ms       &  ms    & \\ \hline
 (8,64)    &14.62     & 190.00   &13.0     & 11.86  & 2.76 & 4.29\\ 
 (64,64)   &63.60     & 1146.20  &18.0     & 58.41  & 5.18 & 11.26\\
 (128,128) &211.19    & 4500.30  &21.2     & 213.76 & 8.17 & 24.94\\
 (128,256) &396.70    & 9122.00  &23.0     & 388.00 & 8.66 & 44.78\\
 (256,512) &1740.40   & ------   &--      & 1718.5 & 21.91& 78.43\\ \hline
\end{tabular}\\\\


{\bf Table 4 : Results Without Function } \newline
\indent{\bf Calling Overhead} \newline
\begin{tabular}{||l|l|l|l|l||} \hline
Resolution&\(T_{ca_a}\)&\(R_{a}\)&\(T_{ca_{as}}\)& \(S_{a}\)  \\ 
  (I,J)   &  ms        &         & ms     &         \\ \hline
 (8,64)   & 4.85       & 1.75    & 85.1   & 11.18   \\ 
 (64,64)  & 26.42      & 5.09    & 522.8  & 17.5    \\
(128,128) & 92.04      & 11.27   & 1985.00& 19.8    \\
(128,256) & 173.45     & 20.01   & 4034.00& 22.15   \\
(256,512) & 714.75     & 32.62   & ------ & -----   \\ \hline
\end{tabular}\\\\


{\bf Table 5 : Functional Breakup } \newline
\indent{ \bf of Table 2 Timings } \newline
\begin{tabular}{||l|l|l|l||} \hline
Resolution&\(T_{ll}\)&\(T_{ll_s}\)& \(S_{ll}\)  \\ 
  (I,J)   &  ms      & ms         &         \\ \hline
 (8,64)   & 0.58     & 15.16      & 26.14   \\ 
 (64,64)  & 3.20     & 92.34      & 28.86   \\
(128,128) & 12.76    & 364.7      & 28.58   \\
(128,256) & 27.57    & 732.5      & 26.56   \\
(256,512) & 99.24    & -----      & -----   \\ \hline
\end{tabular}\\\\


\newpage

{\bf Table 6 : Functional Breakup of Table 2 } \newline
\indent{\bf Timings } \newline
\begin{tabular}{||l|l|l|l|l|l|l||} \hline
Resolution&\(T_{co_g}\)&\(T_{dcl}\)&\(T_{dc}\)&\(T_{dcl_s}\)&
\(S_{dcl}\)&\(S_{dc}\) \\ 
  (I,J)   &  ms     & ms    &   ms   &   ms   &       &     \\ \hline
 (8,64)   & 0.0774  & 0.59  & 0.667  & 13.82  & 23.42 & 20.70\\ 
 (64,64)  & 0.0829  & 3.30  & 3.38   & 77.26  & 23.41 & 22.88\\
(128,128) & 0.0930  & 12.99 & 13.08  & 305.68 & 23.53 & 23.37\\
(128,256) & 0.0910  & 26.29 & 26.38  & 613.21 & 23.32 & 23.24\\
(256,512) & 0.1070  & 108.8 & 108.9  & ------ & ----- & -----\\ \hline
\end{tabular}\\\\


{\bf Table 7 : Functional Breakup of Table 2 } \newline
\indent{\bf Timings } \newline
\begin{tabular}{||l|l|l|l|l|l|l||} \hline
Resolution&\(T_{co_n}\)&\(T_{bc}\)&\(T_{cs}\)&\(T_{bc_s}\)&
\(S_{bc}\)&\(S_{cs}\) \\ 
  (I,J)   &  ms   & ms    &   ms   &   ms   &      &     \\ \hline
 (8,64)   & 2.69  & 1.51  & 4.2    & 5.52   & 3.66 & 1.31\\ 
 (64,64)  & 5.11  & 6.81  & 11.92  & 33.65  & 4.94 & 2.82\\
(128,128) & 8.08  & 17.31 & 25.39  & 128.77 & 7.44 & 5.07\\
(128,256) & 8.57  & 24.41 & 32.98  & 293.8  & 12.04& 8.91\\
(256,512) & 21.8  & 87.35 & 109.15 & ------ & -----& ----\\ \hline
\end{tabular}\\\\


{\bf Table 8 : Functional Breakup of Table 2 } \newline
\indent{\bf Timings } \newline
\begin{tabular}{||l|l|l|l|l|l|l||} \hline
Resolution&\(T_{li}\)&\(T_{li_s}\)&\(S_{li}\)&\(T_{ra}\)&
\(T_{ra_s}\)&\(S_{ra}\) \\ 
  (I,J)   &  ms   & ms     &        &  ms   &  ms   &      \\ \hline
 (8,64)   & 0.274 & 2.16   & 7.88   & 0.85  & 15.03 & 17.68\\ 
 (64,64)  & 1.55  & 13.77  & 8.88   & 5.42  & 89.55 & 16.52\\
(128,128) & 4.04  & 53.17  & 13.16  & 20.93 & 360.56& 17.22\\
(128,256) & 5.55  & 109.65 & 19.75  & 42.80 & 734.56& 17.16\\
(256,512) & 17.2  & ------ & ------ & 90.92 & ------& -----\\ \hline
\end{tabular}\\\\


{\bf Table 9 : Functional Breakup of Table 2 } \newline
\indent{\bf Timings } \newline
\begin{tabular}{||l|l|l|l||} \hline
Resolution&\(T_{int}\)&\(T_{int_s}\)&\(S_{int}\) \\ 
  (I,J)   &  ms   & ms     &        \\ \hline
 (8,64)   & 1.04  & 34.02  & 32.7   \\ 
 (64,64)  & 6.13  & 216.25 & 35.28  \\
(128,128) & 24.02 & 777.20 & 32.15  \\
(128,256) & 46.81 & 1550.9 & 33.12  \\
(256,512) & 311.29& ------ & ------ \\ \hline
\end{tabular}\\\\

\newpage
{\bf Table 10 : Experimental Confirmation of } \newline
\indent{\bf Node-to-Node Communication Time \(T_{co_n}\) } \newline
\begin{tabular}{||l|l|l|l|l|l|l||} \hline
Resolution&  NH  &  SH  &  TH  &\(T_{exp}\) & \(T_{co_ne}\) \\ 
  (I,J)   &      &      &      &     ms     &  ms           \\ \hline
 (8,64)   & 4.0  &  1.5 &  5.5 &   0.422    &    2.32       \\ 
 (64,64)  & 4.2  &  1.5 &  5.7 &   0.971    &    5.53       \\
(128,128) & 3.1  &  2.0 &  5.1 &   1.756    &    8.96       \\
(128,256) & 2.0  &  2.0 &  4.0 &   1.923    &    7.69       \\
(256,512) & 2.0  &  2.0 &  4.0 &   5.390    &   21.56       \\ \hline
\end{tabular}\\\\


   Since the code has 5 levels of subroutine calling, there is approximately
   a factor of 2 difference between the overall timings with and without
   the function calling overhead.  A good production code would have a 
   {\em overall timing} somewhere between these two figures
   \footnote{ It is not
   desirable to eliminate totally the function calling overhead except where
   it is in excess.  For example, one should inline the calls to the actual
   interpolation routine as it is called within an inner loop of a doubly 
   nested loop. }. 

   The {\em inter-nodal distance} seems to be fairly static for both components
   (north, south directions inclusive). On an average this is 
   4 "cshifts" in the Northern direction or 1.5 "cshifts" in the 
   Southern for a resolution of 8 x 64.  This appears to stabilize at 
   2 ``cshifts'' in either direction at higher resolutions as seen from 
   Table 10.   This is close to the assumption that the
   winds generate mostly nearest neighbor communication.

   The {\em ratio of calculation to communication} ( \(R_{o}\) or \(R_{a}\) )
   (Tables 3,4) falls with increasing resolution because of two factors, 
    namely :
   \begin{itemize}
   \item { High latency in the communciation overhead.
   Due to high latency, it takes almost the same time to send and receive
   a message of 256 bytes as a message of size 64 bytes.}
   \item { The decreasing total number of cshifts \(TI\).  The total number
   of ``cshifts'' falls from 5.5 to 4. }
   \end{itemize}
   It is obvious that this algorithm is dominated by overheads other than 
   communication at high resolutions, which will be established shortly 
   \footnote{ Communication 
   overhead is only 1.28 \% at a resolution of 128x256. }.  

   It is also noticed since the function calling overhead is high, the ratio
   \(R_{a}\) is less than \(R_{o}\).  Consequently, the {\em speedup}
   \(S_{a}\) is
   lower than \(S_{o}\).  However, \(S_{o}\) itself is much lower than 
   ideal.  At low resolutions, this is primarily due to the high communication 
   overhead \footnote{ At a resolution of 8 x 64, the ratio \(R_{o}\) is as 
   low as 4.29. }, but does not hold true at higher resolutions.  Hence, a 
   complete picture of speedup will not emerge without a detailed discussion
   at the functional level of the code.  We will first tabulate the functions
   in order of increasing speedup for the resolution of 8x64 (Table 11)
   from the data in 
   Tables 3 to 7 and then analyze the reasons for the performance. We
   will do this similarly for a higher resolution of 128x256 (Table 12).

\newpage
	\paragraph\noindent{
	{\bf Table 11 : Functional Breakup of Timings at a } \newline
        \indent {\bf Resolution of 8 x 64: } \newline
	\begin {tabular}{||l|l|l|l||} \hline
	 Function & Speedup & \%age of \(T_{ca_a} + T_{co}\)
         & \%age of \(T_{ot}\) \\ \hline
   Creation of Segment (cs) &  1.31 & 53.98 \% & 28.1 \%  \\ 
   Buffer Copying in above (bc) &  3.66 & 19.79 \% & 10.3 \%  \\
   Conversion to Local Indices (li) & 7.88 & 3.65 \% & 1.9 \% \\
   Reduction of Address Calculation (ra) & 17.68 & 11.14 \% & 5.8 \% \\
   Degree of Circulation (dc) & 20.70 & 8.83 \% & 4.6 \% \\
   Computation alone in above (dcl) & 23.42 & 7.68 \% & 4.0 \% \\
   Latitude, Longitude Index Calculations (ll) & 26.14 & 7.68 \% & 4.0 \%  \\
   Interpolation (int) & 32.7 & 13.45 \% & 7.0 \% \\ \hline 
	\end{tabular}} \\\\
   
	The ``creation of segment'' function shows a very low speedup 
  primarily because it involves node-to-node communication with high 
  latency.  But, this does not account for all the loss in speedup 
  (speedup of buffer copying is also at a low of 3.66).  This is
  because there are more places in the 32-node run where there is
  {\em buffer copying} \footnote{ Though both the single node and the 32 node
  run the same code, because the single node has no communication to 
  perform, it does not need to perform the buffer copying associated
  with the sending and receiving of partitions. In other words the single
  node program traces a different path through the code.  We will hereafter
  refer to this as a {\em run}. }.  The term buffer copying includes 
  memory to memory copying and the book-keeping associated with it. 
  There is some load imbalance in the memory to memory copying, especially
  at the low resolutions, due to nodes at the poles having to deal with 
  a larger partition (8x5 instead of 8x2) as shown in Figure 2. These 
  nodes also have to handle
  with the pole-wrap problem \footnote{ Fetching of partition by wrapping 
  around the
  pole, needs more book-keeping.} which are absent in other nodes, leading to
  an additional load imbalance.  The load imbalance
  is apparent from the fact that 
  at the next higher resolution (64x64), though
  the grid points have increased by 8,
  the amount of time spent in buffer copying increases by 4. Otherwise
  the buffer copying scales with number of grid points.  In the single node
  we find that the time spent in buffer copying scales well with the number
  of grid points.

  Taking into account the implementation
  and the sizes of buffers involved in both the runs, the single node run
  performs only about 6 times the amount \footnote{ Not to be confused
  with number of places} of buffer copying as the 32 node run.  This is
  slightly
  more than the speedup of 3.66. However, due to the load imbalance the drop 
  in speedup can be expected .  As the buffer copying alone takes 19.79 \% of
  running time (without function calling overhead), this has a significant 
  effect on the overall speedup 
  \footnote{ Inclusive of the communication cost, it takes the lion's share
  of the running time at 53.98 \%.} .

  The speedup of ``conversion to local indices'' is also low.  This 
  could be due to an error in measurement of a very small interval 
  of time (0.27 ms) in the 32-node run.  There is also the possibility
  that due to a load imbalance that affects the nodes corresponding to
  the polar regions  which have to translate indices that cross the poles, 
  that the speedup is low.  This code involves testing of a condition in the
  innermost loop of a doubly nested loop.  This does not occur in the
  single node run where only a straight forward translation to local indices
  takes place.  But, the impact of this function is slight as it
  constitutes only 3.65 \% of the running time (which is the least of all the
  functions).

  The ``reduction of address calculations'' constitutes about 11.14 \%
  of the running time.  The number of times it is called for the 
  32-node run is equal to the number of grid points in a partition.
  There are 16 grid points at a resolution of 8x64.  Therefore the
  approximate time per call is only 53 \(\mu\)s which spans
  8 integer operations.  The measurement of
  such a small time is prone to errors, which may have led to a 
  somewhat moderate speedup of 17 to 18. \footnote{ Even the timings
  for the single node run suffer from the same defect.}

  The ``degree of circulation'' has a moderate speedup of 20.7.  When the
  global maximum operation is left out, the speedup is about 23.4, which
  is fairly high. As the global maximum operation is a constant overhead, 
  its effect at higher resolutions becomes negligible.
  As the code involves testing of conditions and as the quantities involved
  are dynamically varying with both time and space, it is difficult to
  pinpoint how the speedup is affected.  Additionally this function has
  no significance for the single node run.

  The ``calculation of latitude, longitude indices'' shows a good speedup
  of 26.1.  This involves a mix of integer and floating point arithmetic.
 
  The actual ``bicubic interpolation'' shows more than ideal speedup.
  The timing of this routine has the same drawback as the ``reduction of
  address calculations''.  As the time intervals involved are very small,
  about 65 \(\mu\)s, the measurements are prone to errors.  However, it
  may be concluded that this function has a good speedup even after 
  discounting for any errors.
  
  The overall speedup \(S_{a}\) is therefore poor because 
  functions which have poor speedups (\(< 16\)) constitute 57.6 \% of
   running time
  (with no function calling overhead), those with medium speedup 
  (\(16 \le speedup \le 23\)) constitute 20.0 \% and those with good 
  speedups 22.2 \% (\(> 23\)) \footnote{ The classification of poor, medium and
  good speedups may seem arbitrary.  But this distinction serves the purpose
  considering that timings were obtained on a 32 node machine. }.

	\paragraph\noindent{
	{\bf Table 12 : Functional Breakup of Timings at a } \newline
        \indent {\bf Resolution of 128x256: }\newline
	\begin {tabular}{||l|l|l|l||} \hline
	 Function & Speedup & \%age of \(T_{ca_a} + T_{co}\) &
         \%age of \(T_{ot}\) \\ \hline
   Creation of Segment (cs) &  8.91 & 18.1 \% & 8.31 \%  \\ 
   Buffer Copying in above (bc) &  12.04 & 13.40 \% & 6.15 \% \\
   Reduction of Address Calculation (ra) & 17.16 & 23.53 \% & 10.8 \% \\
   Conversion to Local Indices (li) & 19.75 & 3.05 \% & 1.4 \%  \\
   Degree of Circulation (dc) & 23.24 & 14.49 \% & 6.65 \%  \\
   Computation alone in above (dcl) & 23.32 & 14.42 \% & 6.62 \%  \\
   Latitude, Longitude Index Calculations (ll) & 26.56 & 15.14 \% & 6.95 \% \\
   Interpolation (int) & 33.1 & 25.70 \% & 11.8 \% \\ \hline 
	\end{tabular}} \\\\
   
	At this resolution cache miss penalties play an important role as 
   seen from Tables 5 to 9 where there is a jump in speedups from the 128x128
   to 128x256 resolution \footnote{ It was not possible to
   run at a resolution of 256x256 on the single node due to lack of memory. }. 

	Using the same reasoning as in the case of 8x64 resolution for the
   ``creation of segment'' function,
   the single-node run has about 6 times the amount of {\em buffer 
   copying} as the 32 node-run.  However, the speedup of buffer copying
   is 12.04 which is two times more than expected.  This is most
   likely due to caching effects.  It can also be seen that there is an
   increase of 1-1/2 times in \(T_{bc}\) compared with a two-fold increase 
   in the number of grid points from a resolution of 128x128.  This indicates
   that the load imbalance due to unven partition size and pole wrapping
   is becoming less severe compared to a resolution of 8x64.

	The node-to-node communication is still expensive as the speedup
   inclusive of this cost falls to 8.91 from 12.04. However, the adverse
   effect of the ``creation of segment'' function is less on the overall 
   speedup as it only takes 18.1 \% of the running time compared to 53.98 \%
   at the 8x64 resolution. 

  The speedup of ``conversion to local indices'' increases to 19.75.  
  This is certainly due to cache misses, because at a resolution of 
  128x128 the speedup is only 13.2.  At the 128x128 resolution the
  accuracy of the timing measurements should be good as
  there are 512 iterations taking place.  The speedup is not better  
  than 19.75 probably because of the load imbalance problem
  at the nodes corresponding to the poles.

  The measurement of the total time for the ``reduction of 
  calculations'' suffers from the
  drawback discussed earlier.  The speedup for this function shows little
  variance (Table 8).  But since it now constitutes 10.8 \%, it does
  little to contribute to a better speedup.

  The ``degree of circulation'' has a fairly high speedup of 23.24. The 
  global maximum operation now constitutes a negligible amount of the time
  spent in the function, as seen in the difference in speedups \(S_{dcl}\) 
  and \(S_{dc}\).  

  The ``calculation of latitude, longitude indices'' shows a good speedup
  of 26.56.  
 
  The actual ``bicubic interpolation'' as before shows more than
  ideal speedup due to the drawback discussed earlier.

  
  Thus the overall speedup \(S_{a}\) is better because 
  functions which take major part of running time (55.35 \%) have good
  speedups (\(> 23\)).  About 12 \% of the code has medium 
  speedup ( \(16 \le speedup \le 23\) ) and the rest has a poor speedup 
  (\(< 16\)).


  An {\em accuracy} of upto (but not including) the fifth decimal 
  place was obtained between the parallel and serial implementations.

  Table 10 confirms the {\em accuracy of the node-to-node communication
  overhead }.  The results from the secondary experiment, \(T_{co_{ne}}\)
  and \(T_{co_n}\) match closely.

\newpage
\subsection{Bilinear Interpolation}


%   The third approach to selecting data points is adopted for the bilinear
   The reconstruction of segment approach was also adopted for the bilinear
   interpolation.  Table 13 shows the overall results with the function 
   calling overheads.  Table 14 shows the overall results without 
   function calling overhead \footnote{ The function calling overhead 
   does not apply to the timing of the communication routines.} . Tables
   15 to 19 show the functional breakup of the timings in Table 14. We
   obtain very similar results for certain functions in the code, because
   the code is essentially the same excepting for the actual interpolation
   routine.  There is also no "reduction of address calculation" function 
   in this code. \\

{\bf Table 13 : Overall Results } \newline
\begin{tabular}{||l|l|l|l|l|l|l||} \hline
 Resolution&\(T_{ot}\)&\(T_{ot_s}\)&\(S_{o}\)&\(T_{ca_o}\)&\(T_{co}\)
&\(R_{o}\)\\ 
 (I,J)     &  ms      &   ms     &         &   ms       &  ms    & \\ \hline
 (8,64)    &10.50     & 96.88    &9.23     & 8.22   & 2.27 & 3.62\\ 
 (64,64)   &41.12     & 589.6    &14.34    & 36.63  & 4.49 & 8.176\\
 (128,128) &128.46    & 2280.0   &17.74    & 112.03 & 6.43 & 18.97\\
 (128,256) &240.47    & 4717.0   &19.60    & 232.05 & 8.42 & 27.56\\
 (256,512) &955.20    & ------   &--     & 930.68 & 24.53& 37.95\\ \hline
\end{tabular}\\\\


{\bf Table 14 : Results Without Function } \newline
\indent{\bf Calling Overhead} \newline
\begin{tabular}{||l|l|l|l|l||} \hline
Resolution&\(T_{ca_a}\)&\(R_{a}\)&\(T_{ca_{as}}\)& \(S_{a}\)  \\ 
  (I,J)   &  ms        &         & ms     &         \\ \hline
 (8,64)   & 3.37  & 1.48    & 47.31  &  8.39   \\ 
 (64,64)  & 16.17 & 3.60    & 290.3  & 14.05   \\
(128,128) & 52.63 & 8.18    & 1142.0 & 19.34   \\
(128,256) & 98.37 & 11.68   & 2349.1 & 22.03   \\
(256,512) & 381.45& 15.55   & ------ & -----   \\ \hline
\end{tabular}\\\\


{\bf Table 15 : Functional Breakup } \newline
\indent{ \bf of Table 2 Timings } \newline
\begin{tabular}{||l|l|l|l||} \hline
Resolution&\(T_{ll}\)&\(T_{ll_s}\)& \(S_{ll}\)  \\ 
  (I,J)   &  ms      & ms         &         \\ \hline
 (8,64)   & 0.61     & 14.72      & 24.09   \\ 
 (64,64)  & 2.95     & 91.87      & 31.14   \\
(128,128) & 11.19    & 364.8      & 32.60   \\
(128,256) & 22.53    & 733.0      & 32.60   \\
(256,512) & 91.23    & -----      & -----   \\ \hline
\end{tabular}\\\\

\newpage
{\bf Table 16 : Functional Breakup of Table 2 } \newline
\indent{\bf Timings } \newline
\begin{tabular}{||l|l|l|l|l|l|l||} \hline
Resolution&\(T_{co_g}\)&\(T_{dcl}\)&\(T_{dc}\)&\(T_{dcl_s}\)&
\(S_{dcl}\)&\(S_{dc}\) \\ 
  (I,J)   &  ms     & ms    &   ms   &   ms   &       &     \\ \hline
 (8,64)   & 0.0756  & 0.60  & 0.67   & 11.55  & 19.40 & 17.21\\ 
 (64,64)  & 0.0797  & 3.27  & 3.35   & 69.09  & 21.12 & 20.64\\
(128,128) & 0.0907  & 12.88 & 12.97  & 276.50 & 21.47 & 21.31\\
(128,256) & 0.0890  & 25.72 & 25.81  & 556.30 & 21.63 & 21.55\\
(256,512) & 0.0946  & 102.7 & 102.79 & ------ & ----- & -----\\ \hline
\end{tabular}\\\\


{\bf Table 17 : Functional Breakup of Table 2 } \newline
\indent{\bf Timings } \newline
\begin{tabular}{||l|l|l|l|l|l|l||} \hline
Resolution&\(T_{co_n}\)&\(T_{bc}\)&\(T_{cs}\)&\(T_{bc_s}\)&
\(S_{bc}\)&\(S_{cs}\) \\ 
  (I,J)   &  ms   & ms    &   ms   &   ms   &      &     \\ \hline
 (8,64)   & 2.20  & 1.47  & 3.67   & 5.50   & 3.74 & 1.50\\ 
 (64,64)  & 4.41  & 5.98  & 10.39  & 34.59  & 5.78 & 3.33\\
(128,128) & 6.34  & 14.50 & 20.84  & 128.69 & 8.88 & 6.18\\
(128,256) & 8.33  & 24.66 & 32.99  & 296.50 & 12.02& 8.99\\
(256,512) & 24.43 & 87.54 & 111.97 & ------ & -----& ----\\ \hline
\end{tabular}\\\\


{\bf Table 18 : Functional Breakup of Table 2 } \newline
\indent{\bf Timings } \newline
\begin{tabular}{||l|l|l|l||} \hline
Resolution&\(T_{li}\)&\(T_{li_s}\)&\(S_{li}\) \\
  (I,J)   &  ms   & ms     &                  \\ \hline
 (8,64)   & 0.28  & 2.19   & 7.96             \\
 (64,64)  & 1.56  & 13.63  & 8.74             \\
(128,128) & 3.86  & 53.13  & 13.76            \\
(128,256) & 5.44  & 108.9  & 20.02            \\
(256,512) & 17.2  & ------ & ------           \\ \hline
\end{tabular}\\\\


{\bf Table 19 : Functional Breakup of Table 2 } \newline
\indent{\bf Timings } \newline
\begin{tabular}{||l|l|l|l||} \hline
Resolution&\(T_{int}\)&\(T_{int_s}\)&\(S_{int}\) \\ 
  (I,J)   &  ms   & ms     &        \\ \hline
 (8,64)   & 0.410 & 13.36  & 32.59  \\ 
 (64,64)  & 2.40  & 81.11  & 33.8   \\
(128,128) & 10.21 & 318.78 & 31.22  \\
(128,256) & 19.25 & 654.56 & 32.81  \\
(256,512) & 82.81 & ------ & ------ \\ \hline
\end{tabular}\\\\

\newpage
{\bf Table 20 : Experimental Confirmation of } \newline
\indent{\bf Node-to-Node Communication Time \(T_{co_n}\) } \newline
\begin{tabular}{||l|l|l|l|l|l|l||} \hline
Resolution&  NH  &  SH  &  TH  &\(T_{exp}\) & \(T_{co_ne}\) \\ 
  (I,J)   &      &      &      &     ms     &  ms           \\ \hline
 (8,64)   & 3.33 & 1.67 &  5.0 &  0.422     &   2.11        \\ 
 (64,64)  & 3.33 & 1.67 &  5.0 &  0.971     &   4.86        \\
(128,128) & 2.00 & 1.93 & 3.93 &  1.756     &   6.90        \\
(128,256) & 2.0  &  2.0 &  4.0 &  1.923     &   7.692       \\
(256,512) & 2.0  &  2.0 &  4.0 &  5.390     &   21.56       \\ \hline
\end{tabular}\\\\

   Since the code has 5 levels of subroutine calling, there is approximately
   a factor of 2 difference between the overall timings with and without
   the function calling overhead.  A good production code would have a 
   {\em overall timing} somewhere between these two figures.  We would like
   the timings to be closer to the value of \(T_{ca_a} + T_{co}\) than
   \(T_{ot}\).

   The {\em inter-nodal distance} seems to be fairly static for both components
   (north, south directions inclusive). On an average this is 
   3.33 "cshifts" in the Northern direction or 1.67 "cshifts" in the 
   Southern for a resolution of 8 x 64.  This appears to stabilize at 
   2 ``cshifts'' in either direction at higher resolutions as seen from 
   Table 10.   This is close to the assumption that the
   winds generate mostly nearest neighbor communication.

   The {\em ratio of calculation to communication} ( \(R_{o}\) or \(R_{a}\) )
   (Tables 13,14) falls with increasing resolution because of two factors, 
   namely :
   \begin{itemize}
   \item { High latency in the communciation overhead.
   As observed earlier, it takes almost the same time to send and receive
   a message of 256 bytes as a message of size 64 bytes.}
   \item { The decreasing total number of cshifts \(TI\).  The total number
   of ``cshifts'' falls from 5.0 to 4. }
   \end{itemize}
   It is obvious that this algorithm is dominated by overheads other than 
   communication at high resolutions, which will be established shortly 
   \footnote{ Communication 
   overhead is only 2.64 \% at a resolution of 128x256. }.  

   It is also noticed since the function calling overhead is high, the ratio
   \(R_{a}\) is less than \(R_{o}\).  Consequently, the {\em speedup}
   \(S_{a}\) is
   lower than \(S_{o}\).  However, \(S_{o}\) itself is much lower than 
   ideal.  At low resolutions, this is primarily due to the high communication 
   overhead \footnote{ At a resolution of 8 x 64, the ratio \(R_{o}\) is as 
   low as 3.62, whereas at 128x256 \(R_{o}\) reaches a high of 37.95. },
    but does not hold true at higher resolutions.  Hence, a 
   complete picture of speedup will not emerge without a detailed discussion
   at the functional level of the code. \footnote{ Another notable aspect 
   is that the overall speedups for the Bilinear Interpolation is lower than 
   that for Bicubic Interpolation.  This is because the actual interpolation
   routine for the Bilinear Interpolation involves only 6 floating-point 
   multiplications compared to 36 for the Bicubic Interpolation. }. 
   We will first tabulate the functions
   in order of increasing speedup for the resolution of 8x64 from the data in 
   Tables 13 to 19 (Table 21) and then analyze the reasons for the performance.
   We will do this similarly for a higher resolution of 128x256 (Table 22).

	\paragraph\noindent{
	{\bf Table 21 : Functional Breakup of Timings at a } \newline
        \indent {\bf Resolution of 8 x 64: } \newline
	\begin{tabular}{||l|l|l|l||} \hline
	 Function & Speedup & \%age of \(T_{ca_a} + T_{co}\)
         & \%age of \(T_{ot}\) \\ \hline
   Creation of Segment (cs) &  1.50 & 65.07 \% & 34.95 \%  \\ 
   Buffer Copying in above (bc) &  3.74 & 26.06 \% & 14.00 \%  \\
   Conversion to Local Indices (li) & 7.96 & 4.96 \% & 2.67 \% \\
   Degree of Circulation (dc) & 17.21 & 11.88 \% & 6.38 \% \\
   Computation alone in above (dcl) & 19.40 & 10.64 \% & 5.7 \% \\
   Latitude, Longitude Index Calculations (ll) & 24.09 & 10.82 \% & 5.8 \%  \\
   Interpolation (int) & 32.59 & 7.27 \% & 3.90 \% \\ \hline 
	\end{tabular}} \\\\
   
	The ``creation of segment'' function shows a very low speedup 
  primarily because it involves node-to-node communication with high 
  latency.  As discussed in the case of the Bicubic Interpolation, this does 
  not account for all the loss in speedup. There is some load imbalance 
  in the memory to memory copying, especially
  at the low resolutions, due to nodes at the poles having to deal with 
  a larger partition (8x5 instead of 8x2) as shown in Figure 2. These 
  nodes also have to handle with the pole-wrap problem which are absent
  in other nodes, leading to an additional load imbalance.  The load imbalance
  is apparent from the fact that  at the next higher resolution (64x64), though
  the grid points have increased by 8,
  the amount of time spent in buffer copying increases by 4. Otherwise
  the buffer copying scales with number of grid points.  In the single node
  we find that the time spent in buffer copying scales well with the number
  of grid points.

  The speedup of buffer copying is 
  also low, as there are more places in the 32-node run where there is
  buffer copying.   Taking into account the implementation
  and the sizes of buffers involved in both the runs, the single node run
  performs only about 6 times the amount \footnote{ Not to be confused
  with number of places} of buffer copying as the 32 node run.  This is 
  slightly more than the speedup of 3.66, but this is expected due to 
  load imbalance.  As the buffer copying alone takes 
  26.06 \% of running time (without function calling overhead), this has a
  significant effect on the overall speedup 
  \footnote{ Inclusive of the communication cost, it takes the lion's share
  of the running time at 65.07 \%.} .

  The speedup of ``conversion to local indices'' is also low, like the 
  Bicubic Interpolation.  This could be due to an 
  error in measurement of a very small interval 
  of time (0.28 ms) in the 32-node run.  There is also the possibility
  that due to a load imbalance that affects the nodes corresponding to
  the polar regions  which have to translate indices that cross the poles, 
  that the speedup is low.  But, the impact of this function is slight as it
  constitutes only 4.96 \% of the running time (which is the least of all the
  functions).

  The ``degree of circulation'' has a moderate speedup of 17.21.  When the
  global maximum operation is left out, the speedup is about 19.40, which
  is medium. As the global maximum operation is a constant overhead, 
  its effect at higher resolutions becomes negligible.
  As the code involves testing of conditions and as the quantities involved
  are dynamically varying with both time and space, it is difficult to
  pinpoint how the speedup is affected.  Additionally this function has
  no significance for the single node run.

  The ``calculation of latitude, longitude indices'' shows a good speedup
  of 24.09 which involves a mix of integer and floating point arithmetic.
 
  The actual ``bilinear interpolation'' shows more than ideal speedup.
  The timing of this routine has the same drawback as the ``reduction of
  address calculations''.  As the time intervals involved are very small,
  the measurements are prone to errors.  However, it
  may be concluded that this function has a good speedup even after 
  discounting for any errors.
  
  The overall speedup \(S_{a}\) is therefore poor because 
  functions which have poor speedups (\(< 16\)) constitute 70.02 \% of 
  running time
  (with no function calling overhead), those with medium speedup 
  (\(16 \le speedup \le 23\)) constitute 11.88 \% and those with good 
  speedups 18.1 \% (\(> 23\)) \footnote{ The classification of poor, medium and
  good speedups may seem arbitrary.  But this distinction serves the purpose
  considering that timings were obtained on a 32 node machine. }.

	\paragraph\noindent{
	{\bf Table 22 : Functional Breakup of Timings at a } \newline
        \indent {\bf Resolution of 128x256: }\newline
	\begin{tabular}{||l|l|l|l||} \hline
	 Function & Speedup & \%age of \(T_{ca_a} + T_{co}\) &
         \%age of \(T_{ot}\) \\ \hline
   Creation of Segment (cs) &  8.99 & 30.89 \% & 13.72 \%  \\ 
   Buffer Copying in above (bc) &  12.02 & 23.09 \% & 10.25 \% \\
   Conversion to Local Indices (li) & 20.02 & 5.09 \% & 2.20 \%  \\
   Degree of Circulation (dc) & 21.55 & 24.17 \% & 10.73 \%  \\
   Computation alone in above (dcl) & 21.63 & 24.08 \% & 10.70 \%  \\
   Latitude, Longitude Index Calculations (ll) & 32.60 & 21.09 \% & 9.37 \% \\
   Interpolation (int) & 32.81 & 18.02 \% & 8.01 \% \\ \hline 
	\end{tabular}} \\\\

	As discussed in the case of the Bicubic Interpolation, at this
   resolution, the jumps in speedups compared to the 128x128 resolution
   is due to cache misses. \footnote{ It was not possible to
   run at a resolution of 
   256x256 on the single node due to lack of memory. }
   
	Using the same reasoning as in the case of 8x64 resolution,
   the single-node run has about 6 times the amount of buffer 
   copying as the 32 node-run for the "creation of segment" function.
   However, the speedup of buffer copying is 12.04 which is two times more
   than expected.  This is most likely due to caching effects.
   It can also be seen that there is an
   increase of 1.7 times in \(T_{bc}\) compared with a two-fold increase 
   in the number of grid points from a resolution of 128x128.  This indicates
   that the load imbalance due to unven partition size and pole wrapping
   is becoming less severe compared to a resolution of 8x64.

	The node-to-node communication is still expensive as the speedup
   inclusive of this cost falls to 8.99 from 12.02. As the 
   "creation of segment" function still takes 30.89 \% 
   (13.72 \% with function calling overhead) of the
   running time, it affects the overall speedup adversely.

  The speedup of ``conversion to local indices'' increases to 20.02.  
  This is certainly due to cache misses, because at a resolution of 
  128x128 the speedup is only 13.76.  At the 128x128 resolution the
  accuracy of the timing measurements should be good as
  there are 512 iterations taking place.  The speedup is not better  
  than 20.0 because of the load imbalance problem at the nodes 
  corresponding to the poles.
% remove probably !!

  The ``degree of circulation'' has a fairly high speedup of 21.55 The 
  global maximum operation now constitutes a negligible amount of the time
  spent in the function, as seen in the difference in speedups \(S_{dcl}\) 
  and \(S_{dc}\).  

  The ``calculation of latitude, longitude indices'' shows more than ideal
  speedup of 32.60.  This is due to errors in measurement of timings.
 
  The actual ``bicubic interpolation'' as before shows more than
  ideal speedup due to the drawback discussed earlier.

  
  Thus the overall speedup \(S_{o}\) or \(S_{a}\) is better than at a 
  resolution of 8x64. The speedup is only medium because functions 
  which take major part of running time (60.15 \%) have poor or medium
  speedups ( \( \le 23\) ).  About 30.89 \% of the code has poor 
  speedup (\(< 16\))
  , another 29.26 \%  has medium speedup ( \(16 \le speedup \le 23\) )
  and the balance 39.85 \% the rest has a good speedup (\(> 23\)).

  An {\em accuracy} of upto (but not including) the fifth decimal 
  place was obtained between the parallel and serial implementations.

  Table 20 confirms the {\em accuracy of the node-to-node communication
  overhead }.  The results from the secondary experiment, \(T_{co_{ne}}\)
  and \(T_{co_n}\) match closely.


\section {Conclusions}

%	 A reasonable speedup maybe achieved by the message passing version.
%It is apparent that this problem is hard to parallelize because the overhead 
%of the book-keeping and memory management required for the parallel version
%is high.  The communication overhead
%and load imbalance become less important at higher resolutions. However, there
%is scope for improving the overall running time by reducing the function call
% overhead.

	 A reasonable speedup maybe achieved by the message passing version.
This problem is hard to parallelize because of the following reasons.  At low 
resolutions, the speedup is poor
because of very high communication overhead \footnote{
For the bicubic interpolation, \(T_{co}\) takes 57 \% as much time as
 \(T_{ca_a}\).  For the
bilinear interpolation, \(T_{co}\) takes 68 \% as much time as \(T_{ca_a}\), 
which is even worse. }, some load imbalance, and a fair
amount of 
book-keeping and memory management involved.  At higher resolutions the
amount of book-keeping and memory management dominate the communication
\footnote{ \(T_{co}\) takes only about 30 \% as much time as \(T_{bc}\) for 
either the bicubic 
or the bilinear interpolations. \(T_{bc}\) itself takes about 14 \% and
24.9 \% of \(T_{ca_a}\) for the bicubic and bilinear interpolations 
respectively. } and load imbalance overheads.  But, unfortunately, these 
overheads are required for the working of the parallel algorithm and are
difficult to reduce.   Hence, there is only a medium speedup.  If the
algorithm were to be ported to a machine with a low latency communication 
network,  the book keeping and memory managment overheads will still result
in only a medium speedup.  However,
 there is scope for improving the overall running time by reducing the 
function call overhead.

\pagebreak
         
\section {Further Work}

	In the short term further testing needs to be done to verify
the accuracy and stability of the parallel code.  The scaling of the
performance of the parallel code with the number of processors needs
to be observed.

	In the long term more, more scalable decompositions such as
the igloo decomposition have to be studied.  A different approach will
be needed.  The "cshift" operation will no longer be applicable. A key
strategy to make the igloo decomposition feasible would be to find a
near-optimal ratio of the partition dimensions.  This is has to be
ascertained by initial experiments.


\begin{thebibliography}{9}

\bibitem{cm5-ref}
Thinking Machines Co., 1993: CM Fortran Reference Guide {\bf Ch. 12}

\bibitem{cmmd-ref}
Thinking Machines Co., 1993: CMMD User's Guide and Reference Manual

\bibitem{mc-89}
McDonald, A.; Bates, J.R., 1989: Semi-lagrangian Integration of a Gridpoint
Shallow Water Model on the Sphere. {\em J. American Meteor. Soc.}, {\bf 117},
130-137

\bibitem{rob-81}
Robert, A.,  1981: A stable numerical integration scheme for the primitive
meteorological equations. {\em Atmos. Ocean}, {\bf 19}, 35-46

\bibitem{rob-82}
Robert, A.,  1982: A semi-Lagrangian and semi-implicit numerical integration
scheme for the primitive meteorological equations. {\em J. Meteor. Soc. Japan},
{\bf 60}, 319-324

\bibitem{stan-91}
Staniforth, A.; Cot'e, J., 1991: Semi-Lagrangian Integration Schemes for
Atmospheric Models - A Review. {\em Monthly Weather Review}, {\bf 119}, 
2206-2203


\end{thebibliography}

\newpage

\section{Appendix A: Bicubic Interpolation }


  The bicubic interpolation is used where more accuracy is desired than the 
bilinear interpolation.  \footnote{ The bilinear interpolation is more 
economical than the bicubic interpolation and is hence used where more  
accuracy is not required. } The bicubic interpolation used here is a two-step
one-dimensional lagrange polynomial (cubic) interpolation.  


The inputs to the bicubic interpolation consists of sixteen values arranged 
on a 3 cell x 3 cell grid as shown in figure 10. Using the sixteen values
 four cubic interpolations can be carried out using the four 
points available at each value of ii (in the North-South direction).  Then 
another cubic interpolation can be performed on the four resultant values in 
the East-West direction (all of which reside on a straight line which is dely 
units away from the line y = jj + 3) to obtain the desired value. 

Let \(P(y) = a_{0} + a_{1}y + a_{2}y^{2} + a_{3}y^{3} \) be a polynomial of
third degree.  The standard equation for the lagrange polynomial 
interpolation is given by, 

\begin{equation}
P(y) = \sum_{i=0}^{3}f_{i}L_{i}(y)
\end{equation}

 where,

\begin{itemize}
\item {y is the point at which interpolation takes place.}
\item {P(y) represents the desired value.}
\item {\(f_{i}\) are the known values on the grid at points \(y_{i}\).}
\item {\(L_{i}(y)\) are the  Lagrange polynomials such that 
\(L_{i}(y) =  \Pi_{k \neq i, k = 0}^{3} (y - y_{k})/(y_{i} - y_{k}) \)}

\end{itemize}

As shown in figure 10, considering the North-South direction. Let the 
position of point of interpoation be considered relative to the lower right
corner of the {\em grid square} in which it is located.  It is \(\Delta_{y} \)
units away from the lower right corner in North-South direction.

In our usage [figure 10], let

\begin{itemize}
\item { \( y_{i}\) represent jj+2, jj+3, jj+1, jj+4, \(0 \leq i \leq 3\) }
\item { \(f_{i} \) corresponds to the value of \(w_{k}\) at \(y_{i}\),
        \(i = (k-1) mod 4 \), \(1 \leq k \leq 16\)} 
\end{itemize}

As the grid points are evenly spaced, the computation of the Lagrange
polynomials can be simplified as below,

\begin{displaymath}
L_{0} = [(y - y_{1})(y - y_{2})(y - y_{3})] / 
        [(y_{0} - y_{1})(y_{0} - y_{2})(y_{0} - y_{3})]
\end{displaymath}
\begin{displaymath}
L_{0} = [(1 - \Delta_{y})(1 + \Delta_{y})(2 - \Delta_{y})] / (-1)(1)(-2)
      = 0.5 * (2.0 - \Delta_{y}) * (1.0 - \Delta_{y}^{2})
      = \Delta_{y2} * (1.0 - \Delta_{y}^{2})
\end{displaymath}
where, \(\Delta_{y2} = 0.5 * (2.0 - \Delta_{y}) \)


Similarly,
\begin{displaymath}
L_{1} = ( \Delta_{y}^2 + \Delta_{y}) * \Delta_{y2}
\end{displaymath}
\begin{displaymath}
L_{2} = ( \Delta_{y}^2 - \Delta_{y}) * \Delta_{y2} / 3
\end{displaymath}
\begin{displaymath}
L_{3} = ( \Delta_{y}^2 - 1.0) * \Delta_{y} / 6
\end{displaymath}

A similar reasoning holds for East-West direction.  Let the subscript 'y'
 be used for quantities in the North-South direction and 'x' for East-West.
The algorithm for the bicubic interpolation is as follows.

\begin {itemize}
\item Inputs
  \begin{itemize}
  \item  Sixteen values at the grid points as in figure 10
  \item  Relative positions \(\Delta_{x}\), \(\Delta_{y}\)
  \end{itemize}
\item Compute terms \( \Delta_{y2} \) and \( \Delta_{y}^{2} \) - 
2 multiplications
\item Compute terms \( \Delta_{x2} \) and \( \Delta_{x}^{2} \) - 
2 multiplications
\item Compute Lagrange polynomials \(L_{iy}\) - 6 multiplications
\item Compute Lagrange polynomials \(L_{ix}\) - 6 multiplications
\item Perform four lagrange interpolations in North-South direction -
4 x 4 multiplications
\item Perform the final lagrange interpolation in East-West direction -
4 multiplications.  This obtains the desired value in a total of 36 
multiplications.
\item Return the desired value 
\end{itemize}

\newpage

\begin{figure}
\centerline{\psfig{width=6in,height=5in,figure=original_grid.eps}} 
\end{figure} 

\begin{figure}
\centerline{\psfig{width=6in,height=5in,figure=decomp.eps}} 
\end{figure} 

\begin{figure}
\centerline{\psfig{width=6in,height=7in,figure=bilinear.eps}} 
\end{figure} 

\begin{figure}
\centerline{\psfig{width=6in,height=5in,figure=nbuf.eps}} 
\end{figure} 

\begin{figure}
\centerline{\psfig{width=6in,height=5in,figure=comm_rng.eps}} 
\end{figure} 

\begin{figure}
\centerline{\psfig{width=6in,height=5in,figure=defect.eps}} 
\end{figure} 

\begin{figure}
\centerline{\psfig{width=6in,height=5in,figure=newalg.eps}} 
\end{figure} 

\begin{figure}
\centerline{\psfig{width=6in,height=5in,figure=localinda.eps}} 
\end{figure} 

\begin{figure}
\centerline{\psfig{width=6in,height=8in,figure=localindb.eps}} 
\end{figure} 

\begin{figure}
\centerline{\psfig{width=6in,height=8in,figure=localindc.eps}} 
\end{figure} 

\begin{figure}
\centerline{\psfig{width=6in,height=5in,figure=bicubic.eps}} 
\end{figure} 

\end{document}


