module mod_phys_mpi
  
  logical, save :: monocpu
!  logical, save :: parallel
   
  INTEGER,SAVE :: iiphy_begin
  INTEGER,SAVE :: iiphy_end
  INTEGER,SAVE :: jjphy_begin
  INTEGER,SAVE :: jjphy_end
  INTEGER,SAVE :: jjphy_nb
  INTEGER,SAVE :: ijphy_begin
  INTEGER,SAVE :: ijphy_end
  INTEGER,SAVE :: ijphy_nb
  INTEGER,SAVE :: klon_begin
  INTEGER,SAVE :: klon_end
  INTEGER,SAVE :: klon_mpi
  
  INTEGER,SAVE,allocatable,dimension(:) :: jjphy_para_nb
  INTEGER,SAVE,allocatable,dimension(:) :: jjphy_para_begin
  INTEGER,SAVE,allocatable,dimension(:) :: jjphy_para_end
  INTEGER,SAVE,dimension(:),allocatable :: Liste_i
  INTEGER,SAVE,dimension(:),allocatable :: Liste_j
  INTEGER,SAVE,dimension(:),allocatable :: klon_para_nb
  INTEGER,SAVE,dimension(:),allocatable :: klon_para_begin
  INTEGER,SAVE,dimension(:),allocatable :: klon_para_end 

  
  INTEGER,save :: phy_rank
  INTEGER,save :: phy_size
  

  INTEGER,save :: klev, klevp1,klevm1,klon,kfdia,kidia
  INTEGER,save :: klon2
  INTEGER,save :: nbtr 
  INTEGER,save :: kdlon, kflev
!$OMP THREADPRIVATE(klon,kfdia,kidia,kdlon)

 
contains
  
   subroutine InitDimphy

#ifdef CPP_PARA
   use parallel
#endif

   implicit none
#include "dimensions90.h"
#include "paramet90.h" 
    integer :: rank
    integer :: i,j,Pos,Index
    logical,save :: First=.true.
#ifndef CPP_PARA
   logical, save :: pole_nord = .true.
   logical, save :: pole_sud = .true.
#endif

#ifdef CPP_PARA      
      monocpu=.false.
#else
      monocpu = .true.
#endif
      
      if (monocpu) then
        
	phy_rank=0
	phy_size=1
	
#ifdef CPP_PARA      
      else
 
        phy_rank=mpi_rank
	phy_size=mpi_size
#endif
      
      endif

      if (First) then
        allocate(jjphy_para_nb(0:phy_size-1))
        allocate(jjphy_para_begin(0:phy_size-1))
        allocate(jjphy_para_end(0:phy_size-1))
        allocate(klon_para_nb(0:phy_size-1))
        allocate(klon_para_begin(0:phy_size-1))
        allocate(klon_para_end(0:phy_size-1))
      endif
	
      klev=llm
      klevp1=klev+1
      klevm1=klev-1
      klon2=iim*(jjm-1)+2
      
     
!        do i=0,phy_size-1
!          klon_para_nb(i)=(klon2/phy_size)
!          if (i<MOD(klon2,phy_size)) klon_para_nb(i)=klon_para_nb(i)+1
!          if (i==0) then 
!            klon_para_begin(i)=1
!          else 
!            klon_para_begin(i)=klon_para_end(i-1)+1
!          endif
!          klon_para_end(i)=klon_para_begin(i)+klon_para_nb(i)-1
!        enddo
      
      if (First) then
        do i=0,phy_size-1
          klon_para_nb(i)=(klon2/phy_size)
          if (i<MOD(klon2,phy_size)) klon_para_nb(i)=klon_para_nb(i)+1
        enddo
      endif
        
      do i=0,phy_size-1
        if (i==0) then 
          klon_para_begin(i)=1
        else 
          klon_para_begin(i)=klon_para_end(i-1)+1
        endif
        klon_para_end(i)=klon_para_begin(i)+klon_para_nb(i)-1
      enddo
	
      klon=klon_para_nb(phy_rank)
      klon_mpi=klon_para_nb(phy_rank)
      klon_begin=klon_para_begin(phy_rank)
      klon_end=klon_para_end(phy_rank)
	
      nbtr=nqmx-2+1/(nqmx-1)

      kflev=klev
    
      if (.not.first) then
        deallocate(liste_i,liste_j)
      endif
      
      
      allocate(Liste_i(klon))
      allocate(Liste_j(klon))
    
      Index=1
      if (Pole_Nord) then
        Liste_i(1)=1
        Liste_j(1)=1
        iiphy_begin=1
        jjphy_begin=1
	ijphy_begin=1
        Index=Index+1
      endif
    
      Pos=2
      do j=2,jjm
        do i=1,iim
      
	  if (Pos==klon_begin) then
            iiphy_begin=i
            jjphy_begin=j
	    ijphy_begin=j*iip1+i
          endif
       
          if (Pos==klon_end) then
            iiphy_end=i
            jjphy_end=j
	    ijphy_end=j*iip1+i
          endif
        
          if (Pos>=klon_begin .AND. Pos<=klon_end) then
            Liste_i(Index)=i
            Liste_j(Index)=j
            Index=Index+1
          endif
          Pos=Pos+1
        enddo
      enddo
    
      if (pole_sud) then
        Liste_i(Index)=1
        Liste_j(Index)=jjp1
        iiphy_end=1
        jjphy_end=jjp1
	ijphy_end=jjp1*iip1+1
      endif
    
      ijphy_nb=ijphy_end-ijphy_begin+1
      jjphy_nb=jjphy_end-jjphy_begin+1
	
	    
      Pos=2
      rank=0
      jjphy_para_begin(rank)=1
    
      do j=2,jjm
        do i=1,iim
          if (Pos==klon_para_begin(rank)) then
            jjphy_para_begin(rank)=j
          endif
          if (Pos==klon_para_end(rank)) then
            jjphy_para_end(rank)=j
            rank=rank+1
          endif
          Pos=Pos+1
        enddo
      enddo
      jjphy_para_end(rank)=jjp1
      
      First=.false.    
   end subroutine InitDimphy
   

   subroutine GatherField(Fields,Fieldr,ll)
#ifdef CPP_PARA
      USE parallel, ONLY : COMM_LMDZ
#endif
     implicit none
#ifdef CPP_PARA      
     include 'mpif.h'
#endif
     INTEGER :: ll
     REAL, dimension(klon_mpi,ll) :: Fields
     REAL, dimension(klon2,ll) :: Fieldr     
     REAL, dimension(klon2*ll) :: Field_tmp     
     INTEGER, dimension(0:phy_size-1) :: displs
     INTEGER, dimension(0:phy_size-1) :: sendcounts
     INTEGER :: l,Pos,rank,ierr
     INTEGER :: klon_b,klon_e, Nb
     
     Pos=1     
     do rank=0,phy_size-1
       klon_b=klon_para_begin(rank)
       klon_e=klon_para_end(rank)
       Nb=klon_para_nb(rank)
       displs(rank)=Pos-1
       sendcounts(rank)=Nb*ll
       Pos=Pos+Nb*ll
     enddo
         
      if (monocpu) then
        Fieldr(:,:)=Fields(:,:)
      else
#ifdef CPP_PARA
        call MPI_Gatherv(Fields,klon_mpi*ll,MPI_REAL8,Field_tmp,sendcounts,  &
	                 displs,MPI_REAL8,0,COMM_LMDZ,ierr)
        Pos=1     
        do rank=0,phy_size-1
          klon_b=klon_para_begin(rank)
          klon_e=klon_para_end(rank)
          Nb=klon_para_nb(rank)
          do l=1,ll
           Fieldr(klon_b:klon_e,l)=Field_tmp(Pos:Pos+Nb-1)
	   Pos=Pos+Nb
          enddo
        enddo
#endif
     endif
      
   end subroutine GatherField

   subroutine AllGatherField(Fields,Fieldr,ll)
#ifdef CPP_PARA
      USE parallel, ONLY : COMM_LMDZ
#endif
     implicit none
#ifdef CPP_PARA      
     include 'mpif.h'
#endif   
     INTEGER :: ll
     REAL, dimension(klon_mpi,ll) :: Fields
     REAL, dimension(klon2,ll) :: Fieldr     
     
     INTEGER :: ierr
      
      IF (monocpu) THEN
        Fieldr(:,:)=Fields(:,:)
      ELSE 
       call GatherField(Fields,Fieldr,ll)
#ifdef CPP_PARA             
       call MPI_BCAST(Fieldr,klon2*ll,MPI_REAL8,0,COMM_LMDZ,ierr)
#endif 
      ENDIF  
   end subroutine AllGatherField
   
           
   subroutine ScatterField(Fields,Fieldr,ll)
#ifdef CPP_PARA
      USE parallel, ONLY : COMM_LMDZ
#endif
     implicit none
#ifdef CPP_PARA      
     include 'mpif.h'
#endif
     INTEGER :: ll
     REAL, dimension(klon2,ll) :: Fields
     REAL, dimension(klon_mpi,ll) :: Fieldr     
     REAL, dimension(klon2*ll) :: Field_tmp     
     INTEGER, dimension(0:phy_size-1) :: displs
     INTEGER, dimension(0:phy_size-1) :: sendcounts
     INTEGER :: l,Pos,rank,ierr
     INTEGER :: klon_b,klon_e, Nb

     Pos=1    
     do rank=0,phy_size-1
       klon_b=klon_para_begin(rank)
       klon_e=klon_para_end(rank)
       Nb=klon_para_nb(rank)
       displs(rank)=Pos-1
       sendcounts(rank)=Nb*ll
       do l=1,ll
         Field_tmp(Pos:Pos+Nb-1)=Fields(klon_b:klon_e,l)
	 Pos=Pos+Nb
       enddo
     enddo
           
     if (monocpu) then
       Fieldr(:,:)=Fields(:,:)
     else
#ifdef CPP_PARA      
       call MPI_Scatterv(Field_tmp,sendcounts,displs,MPI_REAL8,         &
	                 Fieldr,klon_mpi*ll,MPI_REAL8,0,COMM_LMDZ,ierr)
#endif
     endif
      
        
   end subroutine ScatterField
   
   
end module mod_phys_mpi