!
! $Header$
!
MODULE iophy

  USE phys_output_var_mod

! abd  REAL,private,allocatable,DIMENSION(:),save :: io_lat
! abd  REAL,private,allocatable,DIMENSION(:),save :: io_lon
  REAL,ALLOCATABLE,DIMENSION(:),SAVE :: io_lat
  REAL,ALLOCATABLE,DIMENSION(:),SAVE :: io_lon
  INTEGER, SAVE :: phys_domain_id
  INTEGER, SAVE :: npstn
  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nptabij
  INTEGER, SAVE :: itau_iophy

!$OMP THREADPRIVATE(itau_iophy)
  
  INTERFACE histwrite_phy
    MODULE PROCEDURE histwrite2d_phy,histwrite3d_phy,histwrite2d_phy_old,histwrite3d_phy_old
  END INTERFACE

  INTERFACE histbeg_phy_all
    MODULE PROCEDURE histbeg_phy,histbeg_phy_points
  END INTERFACE


CONTAINS

! ug Routine pour définir itau_iophy depuis phys_output_write_mod:
  SUBROUTINE set_itau_iophy(ito)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: ito
      itau_iophy = ito
  END SUBROUTINE

  SUBROUTINE init_iophy_new(rlat,rlon)
  USE dimphy
  USE mod_phys_lmdz_para
  USE mod_grid_phy_lmdz
  USE ioipsl
  IMPLICIT NONE
  INCLUDE 'dimensions.h'   
    REAL,DIMENSION(klon),INTENT(IN) :: rlon
    REAL,DIMENSION(klon),INTENT(IN) :: rlat

    REAL,DIMENSION(klon_glo)        :: rlat_glo
    REAL,DIMENSION(klon_glo)        :: rlon_glo
    
    INTEGER,DIMENSION(2) :: ddid
    INTEGER,DIMENSION(2) :: dsg
    INTEGER,DIMENSION(2) :: dsl
    INTEGER,DIMENSION(2) :: dpf
    INTEGER,DIMENSION(2) :: dpl
    INTEGER,DIMENSION(2) :: dhs
    INTEGER,DIMENSION(2) :: dhe 
    INTEGER :: i    

    CALL gather(rlat,rlat_glo)
    CALL bcast(rlat_glo)
    CALL gather(rlon,rlon_glo)
    CALL bcast(rlon_glo)
    
!$OMP MASTER  
    ALLOCATE(io_lat(jjm+1-1/(iim*jjm)))
    io_lat(1)=rlat_glo(1)
    io_lat(jjm+1-1/(iim*jjm))=rlat_glo(klon_glo)
    IF ((iim*jjm) > 1) then
      DO i=2,jjm
        io_lat(i)=rlat_glo(2+(i-2)*iim)
      ENDDO
    ENDIF

    ALLOCATE(io_lon(iim))
    io_lon(:)=rlon_glo(2-1/(iim*jjm):iim+1-1/(iim*jjm))

    ddid=(/ 1,2 /)
    dsg=(/ iim, jjm+1-1/(iim*jjm) /)
    dsl=(/ iim, jj_nb /)
    dpf=(/ 1,jj_begin /)
    dpl=(/ iim, jj_end /)
    dhs=(/ ii_begin-1,0 /)
    IF (mpi_rank==mpi_size-1) THEN
      dhe=(/0,0/)
    ELSE
      dhe=(/ iim-ii_end,0 /)  
    ENDIF
    
    CALL flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe, &
                      'APPLE',phys_domain_id)

!$OMP END MASTER
      
  END SUBROUTINE init_iophy_new

  SUBROUTINE init_iophy(lat,lon)
  USE dimphy
  USE mod_phys_lmdz_para
  USE ioipsl
  IMPLICIT NONE
  INCLUDE 'dimensions.h'   
    REAL,DIMENSION(iim),INTENT(IN) :: lon
    REAL,DIMENSION(jjm+1-1/(iim*jjm)),INTENT(IN) :: lat

    INTEGER,DIMENSION(2) :: ddid
    INTEGER,DIMENSION(2) :: dsg
    INTEGER,DIMENSION(2) :: dsl
    INTEGER,DIMENSION(2) :: dpf
    INTEGER,DIMENSION(2) :: dpl
    INTEGER,DIMENSION(2) :: dhs
    INTEGER,DIMENSION(2) :: dhe 

!$OMP MASTER  
    allocate(io_lat(jjm+1-1/(iim*jjm)))
    io_lat(:)=lat(:)
    allocate(io_lon(iim))
    io_lon(:)=lon(:)
   
    ddid=(/ 1,2 /)
    dsg=(/ iim, jjm+1-1/(iim*jjm) /)
    dsl=(/ iim, jj_nb /)
    dpf=(/ 1,jj_begin /)
    dpl=(/ iim, jj_end /)
    dhs=(/ ii_begin-1,0 /)
    if (mpi_rank==mpi_size-1) then
      dhe=(/0,0/)
    else
      dhe=(/ iim-ii_end,0 /)  
    endif
    
    call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe, &
                      'APPLE',phys_domain_id)

!$OMP END MASTER
      
  end SUBROUTINE init_iophy
  
  SUBROUTINE histbeg_phy(name,itau0,zjulian,dtime,nhori,nid_day)
  USE dimphy
  USE mod_phys_lmdz_para
  use ioipsl
  use write_field
  IMPLICIT NONE
  include 'dimensions.h'
    
    character*(*), INTENT(IN) :: name
    integer, INTENT(IN) :: itau0
    REAL,INTENT(IN) :: zjulian
    REAL,INTENT(IN) :: dtime
    integer,intent(out) :: nhori
    integer,intent(out) :: nid_day

!$OMP MASTER    
    if (is_sequential) then
      call histbeg(name,iim,io_lon, jj_nb,io_lat(jj_begin:jj_end), &
                   1,iim,1,jj_nb,itau0, zjulian, dtime, nhori, nid_day)
    else
      call histbeg(name,iim,io_lon, jj_nb,io_lat(jj_begin:jj_end), &
                   1,iim,1,jj_nb,itau0, zjulian, dtime, nhori, nid_day,phys_domain_id)
    endif
!$OMP END MASTER
  
  END SUBROUTINE histbeg_phy

  SUBROUTINE histbeg_phy_points(rlon,rlat,pim,tabij,ipt,jpt, &
             plon,plat,plon_bounds,plat_bounds, &
             nname,itau0,zjulian,dtime,nnhori,nnid_day)
  USE dimphy
  USE mod_phys_lmdz_para
  USE mod_grid_phy_lmdz
  use ioipsl
  use write_field
  IMPLICIT NONE
  include 'dimensions.h' 

    REAL,DIMENSION(klon),INTENT(IN) :: rlon
    REAL,DIMENSION(klon),INTENT(IN) :: rlat
    integer, INTENT(IN) :: itau0
    REAL,INTENT(IN) :: zjulian
    REAL,INTENT(IN) :: dtime
    integer, INTENT(IN) :: pim
    integer, intent(out) :: nnhori
    character(len=20), INTENT(IN) :: nname
    INTEGER, intent(out) :: nnid_day
    integer :: i
    REAL,DIMENSION(klon_glo)        :: rlat_glo
    REAL,DIMENSION(klon_glo)        :: rlon_glo
    INTEGER, DIMENSION(pim), INTENT(IN)  :: tabij
    REAL,DIMENSION(pim), INTENT(IN) :: plat, plon
    INTEGER,DIMENSION(pim), INTENT(IN) :: ipt, jpt
    REAL,DIMENSION(pim,2), intent(out) :: plat_bounds, plon_bounds

    INTEGER, SAVE :: tabprocbeg, tabprocend
!$OMP THREADPRIVATE(tabprocbeg, tabprocend)
    INTEGER :: ip
    INTEGER, PARAMETER :: nip=1
    INTEGER :: npproc
    REAL, allocatable, DIMENSION(:) :: npplat, npplon
    REAL, allocatable, DIMENSION(:,:) :: npplat_bounds, npplon_bounds
    INTEGER, PARAMETER :: jjmp1=jjm+1-1/jjm
    REAL, DIMENSION(iim,jjmp1) :: zx_lon, zx_lat

    CALL gather(rlat,rlat_glo)
    CALL bcast(rlat_glo)
    CALL gather(rlon,rlon_glo)
    CALL bcast(rlon_glo)

!$OMP MASTER
    DO i=1,pim

!    print*,'CFMIP_iophy i tabij lon lat',i,tabij(i),plon(i),plat(i)

     plon_bounds(i,1)=rlon_glo(tabij(i)-1)
     plon_bounds(i,2)=rlon_glo(tabij(i)+1)
     if(plon_bounds(i,2).LE.0..AND.plon_bounds(i,1).GE.0.) THEN
      if(rlon_glo(tabij(i)).GE.0.) THEN
       plon_bounds(i,2)=-1*plon_bounds(i,2)
      endif
     endif
     if(plon_bounds(i,2).GE.0..AND.plon_bounds(i,1).LE.0.) THEN
      if(rlon_glo(tabij(i)).LE.0.) THEN
       plon_bounds(i,2)=-1*plon_bounds(i,2)
      endif
     endif
!
     IF ( tabij(i).LE.iim) THEN
      plat_bounds(i,1)=rlat_glo(tabij(i))
     ELSE
      plat_bounds(i,1)=rlat_glo(tabij(i)-iim)
     ENDIF
     plat_bounds(i,2)=rlat_glo(tabij(i)+iim)
!
!    print*,'CFMIP_iophy point i lon lon_bds',i,plon_bounds(i,1),rlon_glo(tabij(i)),plon_bounds(i,2) 
!    print*,'CFMIP_iophy point i lat lat_bds',i,plat_bounds(i,1),rlat_glo(tabij(i)),plat_bounds(i,2) 
!
    ENDDO
    if (is_sequential) then

     npstn=pim
     IF(.NOT. ALLOCATED(nptabij)) THEN
      ALLOCATE(nptabij(pim))
     ENDIF 
     DO i=1,pim
      nptabij(i)=tabij(i)
     ENDDO

       CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlon_glo,zx_lon)
       if ((iim*jjm).gt.1) then
       DO i = 1, iim
         zx_lon(i,1) = rlon_glo(i+1)
         zx_lon(i,jjmp1) = rlon_glo(i+1)
       ENDDO
       endif
       CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlat_glo,zx_lat)

    DO i=1,pim
!    print*,'CFMIP_iophy i tabij lon lat',i,tabij(i),plon(i),plat(i)

     plon_bounds(i,1)=zx_lon(ipt(i)-1,jpt(i))
     plon_bounds(i,2)=zx_lon(ipt(i)+1,jpt(i))

     if (ipt(i).EQ.1) then
      plon_bounds(i,1)=zx_lon(iim,jpt(i))
      plon_bounds(i,2)=360.+zx_lon(ipt(i)+1,jpt(i))
     endif
 
     if (ipt(i).EQ.iim) then
      plon_bounds(i,2)=360.+zx_lon(1,jpt(i))
     endif

     plat_bounds(i,1)=zx_lat(ipt(i),jpt(i)-1)
     plat_bounds(i,2)=zx_lat(ipt(i),jpt(i)+1)

     if (jpt(i).EQ.1) then
      plat_bounds(i,1)=zx_lat(ipt(i),1)+0.001
      plat_bounds(i,2)=zx_lat(ipt(i),1)-0.001
     endif
 
     if (jpt(i).EQ.jjmp1) then
      plat_bounds(i,1)=zx_lat(ipt(i),jjmp1)+0.001
      plat_bounds(i,2)=zx_lat(ipt(i),jjmp1)-0.001
     endif
!
!    print*,'CFMIP_iophy point i lon lon_bds',i,plon_bounds(i,1),rlon(tabij(i)),plon_bounds(i,2) 
!    print*,'CFMIP_iophy point i lat lat_bds',i,plat_bounds(i,1),rlat(tabij(i)),plat_bounds(i,2) 
!
    ENDDO
!    print*,'iophy is_sequential nname, nnhori, nnid_day=',trim(nname),nnhori,nnid_day
     call histbeg(nname,pim,plon,plon_bounds, & 
                           plat,plat_bounds, &
                           itau0, zjulian, dtime, nnhori, nnid_day)
    else
     npproc=0
     DO ip=1, pim
      tabprocbeg=klon_mpi_begin
      tabprocend=klon_mpi_end
      IF(tabij(ip).GE.tabprocbeg.AND.tabij(ip).LE.tabprocend) THEN
       npproc=npproc+1
       npstn=npproc
      ENDIF 
     ENDDO
!    print*,'CFMIP_iophy mpi_rank npstn',mpi_rank,npstn
     IF(.NOT. ALLOCATED(nptabij)) THEN
      ALLOCATE(nptabij(npstn))
      ALLOCATE(npplon(npstn), npplat(npstn))
      ALLOCATE(npplon_bounds(npstn,2), npplat_bounds(npstn,2))
     ENDIF
     npproc=0
     DO ip=1, pim
      IF(tabij(ip).GE.tabprocbeg.AND.tabij(ip).LE.tabprocend) THEN
       npproc=npproc+1
       nptabij(npproc)=tabij(ip)
!      print*,'mpi_rank npproc ip plon plat tabij=',mpi_rank,npproc,ip, &
!      plon(ip),plat(ip),tabij(ip)
       npplon(npproc)=plon(ip)
       npplat(npproc)=plat(ip)
       npplon_bounds(npproc,1)=plon_bounds(ip,1)
       npplon_bounds(npproc,2)=plon_bounds(ip,2)
       npplat_bounds(npproc,1)=plat_bounds(ip,1)
       npplat_bounds(npproc,2)=plat_bounds(ip,2)
!!!
!!! print qui sert a reordonner les points stations selon l'ordre CFMIP
!!! ne pas enlever
        print*,'iophy_mpi rank ip lon lat',mpi_rank,ip,plon(ip),plat(ip)
!!!
      ENDIF
     ENDDO
     call histbeg(nname,npstn,npplon,npplon_bounds, &
                            npplat,npplat_bounds, &
                            itau0,zjulian,dtime,nnhori,nnid_day,phys_domain_id)
    endif
!$OMP END MASTER

  end SUBROUTINE histbeg_phy_points


  SUBROUTINE histdef2d_old (iff,lpoint,flag_var,nomvar,titrevar,unitvar)

    USE ioipsl
    USE dimphy
    USE mod_phys_lmdz_para

    IMPLICIT NONE

    INCLUDE "dimensions.h"
    INCLUDE "temps.h"
    INCLUDE "clesphys.h"

    INTEGER                          :: iff
    LOGICAL                          :: lpoint
    INTEGER, DIMENSION(nfiles)       :: flag_var
    CHARACTER(LEN=20)                 :: nomvar
    CHARACTER(LEN=*)                 :: titrevar
    CHARACTER(LEN=*)                 :: unitvar

    REAL zstophym

    IF (type_ecri(iff)=='inst(X)'.OR.type_ecri(iff)=='once') THEN
       zstophym=zoutm(iff)
    ELSE
       zstophym=zdtime_moy
    ENDIF

    ! Appel a la lecture des noms et niveau d'ecriture des variables dans output.def
    CALL conf_physoutputs(nomvar,flag_var)

    IF(.NOT.lpoint) THEN  
       IF ( flag_var(iff)<=lev_files(iff) ) THEN
          CALL histdef (nid_files(iff),nomvar,titrevar,unitvar, &
               iim,jj_nb,nhorim(iff), 1,1,1, -99, 32, &
               type_ecri(iff), zstophym,zoutm(iff))                
       ENDIF
    ELSE
       IF ( flag_var(iff)<=lev_files(iff) ) THEN
          CALL histdef (nid_files(iff),nomvar,titrevar,unitvar, &
               npstn,1,nhorim(iff), 1,1,1, -99, 32, &
               type_ecri(iff), zstophym,zoutm(iff))                
       ENDIF
    ENDIF

    ! Set swaero_diag=true if at least one of the concerned variables are defined 
    IF (nomvar=='topswad' .OR. nomvar=='topswai' .OR. nomvar=='solswad' .OR. nomvar=='solswai' ) THEN
       IF  ( flag_var(iff)<=lev_files(iff) ) THEN
          swaero_diag=.TRUE.
       END IF
    END IF
  END SUBROUTINE histdef2d_old



  SUBROUTINE histdef3d_old (iff,lpoint,flag_var,nomvar,titrevar,unitvar)

    USE ioipsl
    USE dimphy
    USE mod_phys_lmdz_para

    IMPLICIT NONE

    INCLUDE "dimensions.h"
    INCLUDE "temps.h"
!    INCLUDE "indicesol.h"
    INCLUDE "clesphys.h"

    INTEGER                          :: iff
    LOGICAL                          :: lpoint
    INTEGER, DIMENSION(nfiles)       :: flag_var
    CHARACTER(LEN=20)                 :: nomvar
    CHARACTER(LEN=*)                 :: titrevar
    CHARACTER(LEN=*)                 :: unitvar

    REAL zstophym

    ! Appel a la lecture des noms et niveau d'ecriture des variables dans output.def
    CALL conf_physoutputs(nomvar,flag_var)

    IF (type_ecri(iff)=='inst(X)'.OR.type_ecri(iff)=='once') THEN
       zstophym=zoutm(iff)
    ELSE
       zstophym=zdtime_moy
    ENDIF

    IF(.NOT.lpoint) THEN
       IF ( flag_var(iff)<=lev_files(iff) ) THEN
          CALL histdef (nid_files(iff), nomvar, titrevar, unitvar, &
               iim, jj_nb, nhorim(iff), klev, levmin(iff), &
               levmax(iff)-levmin(iff)+1, nvertm(iff), 32, type_ecri(iff), &
               zstophym, zoutm(iff))
       ENDIF
    ELSE
       IF ( flag_var(iff)<=lev_files(iff) ) THEN
          CALL histdef (nid_files(iff), nomvar, titrevar, unitvar, &
               npstn,1,nhorim(iff), klev, levmin(iff), &
               levmax(iff)-levmin(iff)+1, nvertm(iff), 32, &
               type_ecri(iff), zstophym,zoutm(iff))
       ENDIF
    ENDIF
  END SUBROUTINE histdef3d_old








  SUBROUTINE histdef2d (iff,var)

    USE ioipsl
    USE dimphy
    USE mod_phys_lmdz_para

    IMPLICIT NONE

    INCLUDE "dimensions.h"
    INCLUDE "temps.h"
    INCLUDE "clesphys.h"

    INTEGER                          :: iff
    TYPE(ctrl_out)                   :: var

    REAL zstophym
    CHARACTER(LEN=20) :: typeecrit


    ! ug On récupère le type écrit de la structure:
    !       Assez moche, à refaire si meilleure méthode...
    IF (INDEX(var%type_ecrit(iff), "once") > 0) THEN
       typeecrit = 'once'
    ELSE IF(INDEX(var%type_ecrit(iff), "t_min") > 0) THEN
       typeecrit = 't_min(X)'
    ELSE IF(INDEX(var%type_ecrit(iff), "t_max") > 0) THEN
       typeecrit = 't_max(X)'
    ELSE IF(INDEX(var%type_ecrit(iff), "inst") > 0) THEN
       typeecrit = 'inst(X)'
    ELSE
       typeecrit = type_ecri_files(iff)
    ENDIF

    IF (typeecrit=='inst(X)'.OR.typeecrit=='once') THEN
       zstophym=zoutm(iff)
    ELSE
       zstophym=zdtime_moy
    ENDIF

    ! Appel a la lecture des noms et niveau d'ecriture des variables dans output.def
    CALL conf_physoutputs(var%name, var%flag)

    IF(.NOT.clef_stations(iff)) THEN  
       IF ( var%flag(iff)<=lev_files(iff) ) THEN
          CALL histdef (nid_files(iff), var%name, var%description, var%unit, &
               iim,jj_nb,nhorim(iff), 1,1,1, -99, 32, &
               typeecrit, zstophym,zoutm(iff))                
       ENDIF
    ELSE
       IF ( var%flag(iff)<=lev_files(iff)) THEN
          CALL histdef (nid_files(iff), var%name, var%description, var%unit, &
               npstn,1,nhorim(iff), 1,1,1, -99, 32, &
               typeecrit, zstophym,zoutm(iff))                
       ENDIF
    ENDIF

    ! Set swaero_diag=true if at least one of the concerned variables are defined 
    IF (var%name=='topswad' .OR. var%name=='topswai' .OR. var%name=='solswad' .OR. var%name=='solswai' ) THEN
       IF  ( var%flag(iff)<=lev_files(iff) ) THEN
          swaero_diag=.TRUE.
       END IF
    END IF
  END SUBROUTINE histdef2d
  SUBROUTINE histdef3d (iff,var)

    USE ioipsl
    USE dimphy
    USE mod_phys_lmdz_para

    IMPLICIT NONE

    INCLUDE "dimensions.h"
    INCLUDE "temps.h"
    INCLUDE "clesphys.h"

    INTEGER                          :: iff
    TYPE(ctrl_out)                   :: var

    REAL zstophym
    CHARACTER(LEN=20) :: typeecrit

    ! ug On récupère le type écrit de la structure:
    !       Assez moche, à refaire si meilleure méthode...
    IF (INDEX(var%type_ecrit(iff), "once") > 0) THEN
       typeecrit = 'once'
    ELSE IF(INDEX(var%type_ecrit(iff), "t_min") > 0) THEN
       typeecrit = 't_min(X)'
    ELSE IF(INDEX(var%type_ecrit(iff), "t_max") > 0) THEN
       typeecrit = 't_max(X)'
    ELSE IF(INDEX(var%type_ecrit(iff), "inst") > 0) THEN
       typeecrit = 'inst(X)'
    ELSE
       typeecrit = type_ecri_files(iff)
    ENDIF


    ! Appel a la lecture des noms et niveau d'ecriture des variables dans output.def
    CALL conf_physoutputs(var%name,var%flag)

    IF (typeecrit=='inst(X)'.OR.typeecrit=='once') THEN
       zstophym=zoutm(iff)
    ELSE
       zstophym=zdtime_moy
    ENDIF

    IF(.NOT.clef_stations(iff)) THEN
       IF ( var%flag(iff)<=lev_files(iff) ) THEN
          CALL histdef (nid_files(iff), var%name, var%description, var%unit, &
               iim, jj_nb, nhorim(iff), klev, levmin(iff), &
               levmax(iff)-levmin(iff)+1, nvertm(iff), 32, typeecrit, &
               zstophym, zoutm(iff))
       ENDIF
    ELSE
       IF ( var%flag(iff)<=lev_files(iff)) THEN
          CALL histdef (nid_files(iff), var%name, var%description, var%unit, &
               npstn,1,nhorim(iff), klev, levmin(iff), &
               levmax(iff)-levmin(iff)+1, nvertm(iff), 32, &
               typeecrit, zstophym,zoutm(iff))
       ENDIF
    ENDIF
  END SUBROUTINE histdef3d

  SUBROUTINE conf_physoutputs(nam_var,flag_var)
!!! Lecture des noms et niveau de sortie des variables dans output.def
    !   en utilisant les routines getin de IOIPSL  
    use ioipsl

    IMPLICIT NONE

    include 'iniprint.h'

    CHARACTER(LEN=20)                :: nam_var
    INTEGER, DIMENSION(nfiles)      :: flag_var

    IF(prt_level>10) WRITE(lunout,*)'Avant getin: nam_var flag_var ',nam_var,flag_var(:)
    CALL getin('flag_'//nam_var,flag_var)
    CALL getin('name_'//nam_var,nam_var)
    IF(prt_level>10) WRITE(lunout,*)'Apres getin: nam_var flag_var ',nam_var,flag_var(:)

  END SUBROUTINE conf_physoutputs


 
  SUBROUTINE histwrite2d_phy_old(nid,lpoint,name,itau,field)
  USE dimphy
  USE mod_phys_lmdz_para
  USE ioipsl
  IMPLICIT NONE
  include 'dimensions.h'
  include 'iniprint.h'
    
    integer,INTENT(IN) :: nid
    logical,INTENT(IN) :: lpoint 
    character*(*), INTENT(IN) :: name
    integer, INTENT(IN) :: itau
    REAL,DIMENSION(:),INTENT(IN) :: field
    REAL,DIMENSION(klon_mpi) :: buffer_omp
    INTEGER, allocatable, DIMENSION(:) :: index2d
    REAL :: Field2d(iim,jj_nb)

    integer :: ip
    REAL,allocatable,DIMENSION(:) :: fieldok


    IF (size(field)/=klon) CALL abort_gcm('iophy::histwrite2d','Field first DIMENSION not equal to klon',1)
    
    CALL Gather_omp(field,buffer_omp)    
!$OMP MASTER
    CALL grid1Dto2D_mpi(buffer_omp,Field2d)
    if(.NOT.lpoint) THEN
     ALLOCATE(index2d(iim*jj_nb))
     ALLOCATE(fieldok(iim*jj_nb))
     IF (prt_level >= 9) write(lunout,*)'Sending ',name,' to IOIPSL'
     CALL histwrite(nid,name,itau,Field2d,iim*jj_nb,index2d)
     IF (prt_level >= 9) write(lunout,*)'Finished sending ',name,' to IOIPSL'
    else
     ALLOCATE(fieldok(npstn))
     ALLOCATE(index2d(npstn))

     if(is_sequential) then
!     klon_mpi_begin=1
!     klon_mpi_end=klon
      DO ip=1, npstn
       fieldok(ip)=buffer_omp(nptabij(ip))
      ENDDO
     else
      DO ip=1, npstn
!     print*,'histwrite2d is_sequential npstn ip name nptabij',npstn,ip,name,nptabij(ip)
       IF(nptabij(ip).GE.klon_mpi_begin.AND. &
          nptabij(ip).LE.klon_mpi_end) THEN
         fieldok(ip)=buffer_omp(nptabij(ip)-klon_mpi_begin+1)
       ENDIF
      ENDDO
     endif
     IF (prt_level >= 9) write(lunout,*)'Sending ',name,' to IOIPSL'
     CALL histwrite(nid,name,itau,fieldok,npstn,index2d)
     IF (prt_level >= 9) write(lunout,*)'Finished sending ',name,' to IOIPSL'
!
    endif
    deallocate(index2d)
    deallocate(fieldok)
!$OMP END MASTER   

 
  end SUBROUTINE histwrite2d_phy_old

  SUBROUTINE histwrite3d_phy_old(nid,lpoint,name,itau,field)
  USE dimphy
  USE mod_phys_lmdz_para

  use ioipsl
  IMPLICIT NONE
  include 'dimensions.h'
  include 'iniprint.h'
    
    integer,INTENT(IN) :: nid
    logical,INTENT(IN) :: lpoint
    character*(*), INTENT(IN) :: name
    integer, INTENT(IN) :: itau
    REAL,DIMENSION(:,:),INTENT(IN) :: field  ! --> field(klon,:)
    REAL,DIMENSION(klon_mpi,size(field,2)) :: buffer_omp
    REAL :: Field3d(iim,jj_nb,size(field,2))
    INTEGER :: ip, n, nlev
    INTEGER, ALLOCATABLE, DIMENSION(:) :: index3d
    REAL,allocatable, DIMENSION(:,:) :: fieldok


    IF (size(field,1)/=klon) CALL abort_gcm('iophy::histwrite3d','Field first DIMENSION not equal to klon',1)
    nlev=size(field,2)

!   print*,'hist3d_phy mpi_rank npstn=',mpi_rank,npstn

!   DO ip=1, npstn
!    print*,'hist3d_phy mpi_rank nptabij',mpi_rank,nptabij(ip)
!   ENDDO

    CALL Gather_omp(field,buffer_omp)
!$OMP MASTER
    CALL grid1Dto2D_mpi(buffer_omp,field3d)
    if(.NOT.lpoint) THEN
     ALLOCATE(index3d(iim*jj_nb*nlev))
     ALLOCATE(fieldok(iim*jj_nb,nlev))
     IF (prt_level >= 9) write(lunout,*)'Sending ',name,' to IOIPSL'
     CALL histwrite(nid,name,itau,Field3d,iim*jj_nb*nlev,index3d)
     IF (prt_level >= 9) write(lunout,*)'Finished sending ',name,' to IOIPSL'
   else
      nlev=size(field,2)
      ALLOCATE(index3d(npstn*nlev))
      ALLOCATE(fieldok(npstn,nlev))

      if(is_sequential) then
!      klon_mpi_begin=1
!      klon_mpi_end=klon
       DO n=1, nlev
       DO ip=1, npstn
        fieldok(ip,n)=buffer_omp(nptabij(ip),n)
       ENDDO
       ENDDO
      else
       DO n=1, nlev
       DO ip=1, npstn
        IF(nptabij(ip).GE.klon_mpi_begin.AND. &
         nptabij(ip).LE.klon_mpi_end) THEN
         fieldok(ip,n)=buffer_omp(nptabij(ip)-klon_mpi_begin+1,n)
        ENDIF
       ENDDO
       ENDDO
      endif
      IF (prt_level >= 9) write(lunout,*)'Sending ',name,' to IOIPSL'
      CALL histwrite(nid,name,itau,fieldok,npstn*nlev,index3d)
      IF (prt_level >= 9) write(lunout,*)'Finished sending ',name,' to IOIPSL'
    endif 
  deallocate(index3d)
  deallocate(fieldok)
!$OMP END MASTER    

  end SUBROUTINE histwrite3d_phy_old




! ug NOUVELLE VERSION DES WRITE AVEC LA BOUCLE DO RENTREE
  SUBROUTINE histwrite2d_phy(var,field, STD_iff)
  USE dimphy
  USE mod_phys_lmdz_para
  USE ioipsl
  


#ifdef CPP_XIOS
!  USE WXIOS
#endif

  IMPLICIT NONE
  INCLUDE 'dimensions.h'
  INCLUDE 'iniprint.h'

    TYPE(ctrl_out), INTENT(IN) :: var
    REAL, DIMENSION(:), INTENT(IN) :: field
    INTEGER, INTENT(IN), OPTIONAL :: STD_iff ! ug RUSTINE POUR LES STD LEVS.....
      
    INTEGER :: iff, iff_beg, iff_end
      
    REAL,DIMENSION(klon_mpi) :: buffer_omp
    INTEGER, allocatable, DIMENSION(:) :: index2d
    REAL :: Field2d(iim,jj_nb)

    INTEGER :: ip
    REAL, ALLOCATABLE, DIMENSION(:) :: fieldok

    IF (prt_level >= 9) WRITE(lunout,*)'Begin histrwrite2d ',var%name

! ug RUSTINE POUR LES STD LEVS.....
      IF (PRESENT(STD_iff)) THEN
            iff_beg = STD_iff
            iff_end = STD_iff
      ELSE
            iff_beg = 1
            iff_end = nfiles
      END IF

  ! On regarde si on est dans la phase de définition ou d'écriture:
  IF(.NOT.vars_defined) THEN
!$OMP MASTER
      !Si phase de définition.... on définit
      DO iff=iff_beg, iff_end
         IF (clef_files(iff)) THEN
            CALL histdef2d(iff, var)
         ENDIF
      ENDDO
!$OMP END MASTER
  ELSE

    !Et sinon on.... écrit
    IF (SIZE(field)/=klon) CALL abort_gcm('iophy::histwrite2d','Field first DIMENSION not equal to klon',1)
    
    CALL Gather_omp(field,buffer_omp)    
!$OMP MASTER
    CALL grid1Dto2D_mpi(buffer_omp,Field2d)
    
! La boucle sur les fichiers:
      DO iff=iff_beg, iff_end
            IF (var%flag(iff) <= lev_files(iff) .AND. clef_files(iff)) THEN 
    
                  IF(.NOT.clef_stations(iff)) THEN
                        ALLOCATE(index2d(iim*jj_nb))
                        ALLOCATE(fieldok(iim*jj_nb))
     
                        CALL histwrite(nid_files(iff),var%name,itau_iophy,Field2d,iim*jj_nb,index2d)
#ifdef CPP_XIOS
!                        IF (iff .EQ. 1) THEN
!                              CALL wxios_write_2D(var%name, Field2d)
!                        ENDIF
#endif
                  ELSE 
                        ALLOCATE(fieldok(npstn))
                        ALLOCATE(index2d(npstn))

                        IF (is_sequential) THEN
                              DO ip=1, npstn
                                    fieldok(ip)=buffer_omp(nptabij(ip))
                              ENDDO
                             ELSE
                              DO ip=1, npstn
                                PRINT*,'histwrite2d is_sequential npstn ip namenptabij',npstn,ip,var%name,nptabij(ip)
                                     IF(nptabij(ip).GE.klon_mpi_begin.AND. &
                                        nptabij(ip).LE.klon_mpi_end) THEN
                                       fieldok(ip)=buffer_omp(nptabij(ip)-klon_mpi_begin+1)
                                     ENDIF
                              ENDDO
                       ENDIF
     
                       CALL histwrite(nid_files(iff),var%name,itau_iophy,fieldok,npstn,index2d)
                  ENDIF
                  
                deallocate(index2d)
                deallocate(fieldok)
            ENDIF !levfiles
      ENDDO
!$OMP END MASTER   
  ENDIF ! vars_defined
  IF (prt_level >= 9) WRITE(lunout,*)'End histrwrite2d ',var%name
  END SUBROUTINE histwrite2d_phy


! ug NOUVELLE VERSION DES WRITE AVEC LA BOUCLE DO RENTREE
  SUBROUTINE histwrite3d_phy(var, field, STD_iff)
  USE dimphy
  USE mod_phys_lmdz_para
  USE ioipsl
  

#ifdef CPP_XIOS
! USE WXIOS
#endif


  IMPLICIT NONE
  INCLUDE 'dimensions.h'
  INCLUDE 'iniprint.h'

    TYPE(ctrl_out), INTENT(IN) :: var
    REAL, DIMENSION(:,:), INTENT(IN) :: field ! --> field(klon,:)
    INTEGER, INTENT(IN), OPTIONAL :: STD_iff ! ug RUSTINE POUR LES STD LEVS.....
      
    INTEGER :: iff

    REAL,DIMENSION(klon_mpi,SIZE(field,2)) :: buffer_omp
    REAL :: Field3d(iim,jj_nb,SIZE(field,2))
    INTEGER :: ip, n, nlev
    INTEGER, ALLOCATABLE, DIMENSION(:) :: index3d
    REAL,ALLOCATABLE, DIMENSION(:,:) :: fieldok

  IF (prt_level >= 9) write(lunout,*)'Begin histrwrite3d ',var%name

  ! On regarde si on est dans la phase de définition ou d'écriture:
  IF(.NOT.vars_defined) THEN
      !Si phase de définition.... on définit
!$OMP MASTER
      DO iff=1, nfiles
        IF (clef_files(iff)) THEN
          CALL histdef3d(iff, var)
        ENDIF
      ENDDO
!$OMP END MASTER
  ELSE
    !Et sinon on.... écrit
    IF (SIZE(field,1)/=klon) CALL abort_gcm('iophy::histwrite3d','Field first DIMENSION not equal to klon',1)
    nlev=SIZE(field,2)


    CALL Gather_omp(field,buffer_omp)
!$OMP MASTER
    CALL grid1Dto2D_mpi(buffer_omp,field3d)


! BOUCLE SUR LES FICHIERS
      DO iff=1, nfiles
            IF (var%flag(iff) <= lev_files(iff) .AND. clef_files(iff)) THEN
                IF (.NOT.clef_stations(iff)) THEN
                        ALLOCATE(index3d(iim*jj_nb*nlev))
                        ALLOCATE(fieldok(iim*jj_nb,nlev))
                        CALL histwrite(nid_files(iff),var%name,itau_iophy,Field3d,iim*jj_nb*nlev,index3d)
#ifdef CPP_XIOS
!                        IF (iff .EQ. 1) THEN
!                              CALL wxios_write_3D(var%name, Field3d(:,:,1:klev))
!                        ENDIF
#endif
                        
                ELSE 
                        nlev=size(field,2)
                        ALLOCATE(index3d(npstn*nlev))
                        ALLOCATE(fieldok(npstn,nlev))

                        IF (is_sequential) THEN
                              DO n=1, nlev
                                    DO ip=1, npstn
                                          fieldok(ip,n)=buffer_omp(nptabij(ip),n)
                                    ENDDO
                              ENDDO
                        ELSE
                              DO n=1, nlev
                                    DO ip=1, npstn
                                                IF(nptabij(ip).GE.klon_mpi_begin.AND. &
                                                      nptabij(ip).LE.klon_mpi_end) THEN
                                                fieldok(ip,n)=buffer_omp(nptabij(ip)-klon_mpi_begin+1,n)
                                          ENDIF
                                    ENDDO
                              ENDDO
                        ENDIF
                        CALL histwrite(nid_files(iff),var%name,itau_iophy,fieldok,npstn*nlev,index3d)
                  ENDIF 
                  deallocate(index3d)
                  deallocate(fieldok)
            ENDIF
      ENDDO
!$OMP END MASTER   
  ENDIF ! vars_defined
  IF (prt_level >= 9) write(lunout,*)'End histrwrite3d ',var%name
  END SUBROUTINE histwrite3d_phy
  
end module iophy