!
! $Header$
!
MODULE cpl_mod
!
! This module excahanges and transforms all fields that should be recieved or sent to 
! coupler. The transformation of the fields are done from the grid 1D-array in phylmd 
! to the regular 2D grid accepted by the coupler. Cumulation of the fields for each 
! timestep is done in here. 
!
! Each type of surface that recevie fields from the coupler have a subroutine named 
! cpl_receive_XXX_fields and each surface that have fields to be sent to the coupler 
! have a subroutine named cpl_send_XXX_fields.
!
!*************************************************************************************

! Use statements
!*************************************************************************************
  USE dimphy, ONLY : klon
  USE mod_phys_lmdz_para
  USE ioipsl
  USE iophy

! The module oasis is always used. Without the cpp key CPP_COUPLE only the parameters 
! in the module are compiled and not the subroutines.
  USE oasis
  USE write_field_phy
  
! Global attributes
!*************************************************************************************
  IMPLICIT NONE
  PRIVATE

  ! All subroutine are public except cpl_send_all and cpl_receive_all
  PUBLIC :: cpl_init, cpl_receive_ocean_fields, cpl_receive_seaice_fields, &
       cpl_send_ocean_fields, cpl_send_seaice_fields, cpl_send_land_fields, &
       cpl_send_landice_fields, gath2cpl
  

! Declaration of module variables
!*************************************************************************************
! variable for coupling period
  INTEGER, SAVE                             :: nexca
  !$OMP THREADPRIVATE(nexca)

! variables for cumulating fields during a coupling periode :
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_sols, cpl_nsol, cpl_rain
  !$OMP THREADPRIVATE(cpl_sols,cpl_nsol,cpl_rain)
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_snow, cpl_evap, cpl_tsol
  !$OMP THREADPRIVATE(cpl_snow,cpl_evap,cpl_tsol)
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_fder, cpl_albe, cpl_taux
  !$OMP THREADPRIVATE(cpl_fder,cpl_albe,cpl_taux)
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_windsp
  !$OMP THREADPRIVATE(cpl_windsp)
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_tauy
  !$OMP THREADPRIVATE(cpl_tauy)
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_rriv2D, cpl_rcoa2D, cpl_rlic2D
  !$OMP THREADPRIVATE(cpl_rriv2D,cpl_rcoa2D,cpl_rlic2D)

! variables read from coupler :
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: read_sst     ! sea surface temperature
  !$OMP THREADPRIVATE(read_sst)
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: read_sit     ! sea ice temperature
  !$OMP THREADPRIVATE(read_sit)
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: read_sic     ! sea ice fraction
  !$OMP THREADPRIVATE(read_sic)
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: read_alb_sic ! albedo at sea ice
  !$OMP THREADPRIVATE(read_alb_sic)
  
! fraction for different surface, saved during whole coupling period
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: pctsrf_sav   
  !$OMP THREADPRIVATE(pctsrf_sav)
  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: unity
  !$OMP THREADPRIVATE(unity)
  INTEGER, SAVE                             :: nidct, nidcs
  !$OMP THREADPRIVATE(nidct,nidcs)

! variables to be sent to the coupler
  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cpl_sols2D, cpl_nsol2D, cpl_rain2D
  !$OMP THREADPRIVATE(cpl_sols2D, cpl_nsol2D, cpl_rain2D)
  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cpl_snow2D, cpl_evap2D, cpl_tsol2D
  !$OMP THREADPRIVATE(cpl_snow2D, cpl_evap2D, cpl_tsol2D)
  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cpl_fder2D, cpl_albe2D
  !$OMP THREADPRIVATE(cpl_fder2D, cpl_albe2D)
  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cpl_taux2D, cpl_tauy2D
  !$OMP THREADPRIVATE(cpl_taux2D, cpl_tauy2D)
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_windsp2D
  !$OMP THREADPRIVATE(cpl_windsp2D)
  

CONTAINS
!
!************************************************************************************
!
  SUBROUTINE cpl_init(dtime, rlon, rlat)

    INCLUDE "dimensions.h"
    INCLUDE "indicesol.h"
    INCLUDE "control.h"
    INCLUDE "temps.h"
    INCLUDE "iniprint.h"

! Input arguments
!*************************************************************************************
    REAL, INTENT(IN)                  :: dtime
    REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat

! Local variables
!*************************************************************************************
    INTEGER                           :: error, sum_error, ig, i
    INTEGER                           :: jf, nhoridct
    INTEGER                           :: nhoridcs
    INTEGER                           :: idtime
    INTEGER                           :: idayref
    INTEGER                           :: npas ! only for OASIS2
    REAL                              :: zjulian
    REAL, DIMENSION(iim,jjm+1)        :: zx_lon, zx_lat
    CHARACTER(len = 20)               :: modname = 'cpl_init'
    CHARACTER(len = 80)               :: abort_message
    CHARACTER(len=80)                 :: clintocplnam, clfromcplnam

!*************************************************************************************
! Calculate coupling period
!
!*************************************************************************************
     
    npas = itaufin/ iphysiq
    nexca = 86400 / dtime
    WRITE(lunout,*)' ##### Ocean couple #####'
    WRITE(lunout,*)' Valeurs des pas de temps'
    WRITE(lunout,*)' npas = ', npas
    WRITE(lunout,*)' nexca = ', nexca
    
!*************************************************************************************
! Allocate variables
!
!*************************************************************************************
    error = 0
    sum_error = 0

    ALLOCATE(unity(klon), stat = error)
    sum_error = sum_error + error
    ALLOCATE(pctsrf_sav(klon,nbsrf), stat = error)
    sum_error = sum_error + error
    ALLOCATE(cpl_sols(klon,2), stat = error) 
    sum_error = sum_error + error
    ALLOCATE(cpl_nsol(klon,2), stat = error)
    sum_error = sum_error + error
    ALLOCATE(cpl_rain(klon,2), stat = error)
    sum_error = sum_error + error
    ALLOCATE(cpl_snow(klon,2), stat = error)
    sum_error = sum_error + error
    ALLOCATE(cpl_evap(klon,2), stat = error)
    sum_error = sum_error + error
    ALLOCATE(cpl_tsol(klon,2), stat = error)
    sum_error = sum_error + error
    ALLOCATE(cpl_fder(klon,2), stat = error)
    sum_error = sum_error + error
    ALLOCATE(cpl_albe(klon,2), stat = error)
    sum_error = sum_error + error
    ALLOCATE(cpl_taux(klon,2), stat = error)
    sum_error = sum_error + error
    ALLOCATE(cpl_windsp(klon,2), stat = error)
    sum_error = sum_error + error
    ALLOCATE(cpl_tauy(klon,2), stat = error)
    sum_error = sum_error + error
    ALLOCATE(cpl_rriv2D(iim,jj_nb), stat=error)
    sum_error = sum_error + error
    ALLOCATE(cpl_rcoa2D(iim,jj_nb), stat=error)
    sum_error = sum_error + error
    ALLOCATE(cpl_rlic2D(iim,jj_nb), stat=error)
    sum_error = sum_error + error
    ALLOCATE(read_sst(iim, jj_nb), stat = error)
    sum_error = sum_error + error
    ALLOCATE(read_sic(iim, jj_nb), stat = error)
    sum_error = sum_error + error
    ALLOCATE(read_sit(iim, jj_nb), stat = error)
    sum_error = sum_error + error
    ALLOCATE(read_alb_sic(iim, jj_nb), stat = error)
    sum_error = sum_error + error

    IF (sum_error /= 0) THEN
       abort_message='Pb allocation variables couplees'
       CALL abort_gcm(modname,abort_message,1)
    ENDIF
!*************************************************************************************
! Initialize the allocated varaibles
!
!*************************************************************************************
    DO ig = 1, klon
       unity(ig) = ig
    ENDDO
    pctsrf_sav = 0.

    cpl_sols = 0.   ; cpl_nsol = 0.  ; cpl_rain = 0.   ; cpl_snow = 0.
    cpl_evap = 0.   ; cpl_tsol = 0.  ; cpl_fder = 0.   ; cpl_albe = 0.
    cpl_taux = 0.   ; cpl_tauy = 0.  ; cpl_rriv2D = 0. ; cpl_rcoa2D = 0.
    cpl_rlic2D = 0. ; cpl_windsp = 0.

!*************************************************************************************
! Initialize coupling
!
!*************************************************************************************
    idtime = INT(dtime)
#ifdef CPP_COUPLE
    CALL inicma
#endif

!*************************************************************************************
! initialize NetCDF output
!
!*************************************************************************************
    IF (is_sequential) THEN
       idayref = day_ini
       CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
       CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
       DO i = 1, iim
          zx_lon(i,1) = rlon(i+1)
          zx_lon(i,jjm+1) = rlon(i+1)
       ENDDO
       CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
       clintocplnam="cpl_atm_tauflx"
       CALL histbeg(clintocplnam, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),&
            1,iim,1,jjm+1, itau_phy,zjulian,dtime,nhoridct,nidct) 
! no vertical axis
       CALL histdef(nidct, 'tauxe','tauxe', &
            "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
       CALL histdef(nidct, 'tauyn','tauyn', &
            "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
       CALL histdef(nidct, 'tmp_lon','tmp_lon', &
            "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
       CALL histdef(nidct, 'tmp_lat','tmp_lat', &
            "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
       DO jf=1,jpflda2o1 + jpflda2o2
          CALL histdef(nidct, cl_writ(jf),cl_writ(jf), &
               "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
       END DO
       CALL histend(nidct)
       CALL histsync(nidct)
       
       clfromcplnam="cpl_atm_sst"
       CALL histbeg(clfromcplnam, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),1,iim,1,jjm+1, &
            0,zjulian,dtime,nhoridcs,nidcs) 
! no vertical axis
       DO jf=1,jpfldo2a
          CALL histdef(nidcs, cl_read(jf),cl_read(jf), &
               "-",iim, jjm+1, nhoridcs, 1, 1, 1, -99, 32, "inst", dtime,dtime)
       END DO
       CALL histend(nidcs)
       CALL histsync(nidcs)

    ENDIF    ! is_sequential
    
  END SUBROUTINE cpl_init
  
!
!*************************************************************************************
!

  SUBROUTINE cpl_receive_all(itime, dtime, pctsrf)
! This subroutine reads from coupler for both ocean and seaice
! 4 fields : read_sst, read_sic, read_sit and read_alb_sic. 

    INCLUDE "indicesol.h"
    INCLUDE "temps.h"
    INCLUDE "iniprint.h"
    INCLUDE "YOMCST.h"
    INCLUDE "dimensions.h"

! Input arguments
!************************************************************************************
    INTEGER, INTENT(IN)                     :: itime
    REAL, INTENT(IN)                        :: dtime
    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf

! Local variables
!************************************************************************************
    INTEGER                                 :: j, ig, il_time_secs
    INTEGER                                 :: itau_w
    INTEGER, DIMENSION(iim*(jjm+1))         :: ndexcs
    CHARACTER(len = 20)                     :: modname = 'cpl_receive_all'
    CHARACTER(len = 80)                     :: abort_message
    REAL, DIMENSION(klon)                   :: read_sic1D
    REAL, DIMENSION(iim,jj_nb,jpfldo2a)  :: tab_read_flds
    REAL, DIMENSION(iim,jj_nb)           :: read_sic

!*************************************************************************************
! Start calculation
! Get fields from coupler
!
!*************************************************************************************

#ifdef CPP_COUPLE
    il_time_secs=(itime-1)*dtime
    CALL fromcpl(il_time_secs, tab_read_flds)
#endif
    
! NetCDF output of received fields
    IF (is_sequential) THEN
       ndexcs(:) = 0
       itau_w = itau_phy + itime
       DO ig = 1, jpfldo2a
          CALL histwrite(nidcs,cl_read(ig),itau_w,tab_read_flds(:,:,ig),iim*(jjm+1),ndexcs)
       END DO
    ENDIF

! Save each field in a 2D array. 
    read_sst(:,:)     = tab_read_flds(:,:,1)  ! Sea surface temperature
    read_sic(:,:)     = tab_read_flds(:,:,2)  ! Sea ice concentration
    read_alb_sic(:,:) = tab_read_flds(:,:,3)  ! Albedo at sea ice
    read_sit(:,:)     = tab_read_flds(:,:,4)  ! Sea ice temperature

!*************************************************************************************
! Temperature and albedo are weighted with the fraction of sea-ice(read-sic)
!
!*************************************************************************************
    DO j = 1, jj_nb
       DO ig = 1, iim
          IF (ABS(1. - read_sic(ig,j)) < 0.00001) THEN
             read_sst(ig,j) = RTT - 1.8
             read_sit(ig,j) = read_sit(ig,j) / read_sic(ig,j)
             read_alb_sic(ig,j) = read_alb_sic(ig,j) / read_sic(ig,j)
          ELSE IF (ABS(read_sic(ig,j)) < 0.00001) THEN
             read_sst(ig,j) = read_sst(ig,j) / (1. - read_sic(ig,j))
             read_sit(ig,j) = read_sst(ig,j)
             read_alb_sic(ig,j) =  0.6
          ELSE
             read_sst(ig,j) = read_sst(ig,j) / (1. - read_sic(ig,j))
             read_sit(ig,j) = read_sit(ig,j) / read_sic(ig,j)
             read_alb_sic(ig,j) = read_alb_sic(ig,j) / read_sic(ig,j)
          ENDIF
       ENDDO
    ENDDO
    
!*************************************************************************************
!  Transform seaice fraction, read_sic into pctsrf_sav
!
!*************************************************************************************
    CALL cpl2gath(read_sic, read_sic1D, klon, unity)

    DO ig = 1, klon
       ! treatment only of ocean and/or seaice points
       IF (pctsrf(ig,is_oce) > epsfra .OR. &
            pctsrf(ig,is_sic) > epsfra) THEN
          pctsrf_sav(ig,is_sic) = (pctsrf(ig,is_oce) + pctsrf(ig,is_sic)) &
               * read_sic1D(ig)
          pctsrf_sav(ig,is_oce) = (pctsrf(ig,is_oce) + pctsrf(ig,is_sic)) &
               - pctsrf_sav(ig,is_sic)
       ENDIF
    ENDDO

!*************************************************************************************
! To avoid round up problems
!
!*************************************************************************************
    WHERE (ABS(pctsrf_sav(:,is_sic)) .LE. 2.*EPSILON(pctsrf_sav(1,is_sic)))
       pctsrf_sav(:,is_sic) = 0.
       pctsrf_sav(:,is_oce) = pctsrf(:,is_oce) + pctsrf(:,is_sic)
    ENDWHERE
    WHERE (ABS(pctsrf_sav(:,is_oce)) .LE. 2.*EPSILON(pctsrf_sav(1,is_oce)))
       pctsrf_sav(:,is_sic) = pctsrf(:,is_oce) + pctsrf(:,is_sic)
       pctsrf_sav(:,is_oce) = 0.
    ENDWHERE
    IF (MINVAL(pctsrf_sav(:,is_oce)) < 0.) THEN
       WRITE(*,*)'Pb fraction ocean inferieure a 0'
       WRITE(*,*)'au point ',MINLOC(pctsrf_sav(:,is_oce))
       WRITE(*,*)'valeur = ',MINVAL(pctsrf_sav(:,is_oce)) 
       abort_message = 'voir ci-dessus'
       CALL abort_gcm(modname,abort_message,1)
    ENDIF
    IF (MINVAL(pctsrf_sav(:,is_sic)) < 0.) THEN
       WRITE(*,*)'Pb fraction glace inferieure a 0'
       WRITE(*,*)'au point ',MINLOC(pctsrf_sav(:,is_sic))
       WRITE(*,*)'valeur = ',MINVAL(pctsrf_sav(:,is_sic)) 
       abort_message = 'voir ci-dessus'
       CALL abort_gcm(modname,abort_message,1)
    ENDIF
        
  END SUBROUTINE cpl_receive_all
!
!*************************************************************************************
! 
  SUBROUTINE cpl_receive_ocean_fields(itime, dtime, knon, knindex, pctsrf, &
       tsurf_new, pctsrf_oce)
!
! This routine reads, if first time step in a coupling period, all fields reveived from
! coupler for all types of surfaces. It returns the fields for the ocean surface which
! are the sea surface temperature and the fraction of ocean. 
! The fields are transformed into 1D arrays with valid points : 
! tsurf_new(1:knon), pctsrf(1:klon).
!
    INCLUDE "indicesol.h"

! Input arguments
!*************************************************************************************
    INTEGER, INTENT(IN)                     :: itime
    REAL, INTENT(IN)                        :: dtime
    INTEGER, INTENT(IN)                     :: knon
    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf

! Output arguments
!*************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)      :: tsurf_new
    REAL, DIMENSION(klon), INTENT(OUT)      :: pctsrf_oce

!*************************************************************************************
! If first time step in a coupling period receive all fields for all types
! of surfaces from coupler : read_sst, read_sit, read_alb_sic and pctsrf_sav. 
! 
!*************************************************************************************

    IF (MOD(itime, nexca) == 1) CALL cpl_receive_all(itime, dtime, pctsrf)


!*************************************************************************************
! Transform read_sst into compressed 1D variable tsurf_new
!
!*************************************************************************************
    CALL cpl2gath(read_sst, tsurf_new, knon, knindex)
    pctsrf_oce(:) = pctsrf_sav(:,is_oce) 


  END SUBROUTINE cpl_receive_ocean_fields
!
!*************************************************************************************
!
  SUBROUTINE cpl_receive_seaice_fields(knon, knindex, &
       tsurf_new, alb_new, pctsrf_sic)
!
! This routine returns the fields for the seaice that have been read from the coupler
! (done earlier with cpl_receive_ocean_fields). These fields are the temperature and 
! albedo at sea ice surface and fraction of sea ice.
! The fields are transformed into  1D arrays with valid points : 
! tsurf_new(1:knon), alb_new(1:knon), pctsrf(1:klon).
!
    INCLUDE "indicesol.h"

! Input arguments
!*************************************************************************************
    INTEGER, INTENT(IN)                     :: knon
    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex

! Output arguments
!*************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)      :: tsurf_new
    REAL, DIMENSION(klon), INTENT(OUT)      :: alb_new
    REAL, DIMENSION(klon), INTENT(OUT)      :: pctsrf_sic


!*************************************************************************************
! Transform fields read from coupler from 2D into compressed 1D variables
!
!*************************************************************************************
    CALL cpl2gath(read_sit, tsurf_new, knon, knindex)
    CALL cpl2gath(read_alb_sic, alb_new, knon, knindex)
    pctsrf_sic(:) = pctsrf_sav(:,is_sic) 

  END SUBROUTINE cpl_receive_seaice_fields

!
!*************************************************************************************
!

  SUBROUTINE cpl_send_ocean_fields(itime, knon, knindex, &
       swdown, lwdown, fluxlat, fluxsens, &
       precip_rain, precip_snow, evap, tsurf, fder, albsol, taux, tauy, windsp)
!
! This subroutine cumulates some fields for each time-step during a coupling 
! period. At last time-step in a coupling period the fields are transformed to the 
! grid accepted by the coupler. No sending to the coupler will be done from here 
! (it is done in cpl_send_seaice_fields).
!
    INCLUDE "indicesol.h"
    INCLUDE "dimensions.h"

! Input arguments
!*************************************************************************************
    INTEGER, INTENT(IN)                     :: itime
    INTEGER, INTENT(IN)                     :: knon
    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
    REAL, DIMENSION(klon), INTENT(IN)       :: swdown, lwdown 
    REAL, DIMENSION(klon), INTENT(IN)       :: fluxlat, fluxsens
    REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow
    REAL, DIMENSION(klon), INTENT(IN)       :: evap, tsurf, fder, albsol
    REAL, DIMENSION(klon), INTENT(IN)       :: taux, tauy, windsp

! Local variables
!*************************************************************************************
    INTEGER                                 :: cpl_index, ig 
    INTEGER                                 :: error, sum_error
    CHARACTER(len = 25)                     :: modname = 'cpl_send_ocean_fields'
    CHARACTER(len = 80)                     :: abort_message

!*************************************************************************************
! Start calculation
! The ocean points are saved with second array index=1
!
!*************************************************************************************
    cpl_index = 1

!*************************************************************************************
! Reset fields to zero in the beginning of a new coupling period 
!
!*************************************************************************************
    IF (MOD(itime, nexca) == 1) THEN
       cpl_sols(:,cpl_index) = 0.0
       cpl_nsol(:,cpl_index) = 0.0
       cpl_rain(:,cpl_index) = 0.0
       cpl_snow(:,cpl_index) = 0.0
       cpl_evap(:,cpl_index) = 0.0
       cpl_tsol(:,cpl_index) = 0.0
       cpl_fder(:,cpl_index) = 0.0
       cpl_albe(:,cpl_index) = 0.0
       cpl_taux(:,cpl_index) = 0.0
       cpl_tauy(:,cpl_index) = 0.0
       cpl_windsp(:,cpl_index) = 0.0
    ENDIF
       
!*************************************************************************************
! Cumulate at each time-step
!
!*************************************************************************************    
    DO ig = 1, knon
       cpl_sols(ig,cpl_index) = cpl_sols(ig,cpl_index) + &
            swdown(ig)      / FLOAT(nexca)
       cpl_nsol(ig,cpl_index) = cpl_nsol(ig,cpl_index) + &
            (lwdown(ig) + fluxlat(ig) +fluxsens(ig)) / FLOAT(nexca)
       cpl_rain(ig,cpl_index) = cpl_rain(ig,cpl_index) + &
            precip_rain(ig) / FLOAT(nexca)
       cpl_snow(ig,cpl_index) = cpl_snow(ig,cpl_index) + &
            precip_snow(ig) / FLOAT(nexca)
       cpl_evap(ig,cpl_index) = cpl_evap(ig,cpl_index) + &
            evap(ig)        / FLOAT(nexca)
       cpl_tsol(ig,cpl_index) = cpl_tsol(ig,cpl_index) + &
            tsurf(ig)       / FLOAT(nexca)
       cpl_fder(ig,cpl_index) = cpl_fder(ig,cpl_index) + &
            fder(ig)        / FLOAT(nexca)
       cpl_albe(ig,cpl_index) = cpl_albe(ig,cpl_index) + &
            albsol(ig)      / FLOAT(nexca)
       cpl_taux(ig,cpl_index) = cpl_taux(ig,cpl_index) + &
            taux(ig)        / FLOAT(nexca)
       cpl_tauy(ig,cpl_index) = cpl_tauy(ig,cpl_index) + &
            tauy(ig)        / FLOAT(nexca)      
       cpl_windsp(ig,cpl_index) = cpl_windsp(ig,cpl_index) + &
            windsp(ig)      / FLOAT(nexca) 
    ENDDO

!*************************************************************************************
! If the time-step corresponds to the end of coupling period the 
! fields are transformed to the 2D grid. 
! No sending to the coupler (it is done from cpl_send_seaice_fields).
!
!*************************************************************************************
    IF (MOD(itime, nexca) == 0) THEN

       IF (.NOT. ALLOCATED(cpl_sols2D)) THEN
          sum_error = 0
          ALLOCATE(cpl_sols2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_nsol2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_rain2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_snow2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_evap2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_tsol2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_fder2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_albe2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_taux2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_tauy2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_windsp2D(iim,jj_nb), stat=error)
          sum_error = sum_error + error
          
          IF (sum_error /= 0) THEN
             abort_message='Pb allocation variables couplees pour l''ecriture'
             CALL abort_gcm(modname,abort_message,1)
          ENDIF
       ENDIF
       

       CALL gath2cpl(cpl_sols(1,cpl_index), cpl_sols2D(1,1,cpl_index), &
            knon, knindex)

       CALL gath2cpl(cpl_nsol(1,cpl_index), cpl_nsol2D(1,1,cpl_index), &
            knon, knindex)

       CALL gath2cpl(cpl_rain(1,cpl_index), cpl_rain2D(1,1,cpl_index), &
            knon, knindex)

       CALL gath2cpl(cpl_snow(1,cpl_index), cpl_snow2D(1,1,cpl_index), &
            knon, knindex)

       CALL gath2cpl(cpl_evap(1,cpl_index), cpl_evap2D(1,1,cpl_index), &
            knon, knindex)

! cpl_tsol2D(:,:,:) not used!
       CALL gath2cpl(cpl_tsol(1,cpl_index), cpl_tsol2D(1,1, cpl_index), &
            knon, knindex)

! cpl_fder2D(:,:,1) not used, only cpl_fder(:,:,2)!
       CALL gath2cpl(cpl_fder(1,cpl_index), cpl_fder2D(1,1,cpl_index), &
            knon, knindex)

! cpl_albe2D(:,:,:) not used!
       CALL gath2cpl(cpl_albe(1,cpl_index), cpl_albe2D(1,1,cpl_index), &
            knon, knindex)

       CALL gath2cpl(cpl_taux(1,cpl_index), cpl_taux2D(1,1,cpl_index), &
            knon, knindex)

       CALL gath2cpl(cpl_tauy(1,cpl_index), cpl_tauy2D(1,1,cpl_index), &
            knon, knindex)

       CALL gath2cpl(cpl_windsp(1,cpl_index), cpl_windsp2D(1,1), &
            knon, knindex)

    ENDIF

  END SUBROUTINE cpl_send_ocean_fields

!
!*************************************************************************************
!

  SUBROUTINE cpl_send_seaice_fields(itime, dtime, knon, knindex, &
       pctsrf, lafin, rlon, rlat, &
       swdown, lwdown, fluxlat, fluxsens, &
       precip_rain, precip_snow, evap, tsurf, fder, albsol, taux, tauy)
!
! This subroutine cumulates some fields for each time-step during a coupling 
! period. At last time-step in a coupling period the fields are transformed to the 
! grid accepted by the coupler. All fields for all types of surfaces are sent to
! the coupler.
!
    INCLUDE "indicesol.h"
    INCLUDE "dimensions.h"

! Input arguments
!*************************************************************************************
    INTEGER, INTENT(IN)                     :: itime
    INTEGER, INTENT(IN)                     :: knon
    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
    REAL, INTENT(IN)                        :: dtime
    REAL, DIMENSION(klon), INTENT(IN)       :: rlon, rlat
    REAL, DIMENSION(klon), INTENT(IN)       :: swdown, lwdown 
    REAL, DIMENSION(klon), INTENT(IN)       :: fluxlat, fluxsens
    REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow
    REAL, DIMENSION(klon), INTENT(IN)       :: evap, tsurf, fder
    REAL, DIMENSION(klon), INTENT(IN)       :: albsol, taux, tauy
    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
    LOGICAL, INTENT(IN)                     :: lafin

! Local variables
!*************************************************************************************
    INTEGER                                 :: cpl_index, ig 
    INTEGER                                 :: error, sum_error
    CHARACTER(len = 25)                     :: modname = 'cpl_send_seaice_fields'
    CHARACTER(len = 80)                     :: abort_message


!*************************************************************************************
! Start calulation
! The sea-ice points are saved with second array index=2
!
!*************************************************************************************
    cpl_index = 2

!*************************************************************************************
! Reset fields to zero in the beginning of a new coupling period 
!
!*************************************************************************************
    IF (MOD(itime, nexca) == 1) THEN
       cpl_sols(:,cpl_index) = 0.0
       cpl_nsol(:,cpl_index) = 0.0
       cpl_rain(:,cpl_index) = 0.0
       cpl_snow(:,cpl_index) = 0.0
       cpl_evap(:,cpl_index) = 0.0
       cpl_tsol(:,cpl_index) = 0.0
       cpl_fder(:,cpl_index) = 0.0
       cpl_albe(:,cpl_index) = 0.0
       cpl_taux(:,cpl_index) = 0.0
       cpl_tauy(:,cpl_index) = 0.0
    ENDIF
       
!*************************************************************************************
! Cumulate at each time-step
!
!*************************************************************************************    
    DO ig = 1, knon
       cpl_sols(ig,cpl_index) = cpl_sols(ig,cpl_index) + &
            swdown(ig)      / FLOAT(nexca)
       cpl_nsol(ig,cpl_index) = cpl_nsol(ig,cpl_index) + &
            (lwdown(ig) + fluxlat(ig) +fluxsens(ig)) / FLOAT(nexca)
       cpl_rain(ig,cpl_index) = cpl_rain(ig,cpl_index) + &
            precip_rain(ig) / FLOAT(nexca)
       cpl_snow(ig,cpl_index) = cpl_snow(ig,cpl_index) + &
            precip_snow(ig) / FLOAT(nexca)
       cpl_evap(ig,cpl_index) = cpl_evap(ig,cpl_index) + &
            evap(ig)        / FLOAT(nexca)
       cpl_tsol(ig,cpl_index) = cpl_tsol(ig,cpl_index) + &
            tsurf(ig)       / FLOAT(nexca)
       cpl_fder(ig,cpl_index) = cpl_fder(ig,cpl_index) + &
            fder(ig)        / FLOAT(nexca)
       cpl_albe(ig,cpl_index) = cpl_albe(ig,cpl_index) + &
            albsol(ig)      / FLOAT(nexca)
       cpl_taux(ig,cpl_index) = cpl_taux(ig,cpl_index) + &
            taux(ig)        / FLOAT(nexca)
       cpl_tauy(ig,cpl_index) = cpl_tauy(ig,cpl_index) + &
            tauy(ig)        / FLOAT(nexca)      
    ENDDO

!*************************************************************************************
! If the time-step corresponds to the end of coupling period the 
! fields are transformed to the 2D grid and all fields are sent to coupler.
!
!*************************************************************************************
    IF (MOD(itime, nexca) == 0) THEN
       IF (.NOT. ALLOCATED(cpl_sols2D)) THEN
          sum_error = 0
          ALLOCATE(cpl_sols2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_nsol2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_rain2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_snow2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_evap2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_tsol2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_fder2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_albe2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_taux2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_tauy2D(iim,jj_nb,2), stat=error)
          sum_error = sum_error + error
          ALLOCATE(cpl_windsp2D(iim,jj_nb), stat=error)
          sum_error = sum_error + error
          
          IF (sum_error /= 0) THEN
             abort_message='Pb allocation variables couplees pour l''ecriture'
             CALL abort_gcm(modname,abort_message,1)
          ENDIF
       ENDIF

       CALL gath2cpl(cpl_sols(1,cpl_index), cpl_sols2D(1,1,cpl_index), &
            knon, knindex)

       CALL gath2cpl(cpl_nsol(1,cpl_index), cpl_nsol2D(1,1,cpl_index), &
            knon, knindex)

       CALL gath2cpl(cpl_rain(1,cpl_index), cpl_rain2D(1,1,cpl_index), &
            knon, knindex)

       CALL gath2cpl(cpl_snow(1,cpl_index), cpl_snow2D(1,1,cpl_index), &
            knon, knindex)

       CALL gath2cpl(cpl_evap(1,cpl_index), cpl_evap2D(1,1,cpl_index), &
            knon, knindex)

! cpl_tsol2D(:,:,:) not used!
       CALL gath2cpl(cpl_tsol(1,cpl_index), cpl_tsol2D(1,1, cpl_index), &
            knon, knindex)

       CALL gath2cpl(cpl_fder(1,cpl_index), cpl_fder2D(1,1,cpl_index), &
            knon, knindex)

! cpl_albe2D(:,:,:) not used!
       CALL gath2cpl(cpl_albe(1,cpl_index), cpl_albe2D(1,1,cpl_index), &
            knon, knindex)

       CALL gath2cpl(cpl_taux(1,cpl_index), cpl_taux2D(1,1,cpl_index), &
            knon, knindex)

       CALL gath2cpl(cpl_tauy(1,cpl_index), cpl_tauy2D(1,1,cpl_index), &
            knon, knindex)

       ! Send all fields
       CALL cpl_send_all(itime, dtime, pctsrf, lafin, rlon, rlat)
    ENDIF

  END SUBROUTINE cpl_send_seaice_fields

!
!*************************************************************************************
!

  SUBROUTINE cpl_send_land_fields(itime, knon, knindex, rriv_in, rcoa_in)
!
! This subroutine cumulates some fields for each time-step during a coupling 
! period. At last time-step in a coupling period the fields are transformed to the 
! grid accepted by the coupler. No sending to the coupler will be done from here 
! (it is done in cpl_send_seaice_fields).
!
    INCLUDE "dimensions.h"

! Input arguments
!*************************************************************************************
    INTEGER, INTENT(IN)                       :: itime
    INTEGER, INTENT(IN)                       :: knon
    INTEGER, DIMENSION(klon), INTENT(IN)      :: knindex
    REAL, DIMENSION(klon), INTENT(IN)         :: rriv_in
    REAL, DIMENSION(klon), INTENT(IN)         :: rcoa_in

! Local variables
!*************************************************************************************
    REAL, DIMENSION(iim,jj_nb)             :: rriv2D
    REAL, DIMENSION(iim,jj_nb)             :: rcoa2D

!*************************************************************************************
! Rearrange fields in 2D variables 
! First initialize to zero to avoid unvalid points causing problems
!
!*************************************************************************************
    rriv2D(:,:) = 0.0
    rcoa2D(:,:) = 0.0
    CALL gath2cpl(rriv_in, rriv2D, knon, knindex)
    CALL gath2cpl(rcoa_in, rcoa2D, knon, knindex)

!*************************************************************************************
! Reset cumulated fields to zero in the beginning of a new coupling period 
!
!*************************************************************************************
    IF (MOD(itime, nexca) == 1) THEN
       cpl_rriv2D(:,:) = 0.0
       cpl_rcoa2D(:,:) = 0.0
    ENDIF

!*************************************************************************************
! Cumulate : Following fields should be cumulated at each time-step
!
!*************************************************************************************    
    cpl_rriv2D(:,:) = cpl_rriv2D(:,:) + rriv2D(:,:) / FLOAT(nexca)
    cpl_rcoa2D(:,:) = cpl_rcoa2D(:,:) + rcoa2D(:,:) / FLOAT(nexca)

  END SUBROUTINE cpl_send_land_fields

!
!*************************************************************************************
!

  SUBROUTINE cpl_send_landice_fields(itime, knon, knindex, rlic_in)
! This subroutine cumulates the field for melting ice for each time-step 
! during a coupling period. This routine will not send to coupler. Sending 
! will be done in cpl_send_seaice_fields.
!
    INCLUDE "dimensions.h"
    
! Input varibales
!*************************************************************************************
    INTEGER, INTENT(IN)                       :: itime
    INTEGER, INTENT(IN)                       :: knon
    INTEGER, DIMENSION(klon), INTENT(IN)      :: knindex
    REAL, DIMENSION(klon), INTENT(IN)         :: rlic_in

! Local varibales
!*************************************************************************************
    REAL, DIMENSION(iim,jj_nb)             :: rlic2D

!*************************************************************************************
! Rearrange field in a 2D variable 
! First initialize to zero to avoid unvalid points causing problems
!
!*************************************************************************************
    rlic2D(:,:) = 0.0
    CALL gath2cpl(rlic_in, rlic2D, knon, knindex)

!*************************************************************************************
! Reset field to zero in the beginning of a new coupling period 
!
!*************************************************************************************
    IF (MOD(itime, nexca) == 1) THEN
       cpl_rlic2D(:,:) = 0.0
    ENDIF

!*************************************************************************************
! Cumulate : Melting ice should be cumulated at each time-step
!
!*************************************************************************************    
    cpl_rlic2D(:,:) = cpl_rlic2D(:,:) + rlic2D(:,:) / FLOAT(nexca)

  END SUBROUTINE cpl_send_landice_fields

!
!*************************************************************************************
!

  SUBROUTINE cpl_send_all(itime, dtime, pctsrf, lafin, rlon, rlat)
! This routine will send fields for all different surfaces to the coupler.
! This subroutine should be executed after calculations by the last surface(sea-ice),
! all calculations at the different surfaces have to be done before. 
!    
! Some includes
!*************************************************************************************
    INCLUDE "indicesol.h"
    INCLUDE "temps.h"
    INCLUDE "dimensions.h"
    
! Input arguments
!*************************************************************************************
    INTEGER, INTENT(IN)                                  :: itime
    REAL, INTENT(IN)                                     :: dtime
    REAL, DIMENSION(klon), INTENT(IN)                    :: rlon, rlat
    REAL, DIMENSION(klon,nbsrf), INTENT(IN)              :: pctsrf
    LOGICAL, INTENT(IN)                                  :: lafin
    
! Local variables
!*************************************************************************************
    INTEGER                                              :: error, sum_error, j
    INTEGER                                              :: itau_w
    INTEGER                                              :: il_time_secs
    INTEGER, DIMENSION(iim*(jjm+1))                      :: ndexct
    REAL                                                 :: Up, Down
    REAL, DIMENSION(iim, jj_nb)                          :: tmp_lon, tmp_lat
    REAL, DIMENSION(iim, jj_nb, 4)                       :: pctsrf2D
    REAL, DIMENSION(iim, jj_nb)                          :: deno
    CHARACTER(len = 20)                                  :: modname = 'cpl_send_all'
    CHARACTER(len = 80)                                  :: abort_message
   
! Variables with fields to coupler
    REAL, DIMENSION(iim, jj_nb)                          :: tmp_taux
    REAL, DIMENSION(iim, jj_nb)                          :: tmp_tauy
    REAL, DIMENSION(iim, jj_nb)                          :: tmp_calv
! Table with all fields to send to coupler
    REAL, DIMENSION(iim, jj_nb, jpflda2o1+jpflda2o2)     :: tab_flds
#ifdef CPP_PARA
    INCLUDE 'mpif.h'
    INTEGER, DIMENSION(MPI_STATUS_SIZE)                  :: status
#endif

! End definitions
!*************************************************************************************
    


!*************************************************************************************
! All fields are stored in a table tab_flds(:,:,:)
! First store the fields 7 to 18 which are already on the right format
!
!*************************************************************************************
    tab_flds(:,:,7)  = cpl_windsp2D(:,:)
    tab_flds(:,:,8)  = cpl_sols2D(:,:,2)
    tab_flds(:,:,9)  = cpl_sols2D(:,:,1)
    tab_flds(:,:,10) = cpl_nsol2D(:,:,2)
    tab_flds(:,:,11) = cpl_nsol2D(:,:,1)
    tab_flds(:,:,12) = cpl_fder2D(:,:,2)
    tab_flds(:,:,13) = cpl_evap2D(:,:,2)
    tab_flds(:,:,14) = cpl_evap2D(:,:,1)
    tab_flds(:,:,17) = cpl_rcoa2D(:,:)
    tab_flds(:,:,18) = cpl_rriv2D(:,:)
    
!*************************************************************************************
! Transform the fraction of sub-surfaces from 1D to 2D array
!
!*************************************************************************************
    pctsrf2D(:,:,:) = 0.
    CALL gath2cpl(pctsrf(:,is_oce), pctsrf2D(:,:,is_oce), klon, unity)
    CALL gath2cpl(pctsrf(:,is_sic), pctsrf2D(:,:,is_sic), klon, unity)
    CALL gath2cpl(pctsrf(:,is_lic), pctsrf2D(:,:,is_lic), klon, unity)

!*************************************************************************************
! Calculate the average calving per latitude
! Store calving in tab_flds(:,:,19)
! 
!*************************************************************************************      
    DO j = 1, jj_nb
       tmp_calv(:,j) = DOT_PRODUCT (cpl_rlic2D(1:iim,j), &
            pctsrf2D(1:iim,j,is_lic)) / REAL(iim)
    ENDDO
    
    
    IF (is_parallel) THEN
       IF (.NOT. is_north_pole) THEN
#ifdef CPP_PARA
          CALL MPI_RECV(Up,1,MPI_REAL_LMDZ,mpi_rank-1,1234,COMM_LMDZ_PHY,status,error)
          CALL MPI_SEND(tmp_calv(1,1),1,MPI_REAL_LMDZ,mpi_rank-1,1234,COMM_LMDZ_PHY,error)
#endif
       ENDIF
       
       IF (.NOT. is_south_pole) THEN
#ifdef CPP_PARA
          CALL MPI_SEND(tmp_calv(1,jj_nb),1,MPI_REAL_LMDZ,mpi_rank+1,1234,COMM_LMDZ_PHY,error)
          CALL MPI_RECV(down,1,MPI_REAL_LMDZ,mpi_rank+1,1234,COMM_LMDZ_PHY,status,error)
#endif
       ENDIF
       
       IF (.NOT. is_north_pole .AND. ii_begin /=1) THEN
          Up=Up+tmp_calv(iim,1)
          tmp_calv(:,1)=Up
       ENDIF
       
       IF (.NOT. is_south_pole .AND. ii_end /= iim) THEN
          Down=Down+tmp_calv(1,jj_nb)
          tmp_calv(:,jj_nb)=Down	 
       ENDIF
    ENDIF
      

    tab_flds(:,:,19) = tmp_calv(:,:)

!*************************************************************************************
! Calculate total flux for snow, rain and wind with weighted addition using the 
! fractions of ocean and seaice.
!
! Store the fields for rain and snow directly in tab_flds(:,:,15) and 
! tab_flds(:,:,16) respectively.
!
!*************************************************************************************    
    tab_flds(:,:,15) = 0.0
    tab_flds(:,:,16) = 0.0
    tmp_taux(:,:)    = 0.0
    tmp_tauy(:,:)    = 0.0
    

    ! fraction oce+seaice
    deno =  pctsrf2D(:,:,is_oce) + pctsrf2D(:,:,is_sic) 
    ! For all valid grid cells containing some fraction of ocean or sea-ice
    WHERE ( deno(:,:) /= 0 )
       tab_flds(:,:,15) = cpl_rain2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
            cpl_rain2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
       tab_flds(:,:,16) = cpl_snow2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
            cpl_snow2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
       
       tmp_taux = cpl_taux2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
            cpl_taux2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
       tmp_tauy = cpl_tauy2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
            cpl_tauy2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
    ENDWHERE

!*************************************************************************************
! Transform the wind components from local atmospheric 2D coordinates to geocentric 
! 3D coordinates. 
! Store the resulting wind components in tab_flds(:,:,1:6)
!*************************************************************************************

! Transform the longitudes and latitudes on 2D arrays
    CALL Grid1DTo2D_mpi(rlon,tmp_lon)
    CALL Grid1DTo2D_mpi(rlat,tmp_lat)
    
    IF (is_sequential) THEN
       IF (is_north_pole) tmp_lon(:,1)     = tmp_lon(:,2)
       IF (is_south_pole) tmp_lon(:,jjm+1) = tmp_lon(:,jjm)
    ENDIF
      
! NetCDF output of the wind before transformation of coordinate system
    IF (is_sequential) THEN
       ndexct(:) = 0
       itau_w = itau_phy + itime
       CALL histwrite(nidct,'tauxe',itau_w,tmp_taux,iim*(jjm+1),ndexct)
       CALL histwrite(nidct,'tauyn',itau_w,tmp_tauy,iim*(jjm+1),ndexct)
       CALL histwrite(nidct,'tmp_lon',itau_w,tmp_lon,iim*(jjm+1),ndexct)
       CALL histwrite(nidct,'tmp_lat',itau_w,tmp_lat,iim*(jjm+1),ndexct)
    ENDIF

! Transform the wind from local atmospheric 2D coordinates to geocentric 
! 3D coordinates 
    CALL atm2geo (iim, jj_nb, tmp_taux, tmp_tauy, tmp_lon, tmp_lat, &
         tab_flds(:,:,1), tab_flds(:,:,2), tab_flds(:,:,3) )
    
    tab_flds(:,:,4)  = tab_flds(:,:,1)
    tab_flds(:,:,5)  = tab_flds(:,:,2)
    tab_flds(:,:,6)  = tab_flds(:,:,3)

!*************************************************************************************
! NetCDF output of all fields just before sending to coupler.
!
!*************************************************************************************
    IF (is_sequential) THEN
       CALL histwrite(nidct,cl_writ(8), itau_w,tab_flds(:,:,8), iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(9), itau_w,tab_flds(:,:,9), iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(10),itau_w,tab_flds(:,:,10),iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(11),itau_w,tab_flds(:,:,11),iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(12),itau_w,tab_flds(:,:,12),iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(13),itau_w,tab_flds(:,:,13),iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(14),itau_w,tab_flds(:,:,14),iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(15),itau_w,tab_flds(:,:,15),iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(16),itau_w,tab_flds(:,:,16),iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(17),itau_w,tab_flds(:,:,17),iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(18),itau_w,tab_flds(:,:,18),iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(19),itau_w,tab_flds(:,:,19),iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(1), itau_w,tab_flds(:,:,1), iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(2), itau_w,tab_flds(:,:,2), iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(3), itau_w,tab_flds(:,:,3), iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(4), itau_w,tab_flds(:,:,4), iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(5), itau_w,tab_flds(:,:,5), iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(6), itau_w,tab_flds(:,:,6), iim*(jjm+1),ndexct)
       CALL histwrite(nidct,cl_writ(7), itau_w,tab_flds(:,:,7), iim*(jjm+1),ndexct)
       CALL histsync(nidct)
    ENDIF


!*************************************************************************************
! Send the table of all fields
!
!*************************************************************************************
#ifdef CPP_COUPLE
    il_time_secs=(itime-1)*dtime
    CALL intocpl(il_time_secs, lafin, tab_flds(:,:,:))
#endif

!*************************************************************************************
! Finish with some dellocate
!
!*************************************************************************************  
    sum_error=0
    DEALLOCATE(cpl_sols2D, cpl_nsol2D, cpl_rain2D, cpl_snow2D, stat=error )
    sum_error = sum_error + error
    DEALLOCATE(cpl_evap2D, cpl_tsol2D, cpl_fder2D, cpl_albe2D, stat=error )
    sum_error = sum_error + error
    DEALLOCATE(cpl_taux2D, cpl_tauy2D, cpl_windsp2D, stat=error )
    sum_error = sum_error + error
    IF (sum_error /= 0) THEN
       abort_message='Pb in deallocation of cpl_xxxx2D coupling variables'
       CALL abort_gcm(modname,abort_message,1)
    ENDIF
    
  END SUBROUTINE cpl_send_all
!
!*************************************************************************************
!
  SUBROUTINE cpl2gath(champ_in, champ_out, knon, knindex)
! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
! au coupleur.
!
! 
! input:         
!   champ_in     champ sur la grille gathere        
!   knon         nombre de points dans le domaine a traiter
!   knindex      index des points de la surface a traiter
!
! output:
!   champ_out    champ sur la grille 2D
!
    INCLUDE "dimensions.h"

! Input
    INTEGER, INTENT(IN)                       :: knon
    REAL, DIMENSION(iim,jj_nb), INTENT(IN)    :: champ_in
    INTEGER, DIMENSION(klon), INTENT(IN)      :: knindex

! Output
    REAL, DIMENSION(klon), INTENT(OUT)        :: champ_out

! Local
    INTEGER                                   :: i, ig
    REAL, DIMENSION(klon)                     :: tamp

!*************************************************************************************
!
! Transform from 2 dimensions (iim,jj_nb) to 1 dimension (klon)
    CALL Grid2Dto1D_mpi(champ_in,tamp)

! Compress from klon to knon
    DO i = 1, knon
       ig = knindex(i)
       champ_out(i) = tamp(ig)
    ENDDO

  END SUBROUTINE cpl2gath
!
!*************************************************************************************
!
  SUBROUTINE gath2cpl(champ_in, champ_out, knon, knindex)
! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
! au coupleur.
!
! input:         
!   champ_in     champ sur la grille gathere        
!   knon         nombre de points dans le domaine a traiter
!   knindex      index des points de la surface a traiter
!
! output:
!   champ_out    champ sur la grille 2D
!
    INCLUDE "dimensions.h"
    
! Input arguments
!*************************************************************************************
    INTEGER, INTENT(IN)                    :: knon
    REAL, DIMENSION(klon), INTENT(IN)      :: champ_in
    INTEGER, DIMENSION(klon), INTENT(IN)   :: knindex

! Output arguments
!*************************************************************************************
    REAL, DIMENSION(iim,jj_nb), INTENT(OUT) :: champ_out

! Local variables
!*************************************************************************************
    INTEGER                                :: i, ig
    REAL, DIMENSION(klon)                  :: tamp

!*************************************************************************************

! Decompress from knon to klon
    tamp = 0.
    DO i = 1, knon
       ig = knindex(i)
       tamp(ig) = champ_in(i)
    ENDDO

! Transform from 1 dimension (klon) to 2 dimensions (iim,jj_nb)
    CALL Grid1Dto2D_mpi(tamp,champ_out)
    
    IF (is_north_pole) champ_out(:,1)=tamp(1)
    IF (is_south_pole) champ_out(:,jj_nb)=tamp(klon)
    
  END SUBROUTINE gath2cpl
!
!*************************************************************************************
!
END MODULE cpl_mod

