PROGRAM start_archive2icosa

  USE xios
  USE mod_wait
  USE netcdf
  
  IMPLICIT NONE
  INCLUDE "mpif.h"
  INTEGER :: rank
  INTEGER :: size
  INTEGER :: ierr

  CHARACTER(len=*),PARAMETER :: id="client"
  INTEGER :: comm
  TYPE(xios_duration) :: dtime
  CHARACTER(len=15) :: calendar_type
  TYPE(xios_context) :: ctx_hdl

  INTEGER :: n,l
  INTEGER :: src_ibegin, src_iend, src_topo_ibegin, src_topo_iend
  INTEGER :: src_jbegin, src_jend, src_topo_jbegin, src_topo_jend
  INTEGER :: src_ni, src_ni_glo, src_topo_ni, src_topo_ni_glo
  INTEGER :: src_nj, src_nj_glo, src_topo_nj, src_topo_nj_glo
  INTEGER :: src_nlev ! number of vertical layers
  INTEGER :: src_nq=1 ! number of tracers
  INTEGER :: src_nt=1 ! number of time steps
  DOUBLE PRECISION,ALLOCATABLE :: lev_values(:) ! vertical axis 
  DOUBLE PRECISION,ALLOCATABLE :: nq_values(:) ! tracer # axis 
  DOUBLE PRECISION,ALLOCATABLE :: src_lon(:) ! mesh center coordinate
  DOUBLE PRECISION,ALLOCATABLE :: src_lat(:)
  DOUBLE PRECISION,ALLOCATABLE :: src_ap(:)
  DOUBLE PRECISION,ALLOCATABLE :: src_bp(:)
  DOUBLE PRECISION,ALLOCATABLE :: src_controle(:)
  DOUBLE PRECISION,ALLOCATABLE :: src_field_2D(:,:)
  DOUBLE PRECISION,ALLOCATABLE :: src_pk(:,:)
  DOUBLE PRECISION,ALLOCATABLE :: src_field_3D(:,:,:)
  DOUBLE PRECISION,ALLOCATABLE :: src_pressure(:,:,:)
  DOUBLE PRECISION,ALLOCATABLE :: src_theta_rhodz(:,:,:)
  DOUBLE PRECISION,ALLOCATABLE :: src_tracers(:,:,:,:)
  DOUBLE PRECISION,ALLOCATABLE :: src_topo_lon(:) ! mesh center coordinate
  DOUBLE PRECISION,ALLOCATABLE :: src_topo_lat(:)
  DOUBLE PRECISION,ALLOCATABLE :: src_topo(:,:)
  
  CHARACTER(LEN=*),PARAMETER :: src_file="start_archive_nc4.nc"
  CHARACTER(LEN=*),PARAMETER :: src_topo_file="topo_Thomas_inv_nc4.nc"
!  CHARACTER(LEN=*),PARAMETER :: src_topo_file="vt1x1inv_nc4.nc"
  CHARACTER(LEN=*),PARAMETER :: dst_coord_file="start_icosa_ref.nc"
  DOUBLE PRECISION,ALLOCATABLE :: dst_lon(:),dst_lat(:)
  DOUBLE PRECISION,ALLOCATABLE :: dst_boundslon(:,:) ! mesh corner coordinates
  DOUBLE PRECISION,ALLOCATABLE :: dst_boundslat(:,:)
  INTEGER :: dst_ibegin !, dst_iend
  INTEGER :: dst_ni, dst_ni_glo
  INTEGER :: dst_nvertex
  INTEGER :: ncid
  INTEGER :: dimids(4)
  INTEGER :: varid
  
  INTEGER :: div, remain
  INTEGER :: ts ! time step #
  DOUBLE PRECISION,PARAMETER :: pi=acos(-1.d0)
  DOUBLE PRECISION :: gravity,kappa,preff

!!! MPI Initialization
  CALL MPI_INIT(ierr)
  CALL init_wait

!!! XIOS Initialization (get the local communicator)
  CALL xios_initialize(id,return_comm=comm)
! get local rank of MPI process
  CALL MPI_COMM_RANK(comm,rank,ierr)
! get total number of MPI processes
  CALL MPI_COMM_SIZE(comm,size,ierr)

!!! Open files and load sizes and coordinates
  ierr=NF90_OPEN(src_topo_file, NF90_NOWRITE, ncid)
  ierr=NF90_INQ_VARID(ncid,"RELIEF",varid)
  ierr=NF90_INQUIRE_VARIABLE(ncid, varid,dimids=dimids)
  write(*,*) "rank=",rank,"dimids=",dimids
  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(1), len=src_topo_ni_glo)
  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(2), len=src_topo_nj_glo)
  write(*,*) "rank=",rank," src_topo_ni_glo=",src_topo_ni_glo ! longitude
  write(*,*) "rank=",rank," src_topo_nj_glo=",src_topo_nj_glo ! latitude

! assume domain splitup with MPI only along latitudes
  src_topo_ni=src_topo_ni_glo
  src_topo_ibegin=0
  src_topo_iend=src_topo_ibegin+src_ni-1
  write(*,*) "rank=",rank," src_topo_ni=",src_topo_ni
  
  src_topo_jbegin=0
  DO n=0,size-1
    src_topo_nj=src_topo_nj_glo/size
    IF (n<MOD(src_topo_nj_glo,size)) src_topo_nj=src_topo_nj+1
    IF (n==rank) exit
    src_topo_jbegin=src_topo_jbegin+src_topo_nj
  ENDDO
  src_topo_jend=src_topo_jbegin+src_topo_nj-1
  write(*,*) "rank=",rank," src_topo_nj=",src_topo_nj, &
             " src_topo_jbegin=",src_topo_jbegin

  ALLOCATE(src_topo_lon(src_topo_ni))
  ALLOCATE(src_topo_lat(src_topo_nj))
  ALLOCATE(src_topo(src_topo_ni,src_topo_nj))

! load src_topo_lon and src_topo_lat
  ierr=NF90_INQ_VARID(ncid,"longitude",varid)
  ierr=NF90_GET_VAR(ncid,varid, src_topo_lon, &
                    start=(/src_topo_ibegin+1/),count=(/src_topo_ni/))
  WRITE(*,*) rank,":src_topo_lon(1:2)=",src_topo_lon(1:2)
  ierr=NF90_INQ_VARID(ncid,"latitude",varid)
  ierr=NF90_GET_VAR(ncid,varid, src_topo_lat, &
                    start=(/src_topo_jbegin+1/),count=(/src_topo_nj/))
  WRITE(*,*) rank,":src_topo_lat(1:2)=",src_topo_lat(1:2)

! from start_archive.nc file
  ierr=NF90_OPEN(src_file, NF90_NOWRITE, ncid)
  ierr=NF90_INQ_VARID(ncid,"temp",varid)
  ierr=NF90_INQUIRE_VARIABLE(ncid, varid,dimids=dimids)
  write(*,*) "rank=",rank,"dimids=",dimids
  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(1), len=src_ni_glo)
  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(2), len=src_nj_glo)
  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(3), len=src_nlev)
  write(*,*) "rank=",rank," src_ni_glo=",src_ni_glo ! longitude
  write(*,*) "rank=",rank," src_nj_glo=",src_nj_glo ! latitude
  write(*,*) "rank=",rank," src_nlev=",src_nlev ! number of vertical layers
  write(*,*) "rank=",rank," src_nq=",src_nq ! number of tracers

! assume domain splitup with MPI only along latitudes
  src_ni=src_ni_glo
  src_ibegin=0
  src_iend=src_ibegin+src_ni-1
  write(*,*) "rank=",rank," src_ni=",src_ni
  
  src_jbegin=0
  DO n=0,size-1
    src_nj=src_nj_glo/size
    IF (n<MOD(src_nj_glo,size)) src_nj=src_nj+1
    IF (n==rank) exit
    src_jbegin=src_jbegin+src_nj
  ENDDO
  src_jend=src_jbegin+src_nj-1
  write(*,*) "rank=",rank," src_nj=",src_nj," src_jbegin=",src_jbegin

  ALLOCATE(src_lon(src_ni))
  ALLOCATE(src_lat(src_nj))
  ALLOCATE(src_field_2D(src_ni,src_nj))
  ALLOCATE(src_pk(src_ni,src_nj))
  ALLOCATE(src_field_3D(src_ni,src_nj,src_nlev))
  ALLOCATE(src_pressure(src_ni,src_nj,src_nlev+1))
  ALLOCATE(src_theta_rhodz(src_ni,src_nj,src_nlev))
  ALLOCATE(src_tracers(src_ni,src_nj,src_nlev,src_nq))

! load src_lon and src_lat
  ierr=NF90_INQ_VARID(ncid,"rlonv",varid)
  ierr=NF90_GET_VAR(ncid,varid, src_lon, &
                    start=(/src_ibegin+1/),count=(/src_ni/))
! convert rad to deg
  src_lon(1:src_ni)=src_lon(1:src_ni)*(180.d0/pi)
  WRITE(*,*) rank,":src_lon=",src_lon
  ierr=NF90_INQ_VARID(ncid,"rlatu",varid)
  ierr=NF90_GET_VAR(ncid,varid, src_lat, &
                    start=(/src_jbegin+1/),count=(/src_nj/))
! convert rad to deg
  src_lat(1:src_nj)=src_lat(1:src_nj)*(180.d0/pi)
  WRITE(*,*) rank,":src_lat=",src_lat

! load ap, bp and controle
  ALLOCATE(src_ap(src_nlev+1),src_bp(src_nlev+1),src_controle(200))
  ierr=NF90_INQ_VARID(ncid,"ap",varid)
  ierr=NF90_GET_VAR(ncid,varid,src_ap)
  WRITE(*,*) rank,":src_ap(1:5)=",src_ap(1:5)
  ierr=NF90_INQ_VARID(ncid,"bp",varid)
  ierr=NF90_GET_VAR(ncid,varid,src_bp)
  WRITE(*,*) rank,":src_bp(1:5)=",src_bp(1:5)
  ierr=NF90_INQ_VARID(ncid,"controle",varid)
  ierr=NF90_GET_VAR(ncid,varid,src_controle)
  gravity=src_controle(8)
  WRITE(*,*) rank,":gravity=",gravity
  kappa=src_controle(10)
  WRITE(*,*) rank,":kappa=",kappa
  preff=src_controle(19)
  WRITE(*,*) rank,":preff=",preff

! destination coordinates
  ierr=NF90_OPEN(dst_coord_file, NF90_NOWRITE, ncid)
  ierr=NF90_INQ_VARID(ncid,"bounds_lon_mesh",varid)
  ierr=NF90_INQUIRE_VARIABLE(ncid, varid,dimids=dimids)
  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(1), len=dst_nvertex)
  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(2), len=dst_ni_glo)
  write(*,*) "rank=",rank," dst_nvertex=",dst_nvertex ! vertex
  write(*,*) "rank=",rank," dst_ni_glo=",dst_ni_glo ! vertex boundaries

! evenly split into MPI domains
  div    = dst_ni_glo/size
  remain = MOD( dst_ni_glo, size )
  IF (rank < remain) THEN
    dst_ni=div+1 ;
    dst_ibegin=rank*(div+1) ;
  ELSE
    dst_ni=div ;
    dst_ibegin= remain * (div+1) + (rank-remain) * div ;
  ENDIF
  write(*,*) "rank=",rank," dst_ni=",dst_ni

  ALLOCATE(dst_lon(dst_ni))
  ALLOCATE(dst_lat(dst_ni))
  ALLOCATE(dst_boundslon(dst_nvertex,dst_ni))
  ALLOCATE(dst_boundslat(dst_nvertex,dst_ni))

! load dst_lon, dst_lat, dst_boundslon and dst_boundslat
  ierr=NF90_INQ_VARID(ncid,"lon_mesh",varid)
  ierr=NF90_GET_VAR(ncid,varid, dst_lon, &
                    start=(/dst_ibegin+1/), count=(/dst_ni/))
  WRITE(*,*) rank,":dst_lon(1:5)=",dst_lon(1:5)
  ierr=NF90_INQ_VARID(ncid,"lat_mesh",varid)
  ierr=NF90_GET_VAR(ncid,varid, dst_lat, &
                    start=(/dst_ibegin+1/), count=(/dst_ni/))
  WRITE(*,*) rank,":dst_lat(1:5)=",dst_lat(1:5)
  ierr=NF90_INQ_VARID(ncid,"bounds_lon_mesh",varid)
  ierr=NF90_GET_VAR(ncid,varid,dst_boundslon, &
                    start=(/1,dst_ibegin+1/), count=(/dst_nvertex,dst_ni/))
  WRITE(*,*) rank,":dst_boundslon(:,1:2)=",dst_boundslon(:,1:2)
  ierr=NF90_INQ_VARID(ncid,"bounds_lat_mesh",varid)
  ierr=NF90_GET_VAR(ncid,varid, dst_boundslat, &
                    start=(/1,dst_ibegin+1/), count=(/dst_nvertex,dst_ni/))
  WRITE(*,*) rank,":dst_boundslat(:,1:2)=",dst_boundslat(:,1:2)


! Initialize XIOS context
  WRITE(*,*) rank,":CALL xios_context_initialize()"
  CALL xios_context_initialize("test",comm)
  CALL xios_get_handle("test",ctx_hdl)
  CALL xios_set_current_context(ctx_hdl)

! Set XIOS calendar and timestep
  CALL xios_get_calendar_type(calendar_type)
  WRITE(*,*) rank,":calendar_type = ", calendar_type
  dtime%second = 3600
  CALL xios_set_timestep(dtime)

! Set axes
  ALLOCATE(lev_values(src_nlev))
  lev_values=(/ (l,l=1,src_nlev) /)
  write(*,*) rank,":lev_values()=",lev_values
  CALL xios_set_axis_attr("lev",n_glo=src_nlev,value=lev_values)
  ALLOCATE(nq_values(src_nq))
  nq_values=(/(l,l=1,src_nq)/)
  CALL xios_set_axis_attr("nq",n_glo=src_nq,value=nq_values)

! Set domains
  CALL xios_set_domain_attr("src_domain_regular", &
                            ni_glo=src_ni_glo, nj_glo=src_nj_glo, &
                            ibegin=src_ibegin, ni=src_ni, &
                            jbegin=src_jbegin, nj=src_nj, &
                            type='rectilinear')
  CALL xios_set_domain_attr("src_domain_regular", &
                             data_dim=2, &
                             data_ibegin=0, data_ni=src_ni, &
                             data_jbegin=0, data_nj=src_nj)
  CALL xios_set_domain_attr("src_domain_regular", &
                            lonvalue_1D=src_lon, &
                            latvalue_1D=src_lat)

  CALL xios_set_domain_attr("src_topo_domain_regular", &
                            ni_glo=src_topo_ni_glo, nj_glo=src_topo_nj_glo, &
                            ibegin=src_topo_ibegin, ni=src_topo_ni, &
                            jbegin=src_topo_jbegin, nj=src_topo_nj, &
                            type='rectilinear')
  CALL xios_set_domain_attr("src_topo_domain_regular", &
                             data_dim=2, &
                             data_ibegin=0, data_ni=src_topo_ni, &
                             data_jbegin=0, data_nj=src_topo_nj)
  CALL xios_set_domain_attr("src_topo_domain_regular", &
                            lonvalue_1D=src_topo_lon, &
                            latvalue_1D=src_topo_lat)

  CALL xios_set_domain_attr("src_domain_regular_clean", &
                            ni_glo=src_ni_glo-1, nj_glo=src_nj_glo, &
                            ibegin=src_ibegin, ni=src_ni-1, &
                            jbegin=src_jbegin, nj=src_nj, &
                            type='rectilinear')
  CALL xios_set_domain_attr("src_domain_regular_clean", &
                             data_dim=2, &
                             data_ibegin=0, data_ni=src_ni-1, &
                             data_jbegin=0, data_nj=src_nj)
  CALL xios_set_domain_attr("src_domain_regular_clean", &
                            lonvalue_1D=src_lon(1:src_ni-1), &
                            latvalue_1D=src_lat)

  CALL xios_set_domain_attr("dst_domain_unstructured", &
                            ni_glo=dst_ni_glo, &
                            ibegin=dst_ibegin, &
                            ni=dst_ni, &
                            type="unstructured")
  CALL xios_set_domain_attr("dst_domain_unstructured", &
                            lonvalue_1D=dst_lon, &
                            latvalue_1D=dst_lat, &
                            bounds_lon_1D=dst_boundslon, &
                            bounds_lat_1D=dst_boundslat, &
                            nvertex=dst_nvertex)

! Finalize XIOS context definition
  WRITE(*,*) rank,":CALL xios_close_context_definition()"
  CALL xios_close_context_definition()
  CALL xios_get_handle("test",ctx_hdl)
  CALL xios_set_current_context(ctx_hdl)

! Temporal loop
  DO ts=1,src_nt
    WRITE(*,*) rank,":ts=",ts
    ! Update calendar
    CALL xios_update_calendar(ts)

    ! Topography
    CALL xios_recv_field("RELIEF",src_topo)
    WRITE(*,*) rank,":topo(1:2,1:3)=",src_topo(1:2,1:3)
    ! Send surface geopotential
    CALL xios_send_field("topo",src_topo(:,:)*gravity)

    ! Surface pressure:
    !! get data using XIOS:
    CALL xios_recv_field("src_ps",src_field_2D)
    WRITE(*,*) rank,":src_ps(1:2,1:3)=",src_field_2D(1:2,1:3)
    !! write data using XIOS
    CALL xios_send_field("ps_clean",src_field_2D(1:src_ni-1,1:src_nj))

    ! compute inter-layer pressures
    DO l=1,src_nlev+1
      src_pressure(:,:,l) = src_ap(l)+src_bp(l)*src_field_2D(:,:)
    ENDDO
    
    ! surface temperature:
    CALL xios_recv_field("src_tsurf",src_field_2D)
    WRITE(*,*) rank,":src_tsurf(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("tsurf_clean",src_field_2D(1:src_ni-1,1:src_nj))
    
    ! Temperature:
    CALL xios_recv_field("src_temp",src_field_3D)
    ! compute theta_rhodz
    DO l=1,src_nlev
      src_pk(:,:)=((.5/preff)*(src_pressure(:,:,l)+src_pressure(:,:,l+1)))**kappa
      src_theta_rhodz(:,:,l) = src_field_3D(:,:,l) * &
      ((src_pressure(:,:,l)-src_pressure(:,:,l+1))/gravity)/src_pk(:,:)
    ENDDO
    CALL xios_send_field("theta_rhodz_clean", &
                         src_theta_rhodz(1:src_ni-1,1:src_nj,1:src_nlev))

    ! zonal wind
    CALL xios_recv_field("src_u",src_field_3D)
    CALL xios_send_field("u_clean", &
                          src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
    
    ! meridional wind
    CALL xios_recv_field("src_v",src_field_3D)
    CALL xios_send_field("v_clean", &
                          src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
    ! tracers
    !CALL xios_recv_field("src_tracer001",src_tracers(:,:,:,1))
    CALL xios_recv_field("src_co2",src_tracers(:,:,:,1))
    CALL xios_send_field("tracers_clean", &
                          src_tracers(1:src_ni-1,1:src_nj,1:src_nlev,1:src_nq))

    ! subsurface temperatures
    CALL xios_recv_field("src_Tsoil01",src_field_2D)
    WRITE(*,*) rank,":src_Tsoil01(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("Tsoil01_clean",src_field_2D(1:src_ni-1,1:src_nj))
    CALL xios_recv_field("src_Tsoil02",src_field_2D)
    CALL xios_send_field("Tsoil02_clean",src_field_2D(1:src_ni-1,1:src_nj))
    CALL xios_recv_field("src_Tsoil03",src_field_2D)
    CALL xios_send_field("Tsoil03_clean",src_field_2D(1:src_ni-1,1:src_nj))
    CALL xios_recv_field("src_Tsoil04",src_field_2D)
    CALL xios_send_field("Tsoil04_clean",src_field_2D(1:src_ni-1,1:src_nj))
    CALL xios_recv_field("src_Tsoil05",src_field_2D)
    CALL xios_send_field("Tsoil05_clean",src_field_2D(1:src_ni-1,1:src_nj))
    CALL xios_recv_field("src_Tsoil06",src_field_2D)
    CALL xios_send_field("Tsoil06_clean",src_field_2D(1:src_ni-1,1:src_nj))
    CALL xios_recv_field("src_Tsoil07",src_field_2D)
    CALL xios_send_field("Tsoil07_clean",src_field_2D(1:src_ni-1,1:src_nj))
    CALL xios_recv_field("src_Tsoil08",src_field_2D)
    CALL xios_send_field("Tsoil08_clean",src_field_2D(1:src_ni-1,1:src_nj))
    CALL xios_recv_field("src_Tsoil09",src_field_2D)
    CALL xios_send_field("Tsoil09_clean",src_field_2D(1:src_ni-1,1:src_nj))
    CALL xios_recv_field("src_Tsoil10",src_field_2D)
    CALL xios_send_field("Tsoil10_clean",src_field_2D(1:src_ni-1,1:src_nj))
    CALL xios_recv_field("src_Tsoil11",src_field_2D)
    CALL xios_send_field("Tsoil11_clean",src_field_2D(1:src_ni-1,1:src_nj))

    ! Albedo
    CALL xios_recv_field("src_albe",src_field_2D)
    WRITE(*,*) rank,":src_albe(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("ALBE_clean",src_field_2D(1:src_ni-1,1:src_nj))
    
    ! SW flux at the surface
    CALL xios_recv_field("src_solsw",src_field_2D)
    WRITE(*,*) rank,":src_solsw(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("solsw_clean",src_field_2D(1:src_ni-1,1:src_nj))
    
    ! LW flux at the surface
    CALL xios_recv_field("src_sollw",src_field_2D)
    WRITE(*,*) rank,":src_sollw(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("sollw_clean",src_field_2D(1:src_ni-1,1:src_nj))
    
    ! fder "Derive de flux"
    CALL xios_recv_field("src_fder",src_field_2D)
    WRITE(*,*) rank,":src_fder(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("fder_clean",src_field_2D(1:src_ni-1,1:src_nj))
    
    ! dlw "Derivee flux IR"
    CALL xios_recv_field("src_dlw",src_field_2D)
    WRITE(*,*) rank,":src_dlw(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("dlw_clean",src_field_2D(1:src_ni-1,1:src_nj))
    
    ! sollwdown "Flux IR vers le bas a la surface"
    CALL xios_recv_field("src_sollwdown",src_field_2D)
    WRITE(*,*) rank,":src_sollwdown(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("sollwdown_clean",src_field_2D(1:src_ni-1,1:src_nj))

    ! RADS "Net flux at surface"
    CALL xios_recv_field("src_RADS",src_field_2D)
    WRITE(*,*) rank,":src_RADS(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("RADS_clean",src_field_2D(1:src_ni-1,1:src_nj))

    ! ZMEA "zmea Orographie sous-maille"
    CALL xios_recv_field("src_ZMEA",src_field_2D)
    WRITE(*,*) rank,":src_ZMEA(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("ZMEA_clean",src_field_2D(1:src_ni-1,1:src_nj))
    
    ! ZSTD "zstd Orographie sous-maille"
    CALL xios_recv_field("src_ZSTD",src_field_2D)
    WRITE(*,*) rank,":src_ZSTD(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("ZSTD_clean",src_field_2D(1:src_ni-1,1:src_nj))
    
    ! ZSIG "zsig Orographie sous-maille"
    CALL xios_recv_field("src_ZSIG",src_field_2D)
    WRITE(*,*) rank,":src_ZSIG(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("ZSIG_clean",src_field_2D(1:src_ni-1,1:src_nj))
    
    ! ZGAM "zgam Orographie sous-maille"
    CALL xios_recv_field("src_ZGAM",src_field_2D)
    WRITE(*,*) rank,":src_ZGAM(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("ZGAM_clean",src_field_2D(1:src_ni-1,1:src_nj))
    
    ! ZTHE "zthe Orographie sous-maille"
    CALL xios_recv_field("src_ZTHE",src_field_2D)
    WRITE(*,*) rank,":src_ZTHE(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("ZTHE_clean",src_field_2D(1:src_ni-1,1:src_nj))
    
    ! ZPIC "zpic Orographie sous-maille"
    CALL xios_recv_field("src_ZPIC",src_field_2D)
    WRITE(*,*) rank,":src_ZPIC(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("ZPIC_clean",src_field_2D(1:src_ni-1,1:src_nj))
    
    ! ZVAL "zval Orographie sous-maille"
    CALL xios_recv_field("src_ZVAL",src_field_2D)
    WRITE(*,*) rank,":src_ZVAL(1:2,1:3)=",src_field_2D(1:2,1:3)
    CALL xios_send_field("ZVAL_clean",src_field_2D(1:src_ni-1,1:src_nj))
    
    ! TANCIEN => not necessary
  ENDDO ! of DO ts=1,src_nt
  
!! Finalize
  write(*,*) rank,":Finalize: call xios_context_finalize"
  CALL xios_context_finalize()

  write(*,*) rank,":Finalize: call MPI_COMM_FREE"
  CALL MPI_COMM_FREE(comm, ierr)

  write(*,*) rank,":Finalize: call xios_finalize"
  CALL xios_finalize()

  if (rank==0) then
    ! add a couple of things in the "startphy_icosa.nc" file
    write(*,*) rank,"Write controle() to startphy_icosa.nc"
    ierr=NF90_OPEN("startphy_icosa.nc",NF90_WRITE,ncid)
    ierr=NF90_REDEF(ncid) ! switch to define mode
    ierr=NF90_DEF_DIM(ncid,"index",100,dimids(1))
    ierr=NF90_DEF_VAR(ncid,"controle",NF90_DOUBLE,dimids(1),varid)
    ierr=NF90_ENDDEF(ncid) ! switch out of define mode
    ierr=NF90_PUT_VAR(ncid,varid,src_controle(101:200))
    if (ierr.ne.NF90_NOERR) then
      write(*,*) "NetCDF Error:",NF90_STRERROR(ierr)
    endif
    ierr=NF90_CLOSE(ncid)
    
    ! add a couple of things in the "start_icosa.nc" file
    ierr=NF90_OPEN("start_icosa.nc",NF90_WRITE,ncid)
    ierr=NF90_REDEF(ncid)
    ierr=NF90_DEF_DIM(ncid,"nvertex_u",2,dimids(1))
    ierr=NF90_DEF_VAR(ncid,"iteration",NF90_FLOAT,varid)
    ierr=NF90_ENDDEF(ncid)
    if (ierr.ne.NF90_NOERR) then
      write(*,*) "NetCDF Error:",NF90_STRERROR(ierr)
    endif
    ierr=NF90_PUT_VAR(ncid,varid,0) ! set "iteration" value to 0
    ierr=NF90_CLOSE(ncid)
    
  endif ! of if (rank==0)

  write(*,*) rank,":Finalize: call MPI_FINALIZE"
  CALL MPI_FINALIZE(ierr)

  write(*,*) rank,":my_remap: all is well that ends well!"

END PROGRAM start_archive2icosa