C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/fizhi/AtoC.F,v 1.5 2004/07/26 18:45:17 molod Exp $
C $Name:  $

#include "FIZHI_OPTIONS.h"
       subroutine AtoC(myThid,fieldin1,fieldin2,mask,im1,im2,jm1,jm2,
     .  numlevs,Nsx,Nsy,idim1,idim2,jdim1,jdim2,fieldout1,fieldout2)
c----------------------------------------------------------------------
c  Subroutine AtoC - Routine to map a velocity component quantity
c         from the A-Grid to the C-Grid. 
c         This includes doing an exchange to fill the halo region, and 
c         then a linear average with the appropriate topography mask.
c         Also: Set up "bi, bj loop" here.
c
c  Input: myThid
c         fieldin1    Field on a-grid to move to a-grid (1st component)
c         fieldin2    Field on a-grid to move to a-grid (2nd component)
c         mask        Topography [0,1] mask - 1 to indicate above ground
c         im1,im2,jm1,jm2   Indeces in x and y for computations
c         numlevs     Number of vertical levels
c         Nsx, Nsy    
c         idim1,idim2,jdim1,jdim2   Span of fields in x and y
c
c Output: fieldout1   Field mapped to C-Grid (1st component)
c         fieldout2   Field mapped to C-Grid (2nd component)
c
c Call:   exch_uv_agrid_xyz_RL - exchange on a-grid
c-----------------------------------------------------------------------
       implicit none
#include "EEPARAMS.h"
                                                                                     
       integer myThid, numlevs
       integer Nsx, Nsy
       integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2
       _RS mask(im1:im2,jm1:jm2,numlevs,Nsx,Nsy)
       _RL fieldin1(im1:im2,jm1:jm2,numlevs,Nsx,Nsy)
       _RL fieldin2(im1:im2,jm1:jm2,numlevs,Nsx,Nsy)
       _RL fieldout1(im1:im2,jm1:jm2,numlevs,Nsx,Nsy)
       _RL fieldout2(im1:im2,jm1:jm2,numlevs,Nsx,Nsy)

       integer i, j, L, bi, bj

c Call the exchange routine to fill in the halo regions
       call exch_uv_agrid_xyz_RL (fieldin1,fieldin2,myThid)

c Now take average 
       do bj = myByLo(myThid), myByHi(myThid)
       do bi = myBxLo(myThid), myBxHi(myThid)

        do L = 1,numlevs
         do j = jdim1,jdim2
         do i = idim1,idim2
          if( (mask(i-1,j,L,bi,bj).ne.0.) .or.
     .                                  (mask(i,j,L,bi,bj).ne.0.) ) then
           fieldout1(i,j,L,bi,bj) = 
     .       ( fieldin1(i-1,j,L,bi,bj)*mask(i-1,j,L,bi,bj) +
     .         fieldin1(i,j,L,bi,bj)*mask(i,j,L,bi,bj) ) /
     .           ( mask(i-1,j,L,bi,bj) + mask(i,j,L,bi,bj) )
          else
           fieldout1(i,j,L,bi,bj) = 0.
          endif
          if( (mask(i,j-1,L,bi,bj).ne.0.) .or.
     .                                  (mask(i,j,L,bi,bj).ne.0.) ) then
           fieldout2(i,j,L,bi,bj) = 
     .       ( fieldin2(i,j-1,L,bi,bj)*mask(i,j-1,L,bi,bj) +
     .         fieldin2(i,j,L,bi,bj)*mask(i,j,L,bi,bj) ) /
     .           ( mask(i,j,L,bi,bj) + mask(i,j-1,L,bi,bj) )
          else
           fieldout2(i,j,L,bi,bj) = 0.
          endif
         enddo
         enddo
        enddo

       enddo
       enddo

       return
       end