#define CPL_DEBUG
#define FIX_FOR_EDGE_WINDS
#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: CPL_MPMICE
C     !INTERFACE:
      SUBROUTINE CPL_MPMICE( myTime, myIter, myThid )

C     !DESCRIPTION: \bv
C     *==================================================================
C     | SUBROUTINE cpl_mpmice
C     | o Couple MITgcm ocean model with MPMice sea ice model
C     *==================================================================
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#ifdef ALLOW_EXF
# include "EXF_OPTIONS.h"
# include "EXF_FIELDS.h"
#endif
#ifdef ALLOW_SEAICE
# include "SEAICE_OPTIONS.h"
# include "SEAICE.h"
#endif

      LOGICAL  DIFFERENT_MULTIPLE
      EXTERNAL DIFFERENT_MULTIPLE

C     !LOCAL VARIABLES:
C     mytime - time counter for this thread (seconds)
C     myiter - iteration counter for this thread
C     mythid - thread number for this instance of the routine.
      _RL     mytime
      INTEGER myiter, mythid 
CEOP

#ifdef ALLOW_CPL_MPMICE
# include "EESUPPORT.h"
# include "CPL_MPMICE.h"
      COMMON /CPL_MPI_ID/
     &        myworldid, local_ocean_leader, local_ice_leader
      integer myworldid, local_ocean_leader, local_ice_leader
      integer mpistatus(MPI_STATUS_SIZE), mpierr
      integer xfer_gridsize(2)
      integer i, j, bi, bj, buffsize, idx
      Real*8  xfer_scalar
      Real*8  xfer_array(Nx,Ny)
      Real*8  xfer_bc_tracer(2*(Nx+Ny)-4)
      Real*8  xfer_bc_veloc(2*(Nx+Ny)-6)
      _RL     local(1:sNx,1:sNy,nSx,nSy)

      COMMON /FFIELDS_tmp/
     & fu_tmp, fv_tmp, Qnet_tmp, Qsw_tmp, EmPmR_tmp, saltFlux_tmp
      _RS  fu_tmp       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  fv_tmp       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  Qnet_tmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  Qsw_tmp      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  EmPmR_tmp    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)

#ifdef CPL_DEBUG
      _RL ScatArray(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          ScatArray(i,j,bi,bj) = 0.0 _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#endif

      IF( myTime .EQ. startTime ) THEN

C     Send deltatimestep
       _BEGIN_MASTER( myThid )
       IF ( myworldid .EQ. local_ocean_leader ) THEN
        xfer_scalar = deltat
        buffsize = 1
        CALL MPI_SEND(xfer_scalar,buffsize,MPI_DOUBLE_PRECISION,
     &     local_ice_leader,TimeIntervalTag,MPI_COMM_WORLD,mpierr)
#ifdef CPL_DEBUG
        print*,'MITgcm send TimeInterval', xfer_scalar
#endif
       ENDIF
       _END_MASTER( myThid )

C     Send grid dimensions (Nx,Ny)
       _BEGIN_MASTER( myThid )
       IF ( myworldid .EQ. local_ocean_leader ) THEN
        xfer_gridsize(1)=Nx
        xfer_gridsize(2)=Ny
        buffsize = 2
        CALL MPI_SEND(xfer_gridsize,buffsize,MPI_INTEGER,
     &     local_ice_leader,OceanGridsizeTag,MPI_COMM_WORLD,mpierr)
#ifdef CPL_DEBUG
        print*,'MITgcm send OceanGridsize', xfer_gridsize
#endif
       ENDIF
       _END_MASTER( myThid )

C     Send ice area
       DO bj=1,nSy
        DO bi=1,nSx
         DO j=1,sNy
          DO i=1,sNx
           local(i,j,bi,bj) = AREA(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
       CALL GATHER_2D( xfer_array, local, myThid )
       _BEGIN_MASTER( myThid )
       IF ( myworldid .EQ. local_ocean_leader ) THEN
        buffsize = Nx*Ny
        CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &     local_ice_leader,AreaTag,MPI_COMM_WORLD,mpierr)
       ENDIF
       _END_MASTER( myThid )
#ifdef CPL_DEBUG
       CALL PLOT_FIELD_XYRL( AREA, 'AREA', myIter, myThid )
#endif

C     Send ice thickness
       DO bj=1,nSy
        DO bi=1,nSx
         DO j=1,sNy
          DO i=1,sNx
           local(i,j,bi,bj) = HEFF(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
       CALL GATHER_2D( xfer_array, local, myThid )
       _BEGIN_MASTER( myThid )
       IF ( myworldid .EQ. local_ocean_leader ) THEN
        buffsize = Nx*Ny
        CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &     local_ice_leader,HeffTag,MPI_COMM_WORLD,mpierr)
       ENDIF
       _END_MASTER( myThid )
#ifdef CPL_DEBUG
       CALL PLOT_FIELD_XYRL( HEFF, 'HEFF', myIter, myThid )
#endif

C     Send ice salinity
       DO bj=1,nSy
        DO bi=1,nSx
         DO j=1,sNy
          DO i=1,sNx
           local(i,j,bi,bj) = HSALT(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
       CALL GATHER_2D( xfer_array, local, myThid )
       _BEGIN_MASTER( myThid )
       IF ( myworldid .EQ. local_ocean_leader ) THEN
        buffsize = Nx*Ny
        CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &     local_ice_leader,HsaltTag,MPI_COMM_WORLD,mpierr)
       ENDIF
       _END_MASTER( myThid )
#ifdef CPL_DEBUG
        CALL PLOT_FIELD_XYRL( HSALT, 'HSALT', myIter, myThid )
#endif

C     Send snow thickness
       DO bj=1,nSy
        DO bi=1,nSx
         DO j=1,sNy
          DO i=1,sNx
           local(i,j,bi,bj) = HSNOW(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
       CALL GATHER_2D( xfer_array, local, myThid )
       _BEGIN_MASTER( myThid )
       IF ( myworldid .EQ. local_ocean_leader ) THEN
        buffsize = Nx*Ny
        CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &     local_ice_leader,HsnowTag,MPI_COMM_WORLD,mpierr)
       ENDIF
       _END_MASTER( myThid )
#ifdef CPL_DEBUG
        CALL PLOT_FIELD_XYRL( HSNOW, 'HSNOW', myIter, myThid )
#endif

      ENDIF ! ( myTime .EQ. startTime )

C     Send ocean model time
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       xfer_scalar = myTime
       buffsize = 1
       CALL MPI_SEND(xfer_scalar,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,OceanTimeTag,MPI_COMM_WORLD,mpierr)
#ifdef CPL_DEBUG
       print*,'MITgcm send OceanTime', xfer_scalar
#endif
      ENDIF
      _END_MASTER( myThid )

C     Send boundary ice area
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = AREA(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      idx = 0
      DO i = 1, Nx
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(i,1)
      ENDDO
      DO j = 2, Ny
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(Nx,j)
      ENDDO
      DO i = (Nx-1), 1, -1
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(i,Ny)
      ENDDO
      DO j = (Ny-1), 2, -1
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(1,j)
      ENDDO
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = 2*(Nx+Ny)-4
       print*,'MITgcm is about to send AreaBcTag',buffsize
       CALL MPI_SEND(xfer_bc_tracer,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,AreaBcTag,MPI_COMM_WORLD,mpierr)
       print*,'MITgcm has sent AreaBcTag',buffsize
      ENDIF
      _END_MASTER( myThid )

C     Send boundary ice thickness
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = HEFF(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      idx = 0
      DO i = 1, Nx
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(i,1)
      ENDDO
      DO j = 2, Ny
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(Nx,j)
      ENDDO
      DO i = (Nx-1), 1, -1
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(i,Ny)
      ENDDO
      DO j = (Ny-1), 2, -1
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(1,j)
      ENDDO
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = 2*(Nx+Ny)-4
       CALL MPI_SEND(xfer_bc_tracer,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,HeffBcTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send boundary ice salinity
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = HSALT(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      idx = 0
      DO i = 1, Nx
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(i,1)
      ENDDO
      DO j = 2, Ny
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(Nx,j)
      ENDDO
      DO i = (Nx-1), 1, -1
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(i,Ny)
      ENDDO
      DO j = (Ny-1), 2, -1
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(1,j)
      ENDDO
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = 2*(Nx+Ny)-4
       CALL MPI_SEND(xfer_bc_tracer,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,HsaltBcTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send boundary snow thickness
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = HSNOW(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      idx = 0
      DO i = 1, Nx
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(i,1)
      ENDDO
      DO j = 2, Ny
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(Nx,j)
      ENDDO
      DO i = (Nx-1), 1, -1
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(i,Ny)
      ENDDO
      DO j = (Ny-1), 2, -1
       idx = idx + 1
       xfer_bc_tracer(idx) = xfer_array(1,j)
      ENDDO
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = 2*(Nx+Ny)-4
       CALL MPI_SEND(xfer_bc_tracer,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,HsnowBcTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send boundary u ice
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = UICE(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      idx = 0
      DO i = 2, Nx
       idx = idx + 1
       xfer_bc_veloc(idx) = xfer_array(i,1)
      ENDDO
      DO j = 2, Ny
       idx = idx + 1
       xfer_bc_veloc(idx) = xfer_array(Nx,j)
      ENDDO
      DO i = (Nx-1), 2, -1
       idx = idx + 1
       xfer_bc_veloc(idx) = xfer_array(i,Ny)
      ENDDO
      DO j = (Ny-1), 2, -1
       idx = idx + 1
       xfer_bc_veloc(idx) = xfer_array(2,j)
      ENDDO
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = 2*(Nx+Ny)-6
       CALL MPI_SEND(xfer_bc_veloc,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,UiceBcTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send boundary v ice
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = VICE(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      idx = 0
      DO i = 1, Nx
       idx = idx + 1
       xfer_bc_veloc(idx) = xfer_array(i,2)
      ENDDO
      DO j = 3, Ny
       idx = idx + 1
       xfer_bc_veloc(idx) = xfer_array(Nx,j)
      ENDDO
      DO i = (Nx-1), 1, -1
       idx = idx + 1
       xfer_bc_veloc(idx) = xfer_array(i,Ny)
      ENDDO
      DO j = (Ny-1), 3, -1
       idx = idx + 1
       xfer_bc_veloc(idx) = xfer_array(1,j)
      ENDDO
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = 2*(Nx+Ny)-6
       CALL MPI_SEND(xfer_bc_veloc,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,ViceBcTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send u-wind velocity 
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = uwind(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
#ifdef FIX_FOR_EDGE_WINDS
       DO j=1,Ny
        xfer_array(Nx,j)=xfer_array(Nx-1,j)
       ENDDO
#endif
       buffsize = Nx*Ny
       CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,UwindTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send v-wind velocity 
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = vwind(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
#ifdef FIX_FOR_EDGE_WINDS
       DO i=1,Nx
        xfer_array(i,Ny)=xfer_array(i,Ny-1)
       ENDDO
#endif
       buffsize = Nx*Ny
       CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,VwindTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send downward longwave radiation
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = lwdown(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,LwDownTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send downward shortwave radiation
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = swdown(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,SwDownTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send air temperature
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = atemp(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,AtempTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send humidity
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = aqh(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,AqhTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send precipitation
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = precip(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,PrecipTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send ocean surface temperature 
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = theta(i,j,1,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,SstTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send ocean surface salinity
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = salt(i,j,1,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,SssTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send surface u current
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = uVel(i,j,1,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,UvelTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )

C     Send surface v current
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          local(i,j,bi,bj) = vVel(i,j,1,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GATHER_2D( xfer_array, local, myThid )
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_SEND(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,VvelTag,MPI_COMM_WORLD,mpierr)
      ENDIF
      _END_MASTER( myThid )
#ifdef CPL_DEBUG
        CALL PLOT_FIELD_XYZRL( vVel, 'vVel(k=1)', 1, myIter, myThid )
#endif

C     Receive ice model time
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = 1
       CALL MPI_RECV(xfer_scalar,1,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,IceTimeTag,MPI_COMM_WORLD,mpistatus,mpierr)
#ifdef CPL_DEBUG
       print*,'MITgcm receive IceTime', xfer_scalar
#endif
      ENDIF
      _END_MASTER( myThid )

C     Receive ice area                     Nx*Ny       Real*8  
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_RECV(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,AreaTag,MPI_COMM_WORLD,mpistatus,mpierr)
      ENDIF
      _END_MASTER( myThid )
      CALL SCATTER_2D( xfer_array, local, myThid )
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          AREA(i,j,bi,bj) = local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
    
#ifdef CPL_DEBUG
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          ScatArray(i,j,bi,bj) = local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL PLOT_FIELD_XYRL( ScatArray, 'ice area', myIter, myThid )
#endif

C     Receive ice thickness
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_RECV(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,HeffTag,MPI_COMM_WORLD,mpistatus,mpierr)
      ENDIF
      _END_MASTER( myThid )
      CALL SCATTER_2D( xfer_array, local, myThid )
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          HEFF(i,j,bi,bj) = local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#ifdef CPL_DEBUG
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          ScatArray(i,j,bi,bj) = local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL PLOT_FIELD_XYRL( ScatArray, 'ice thickness', myIter, myThid )
#endif

C     Receive ice salinity
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_RECV(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,HsaltTag,MPI_COMM_WORLD,mpistatus,mpierr)
      ENDIF
      _END_MASTER( myThid )
      CALL SCATTER_2D( xfer_array, local, myThid )
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          HSALT(i,j,bi,bj) = local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#ifdef CPL_DEBUG
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          ScatArray(i,j,bi,bj) = local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL PLOT_FIELD_XYRL( ScatArray, 'ice salinity', myIter, myThid )
#endif

C     Receive snow thickness
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_RECV(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,HsnowTag,MPI_COMM_WORLD,mpistatus,mpierr)
      ENDIF
      _END_MASTER( myThid )
      CALL SCATTER_2D( xfer_array, local, myThid )
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          HSNOW(i,j,bi,bj) = local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#ifdef CPL_DEBUG
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          ScatArray(i,j,bi,bj) = local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL PLOT_FIELD_XYRL( ScatArray, 'ice thickness', myIter, myThid )
#endif

C     Receive u surface stress
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_RECV(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,UstressTag,MPI_COMM_WORLD,mpistatus,mpierr)
      ENDIF
      _END_MASTER( myThid )
      CALL SCATTER_2D( xfer_array, local, myThid )
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          fu(i,j,bi,bj) = AREA(i,j,bi,bj) * local(i,j,bi,bj)
     &                (1.-AREA(i,j,bi,bj) * fu_tmp(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#ifdef CPL_DEBUG
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          ScatArray(i,j,bi,bj) = local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL PLOT_FIELD_XYRL( ScatArray, 'u stress', myIter, myThid )
#endif

C     Receive v surface stress
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_RECV(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,VstressTag,MPI_COMM_WORLD,mpistatus,mpierr)
      ENDIF
      _END_MASTER( myThid )
      CALL SCATTER_2D( xfer_array, local, myThid )
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          fv(i,j,bi,bj) = AREA(i,j,bi,bj) * local(i,j,bi,bj)
     &                (1.-AREA(i,j,bi,bj) * fv_tmp(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#ifdef CPL_DEBUG
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          ScatArray(i,j,bi,bj) = local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL PLOT_FIELD_XYRL( ScatArray, 'v stress', myIter, myThid )
#endif

C     Receive residual shortwave
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_RECV(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,SwResidTag,MPI_COMM_WORLD,mpistatus,mpierr)
      ENDIF
      _END_MASTER( myThid )
      CALL SCATTER_2D( xfer_array, local, myThid )
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          Qsw(i,j,bi,bj) = -AREA(i,j,bi,bj) * local(i,j,bi,bj) +
     &                  (1.-AREA(i,j,bi,bj) * Qsw_tmp(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#ifdef CPL_DEBUG
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          ScatArray(i,j,bi,bj) = local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL PLOT_FIELD_XYRL( ScatArray, 'shortwave', myIter, myThid )
#endif

C     Receive heat flux
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_RECV(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,HeatFluxTag,MPI_COMM_WORLD,mpistatus,mpierr)
      ENDIF
      _END_MASTER( myThid )
      CALL SCATTER_2D( xfer_array, local, myThid )
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          fv(i,j,bi,bj) = Qsw(i,j,bi,bj) -
     &                   AREA(i,j,bi,bj) * local(i,j,bi,bj) +
     &               (1.-AREA(i,j,bi,bj) * Qnet_tmp(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#ifdef CPL_DEBUG
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          ScatArray(i,j,bi,bj) = local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL PLOT_FIELD_XYRL( ScatArray, 'heat flux', myIter, myThid )
#endif

C     Receive freshwater flux
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_RECV(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,WaterFluxTag,MPI_COMM_WORLD,mpistatus,mpierr)
      ENDIF
      _END_MASTER( myThid )
      CALL SCATTER_2D( xfer_array, local, myThid )
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          EmPmR(i,j,bi,bj) = - rhoConstFresh *
     &                    AREA(i,j,bi,bj) * local(i,j,bi,bj) +
     &                (1.-AREA(i,j,bi,bj) * EmPmR_tmp(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#ifdef CPL_DEBUG
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          ScatArray(i,j,bi,bj) = local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL PLOT_FIELD_XYRL( ScatArray, 'freshwater', myIter, myThid )
#endif

C     Receive salt flux
      _BEGIN_MASTER( myThid )
      IF ( myworldid .EQ. local_ocean_leader ) THEN
       buffsize = Nx*Ny
       CALL MPI_RECV(xfer_array,buffsize,MPI_DOUBLE_PRECISION,
     &    local_ice_leader,SaltFluxTag,MPI_COMM_WORLD,mpistatus,mpierr)
      ENDIF
      _END_MASTER( myThid )
      CALL SCATTER_2D( xfer_array, local, myThid )
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          saltFlux(i,j,bi,bj) = - AREA(i,j,bi,bj) * local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#ifdef CPL_DEBUG
      DO bj=1,nSy
       DO bi=1,nSx
        DO j=1,sNy
         DO i=1,sNx
          ScatArray(i,j,bi,bj) = local(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL PLOT_FIELD_XYRL( ScatArray, 'salt flux', myIter, myThid )
#endif

#endif /* ALLOW_CPL_MPMICE */

      RETURN
      END