!
! $Id: calfis_p.F 1823 2013-07-31 10:38:37Z jghattas $
!
C
C
      SUBROUTINE calfis_p(lafin,
     $                  jD_cur, jH_cur,
     $                  pucov,
     $                  pvcov,
     $                  pteta,
     $                  pq,
     $                  pmasse,
     $                  pps,
     $                  pp,
     $                  ppk,
     $                  pphis,
     $                  pphi,
     $                  pducov,
     $                  pdvcov,
     $                  pdteta,
     $                  pdq,
     $                  flxw,
     $                  clesphy0,
     $                  pdufi,
     $                  pdvfi,
     $                  pdhfi,
     $                  pdqfi,
     $                  pdpsfi)
#ifdef CPP_PHYS
! If using physics
      USE dimphy
      USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root 
      USE mod_interface_dyn_phys
      USE IOPHY
#endif
      USE parallel_lmdz, ONLY : omp_chunk, using_mpi
      USE Write_Field
      Use Write_field_p
      USE Times
      USE infotrac
      USE control_mod

      IMPLICIT NONE
c=======================================================================
c
c   1. rearrangement des tableaux et transformation
c      variables dynamiques  >  variables physiques
c   2. calcul des termes physiques
c   3. retransformation des tendances physiques en tendances dynamiques
c
c   remarques:
c   ----------
c
c    - les vents sont donnes dans la physique par leurs composantes 
c      naturelles.
c    - la variable thermodynamique de la physique est une variable
c      intensive :   T 
c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
c    - les deux seules variables dependant de la geometrie necessaires
c      pour la physique sont la latitude pour le rayonnement et 
c      l'aire de la maille quand on veut integrer une grandeur 
c      horizontalement.
c    - les points de la physique sont les points scalaires de la 
c      la dynamique; numerotation:
c          1 pour le pole nord
c          (jjm-1)*iim pour l'interieur du domaine
c          ngridmx pour le pole sud
c      ---> ngridmx=2+(jjm-1)*iim
c
c     Input :
c     -------
c       ecritphy        frequence d'ecriture (en jours)de histphy
c       pucov           covariant zonal velocity
c       pvcov           covariant meridional velocity 
c       pteta           potential temperature
c       pps             surface pressure
c       pmasse          masse d'air dans chaque maille
c       pts             surface temperature  (K)
c       callrad         clef d'appel au rayonnement
c
c    Output :
c    --------
c        pdufi          tendency for the natural zonal velocity (ms-1)
c        pdvfi          tendency for the natural meridional velocity 
c        pdhfi          tendency for the potential temperature
c        pdtsfi         tendency for the surface temperature
c
c        pdtrad         radiative tendencies  \  both input
c        pfluxrad       radiative fluxes      /  and output
c
c=======================================================================
c
c-----------------------------------------------------------------------
c
c    0.  Declarations :
c    ------------------

#include "dimensions.h"
#include "paramet.h"
#include "temps.h"

      INTEGER ngridmx
      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )

#include "comconst.h"
#include "comvert.h"
#include "comgeom2.h"
#include "iniprint.h"
#ifdef CPP_MPI
      include 'mpif.h'
#endif
c    Arguments :
c    -----------
      LOGICAL  lafin
!      REAL heure
      REAL, intent(in):: jD_cur, jH_cur
      REAL pvcov(iip1,jjm,llm)
      REAL pucov(iip1,jjp1,llm)
      REAL pteta(iip1,jjp1,llm)
      REAL pmasse(iip1,jjp1,llm)
      REAL pq(iip1,jjp1,llm,nqtot)
      REAL pphis(iip1,jjp1)
      REAL pphi(iip1,jjp1,llm)
c
      REAL pdvcov(iip1,jjm,llm)
      REAL pducov(iip1,jjp1,llm)
      REAL pdteta(iip1,jjp1,llm)
      REAL pdq(iip1,jjp1,llm,nqtot)
      REAL flxw(iip1,jjp1,llm)  ! Flux de masse verticale sur la grille dynamique
c
      REAL pps(iip1,jjp1)
      REAL pp(iip1,jjp1,llmp1)
      REAL ppk(iip1,jjp1,llm)
c
      REAL pdvfi(iip1,jjm,llm)
      REAL pdufi(iip1,jjp1,llm)
      REAL pdhfi(iip1,jjp1,llm)
      REAL pdqfi(iip1,jjp1,llm,nqtot)
      REAL pdpsfi(iip1,jjp1)

      INTEGER        longcles
      PARAMETER    ( longcles = 20 )
      REAL clesphy0( longcles )

#ifdef CPP_PHYS
! Ehouarn: for now calfis_p needs some informations from physics to compile
c    Local variables :
c    -----------------

      INTEGER i,j,l,ig0,ig,iq,iiq
      REAL,ALLOCATABLE,SAVE :: zpsrf(:)
      REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
      REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
c
      REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
      REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
c
      REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
      REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
c
      REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
      REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
      REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
      REAL,SAVE,ALLOCATABLE ::  flxwfi(:,:)     ! Flux de masse verticale sur la grille physiq

c
      REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
      REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
      REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:) 
      REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
      REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:)
      REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:)
      REAL,SAVE,ALLOCATABLE ::  flxwfi_omp(:,:)     ! Flux de masse verticale sur la grille physiq

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Introduction du splitting (FH)
! Question pour Yann :
! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent
! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il
! soit allocatable (plutot par exemple que de passer une dimension
! dépendant du process en argument des routines) et que, du coup,
! le SAVE évite d'avoir à refaire l'allocation à chaque appel.
! Tu confirmes ?
! J'ai suivi le même principe pour les zdufic_omp
! Mais c'est surement bien que tu controles.
! 

      REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:)
      REAL jH_cur_split,zdt_split
      LOGICAL debut_split,lafin_split
      INTEGER isplit
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
c$OMP+                 presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
c$OMP+                 zqfi_omp,zdufi_omp,zdvfi_omp,
c$OMP+                 zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp,
c$OMP+                 zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp)       

      LOGICAL,SAVE :: first_omp=.true.
c$OMP THREADPRIVATE(first_omp)
      
      REAL zsin(iim),zcos(iim),z1(iim)
      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
      REAL unskap, pksurcp
c
cIM diagnostique PVteta, Amip2
      INTEGER ntetaSTD
      PARAMETER(ntetaSTD=3)
      REAL rtetaSTD(ntetaSTD)
      DATA rtetaSTD/350., 380., 405./ ! Earth-specific values, beware !!
      REAL PVteta(klon,ntetaSTD)
      
      
      REAL SSUM

      LOGICAL firstcal, debut
      DATA firstcal/.true./
      SAVE firstcal,debut
c$OMP THREADPRIVATE(firstcal,debut)
      
      REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
      INTEGER :: ierr
#ifdef CPP_MPI
      INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
#else
      INTEGER,dimension(1,4) :: Status
#endif
      INTEGER, dimension(4) :: Req
      REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
      integer :: k,kstart,kend
      INTEGER :: offset  

      LOGICAL tracerdyn
c
c-----------------------------------------------------------------------
c
c    1. Initialisations :
c    --------------------
c

      klon=klon_mpi
      
      PVteta(:,:)=0.
            
c
      IF ( firstcal )  THEN
        debut = .TRUE.
        IF (ngridmx.NE.2+(jjm-1)*iim) THEN
         write(lunout,*) 'STOP dans calfis'
         write(lunout,*) 
     &   'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
         write(lunout,*) '  ngridmx  jjm   iim   '
         write(lunout,*) ngridmx,jjm,iim
         STOP
        ENDIF
c$OMP MASTER
      ALLOCATE(zpsrf(klon))
      ALLOCATE(zplev(klon,llm+1),zplay(klon,llm))
      ALLOCATE(zphi(klon,llm),zphis(klon))
      ALLOCATE(zufi(klon,llm), zvfi(klon,llm))
      ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
      ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
      ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
      ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm))
      ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
      ALLOCATE(zdpsrf(klon))
      ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
      ALLOCATE(flxwfi(klon,llm))
c$OMP END MASTER
c$OMP BARRIER	  
      ELSE
          debut = .FALSE.
      ENDIF

c
c
c-----------------------------------------------------------------------
c   40. transformation des variables dynamiques en variables physiques:
c   ---------------------------------------------------------------

c   41. pressions au sol (en Pascals)
c   ----------------------------------

c$OMP MASTER
      call start_timer(timer_physic)
c$OMP END MASTER

c$OMP MASTER             
!CDIR ON_ADB(index_i)
!CDIR ON_ADB(index_j) 
      do ig0=1,klon
        i=index_i(ig0)
        j=index_j(ig0)
        zpsrf(ig0)=pps(i,j)
      enddo
c$OMP END MASTER


c   42. pression intercouches :
c
c   -----------------------------------------------------------------
c     .... zplev  definis aux (llm +1) interfaces des couches  ....
c     .... zplay  definis aux (  llm )    milieux des couches  .... 
c   -----------------------------------------------------------------

c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
c
       unskap   = 1./ kappa
c
c      print *,omp_rank,'klon--->',klon
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 1, llmp1
!CDIR ON_ADB(index_i)
!CDIR ON_ADB(index_j) 
        do ig0=1,klon
          i=index_i(ig0)
          j=index_j(ig0)
          zplev( ig0,l ) = pp(i,j,l)
        enddo
      ENDDO
c$OMP END DO NOWAIT
c
c

c   43. temperature naturelle (en K) et pressions milieux couches .
c   ---------------------------------------------------------------
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l=1,llm
!CDIR ON_ADB(index_i)
!CDIR ON_ADB(index_j) 
        do ig0=1,klon
          i=index_i(ig0)
          j=index_j(ig0)
          pksurcp        = ppk(i,j,l) / cpp
          zplay(ig0,l)   = preff * pksurcp ** unskap
          ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
        enddo

      ENDDO
c$OMP END DO NOWAIT

c   43.bis traceurs
c   ---------------
c

      DO iq=1,nqtot
         iiq=niadv(iq)
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
         DO l=1,llm
!CDIR ON_ADB(index_i)
!CDIR ON_ADB(index_j) 
           do ig0=1,klon
             i=index_i(ig0)
             j=index_j(ig0)
             zqfi(ig0,l,iq)  = pq(i,j,l,iiq)
           enddo
         ENDDO
c$OMP END DO NOWAIT	 
      ENDDO


c   Geopotentiel calcule par rapport a la surface locale:
c   -----------------------------------------------------

      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)

      CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)

c$OMP BARRIER

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l=1,llm
	 DO ig=1,klon
	   zphi(ig,l)=zphi(ig,l)-zphis(ig)
	 ENDDO
      ENDDO
c$OMP END DO NOWAIT
      

c
c   45. champ u:
c   ------------

      kstart=1
      kend=klon
      
      if (is_north_pole) kstart=2
      if (is_south_pole) kend=klon-1
      
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
      DO l=1,llm
!CDIR ON_ADB(index_i)
!CDIR ON_ADB(index_j) 
!CDIR SPARSE
        do ig0=kstart,kend
          i=index_i(ig0)
          j=index_j(ig0)
          if (i==1) then
            zufi(ig0,l)= 0.5 *(  pucov(iim,j,l)/cu(iim,j)
     $                         + pucov(1,j,l)/cu(1,j) )
          else
            zufi(ig0,l)= 0.5*(  pucov(i-1,j,l)/cu(i-1,j) 
     $                       + pucov(i,j,l)/cu(i,j) )
          endif
        enddo
      ENDDO
c$OMP END DO NOWAIT

c   46.champ v:
c   -----------

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l=1,llm
!CDIR ON_ADB(index_i)
!CDIR ON_ADB(index_j) 
        DO ig0=kstart,kend
          i=index_i(ig0)
          j=index_j(ig0)
          zvfi(ig0,l)= 0.5 *(  pvcov(i,j-1,l)/cv(i,j-1) 
     $                       + pvcov(i,j,l)/cv(i,j) )
    
         ENDDO
      ENDDO
c$OMP END DO NOWAIT

c   47. champs de vents aux pole nord   
c   ------------------------------
c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]

      if (is_north_pole) then
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
        DO l=1,llm

           z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
           DO i=2,iim
              z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
           ENDDO
  
           DO i=1,iim
              zcos(i)   = COS(rlonv(i))*z1(i)
              zsin(i)   = SIN(rlonv(i))*z1(i)
           ENDDO
  
           zufi(1,l)  = SSUM(iim,zcos,1)/pi
           zvfi(1,l)  = SSUM(iim,zsin,1)/pi
  
        ENDDO
c$OMP END DO NOWAIT      
      endif


c   48. champs de vents aux pole sud:
c   ---------------------------------
c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]

      if (is_south_pole) then
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
        DO l=1,llm
  
         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
           DO i=2,iim
             z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
	   ENDDO
  
           DO i=1,iim
              zcos(i)    = COS(rlonv(i))*z1(i)
              zsin(i)    = SIN(rlonv(i))*z1(i)
	   ENDDO
  
           zufi(klon,l)  = SSUM(iim,zcos,1)/pi
           zvfi(klon,l)  = SSUM(iim,zsin,1)/pi
        ENDDO
c$OMP END DO NOWAIT       
      endif


      IF (is_sequential.and.(planet_type=="earth")) THEN
#ifdef CPP_PHYS
! PVtheta calls tetalevel, which is in the physics
cIM calcul PV a teta=350, 380, 405K
        CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
     $           ztfi,zplay,zplev,
     $           ntetaSTD,rtetaSTD,PVteta)
c
#endif
      ENDIF

c On change de grille, dynamique vers physiq, pour le flux de masse verticale
      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)

c-----------------------------------------------------------------------
c   Appel de la physique:
c   ---------------------


c$OMP BARRIER
      if (first_omp) then
        klon=klon_omp

        allocate(zplev_omp(klon,llm+1))
        allocate(zplay_omp(klon,llm))
        allocate(zphi_omp(klon,llm))
        allocate(zphis_omp(klon))
        allocate(presnivs_omp(llm))
        allocate(zufi_omp(klon,llm))
        allocate(zvfi_omp(klon,llm))
        allocate(ztfi_omp(klon,llm))
        allocate(zqfi_omp(klon,llm,nqtot))
        allocate(zdufi_omp(klon,llm))
        allocate(zdvfi_omp(klon,llm))
        allocate(zdtfi_omp(klon,llm))
        allocate(zdqfi_omp(klon,llm,nqtot))
        allocate(zdufic_omp(klon,llm))
        allocate(zdvfic_omp(klon,llm))
        allocate(zdtfic_omp(klon,llm))
        allocate(zdqfic_omp(klon,llm,nqtot))
        allocate(zdpsrf_omp(klon))
        allocate(flxwfi_omp(klon,llm))
	first_omp=.false.
      endif
       
	   
      klon=klon_omp
      offset=klon_omp_begin-1
      
      do l=1,llm+1
        do i=1,klon
          zplev_omp(i,l)=zplev(offset+i,l)
	enddo 
      enddo
	  
       do l=1,llm
        do i=1,klon  
	  zplay_omp(i,l)=zplay(offset+i,l)
	enddo 
      enddo
	
      do l=1,llm
        do i=1,klon
	  zphi_omp(i,l)=zphi(offset+i,l)
	enddo 
      enddo
	
      do i=1,klon
	zphis_omp(i)=zphis(offset+i)
      enddo 
     
	
      do l=1,llm
        presnivs_omp(l)=presnivs(l)
      enddo 
	
      do l=1,llm
        do i=1,klon
	  zufi_omp(i,l)=zufi(offset+i,l)
	enddo 
      enddo
	
      do l=1,llm
        do i=1,klon
	  zvfi_omp(i,l)=zvfi(offset+i,l)
	enddo 
      enddo
	
      do l=1,llm
        do i=1,klon
	  ztfi_omp(i,l)=ztfi(offset+i,l)
	enddo 
      enddo
	
      do iq=1,nqtot
        do l=1,llm
          do i=1,klon
            zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
	  enddo
	enddo 
      enddo
	
      do l=1,llm
        do i=1,klon
	  zdufi_omp(i,l)=zdufi(offset+i,l)
	enddo 
      enddo
	
      do l=1,llm
        do i=1,klon
	  zdvfi_omp(i,l)=zdvfi(offset+i,l)
	enddo 
      enddo
	
      do l=1,llm
        do i=1,klon
          zdtfi_omp(i,l)=zdtfi(offset+i,l)
	enddo 
      enddo
	
      do iq=1,nqtot
        do l=1,llm
          do i=1,klon
	    zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
	  enddo 
        enddo
      enddo
      	
      do i=1,klon
	zdpsrf_omp(i)=zdpsrf(offset+i)
      enddo 

      do l=1,llm
        do i=1,klon
          flxwfi_omp(i,l)=flxwfi(offset+i,l)
	enddo 
      enddo
      
c$OMP BARRIER
      
!$OMP MASTER
!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
!$OMP END MASTER
      zdt_split=dtphys/nsplit_phys
      zdufic_omp(:,:)=0.
      zdvfic_omp(:,:)=0.
      zdtfic_omp(:,:)=0.
      zdqfic_omp(:,:,:)=0.

#ifdef CPP_PHYS
      do isplit=1,nsplit_phys

         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
         debut_split=debut.and.isplit==1
         lafin_split=lafin.and.isplit==nsplit_phys

      if (planet_type=="earth") then

      CALL physiq (klon,
     .             llm,
     .             debut_split,
     .             lafin_split,
     .             jD_cur,
     .             jH_cur_split,
     .             zdt_split,
     .             zplev_omp,
     .             zplay_omp,
     .             zphi_omp,
     .             zphis_omp,
     .             presnivs_omp,
     .             clesphy0,
     .             zufi_omp,
     .             zvfi_omp,
     .             ztfi_omp,
     .             zqfi_omp,
c#ifdef INCA
     .             flxwfi_omp,
c#endif
     .             zdufi_omp,
     .             zdvfi_omp,
     .             zdtfi_omp,
     .             zdqfi_omp,
     .             zdpsrf_omp,
cIM diagnostique PVteta, Amip2          
     .             pducov,
     .             PVteta)

      else if ( planet_type=="generic" ) then

      CALL physiq (klon,     !! ngrid
     .             llm,            !! nlayer
     .             nqtot,          !! nq
     .             tname,          !! tracer names from dynamical core (given in infotrac)
     .             debut_split,    !! firstcall 
     .             lafin_split,    !! lastcall
     .             jD_cur,         !! pday. see leapfrog_p 
     .             jH_cur_split,   !! ptime "fraction of day"
     .             zdt_split,      !! ptimestep
     .             zplev_omp,  !! pplev
     .             zplay_omp,  !! pplay
     .             zphi_omp,   !! pphi
     .             zufi_omp,   !! pu
     .             zvfi_omp,   !! pv
     .             ztfi_omp,   !! pt
     .             zqfi_omp,   !! pq
     .             flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
     .             zdufi_omp,  !! pdu
     .             zdvfi_omp,  !! pdv
     .             zdtfi_omp,  !! pdt
     .             zdqfi_omp,  !! pdq
     .             zdpsrf_omp, !! pdpsrf
     .             tracerdyn)      !! tracerdyn <-- utilite ???

      endif ! of if (planet_type=="earth")


         zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
         zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
         ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
         zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split

         zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
         zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
         zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
         zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)

      enddo

#endif
! of #ifdef CPP_PHYS

      zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
      zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
      zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
      zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
c$OMP BARRIER

      do l=1,llm+1
        do i=1,klon
          zplev(offset+i,l)=zplev_omp(i,l)
	enddo 
      enddo
	  
       do l=1,llm
        do i=1,klon  
	  zplay(offset+i,l)=zplay_omp(i,l)
	enddo 
      enddo
	
      do l=1,llm
        do i=1,klon
	  zphi(offset+i,l)=zphi_omp(i,l)
	enddo 
      enddo
	

      do i=1,klon
	zphis(offset+i)=zphis_omp(i)
      enddo 
     
	
      do l=1,llm
        presnivs(l)=presnivs_omp(l)
      enddo 
	
      do l=1,llm
        do i=1,klon
	  zufi(offset+i,l)=zufi_omp(i,l)
	enddo 
      enddo
	
      do l=1,llm
        do i=1,klon
	  zvfi(offset+i,l)=zvfi_omp(i,l)
	enddo 
      enddo
	
      do l=1,llm
        do i=1,klon
	  ztfi(offset+i,l)=ztfi_omp(i,l)
	enddo 
      enddo
	
      do iq=1,nqtot
        do l=1,llm
          do i=1,klon
            zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
	  enddo
	enddo 
      enddo
	
      do l=1,llm
        do i=1,klon
	  zdufi(offset+i,l)=zdufi_omp(i,l)
	enddo 
      enddo
	
      do l=1,llm
        do i=1,klon
	  zdvfi(offset+i,l)=zdvfi_omp(i,l)
	enddo 
      enddo
	
      do l=1,llm
        do i=1,klon
          zdtfi(offset+i,l)=zdtfi_omp(i,l)
	enddo 
      enddo
	
      do iq=1,nqtot
        do l=1,llm
          do i=1,klon
	    zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
	  enddo 
        enddo
      enddo
      	
      do i=1,klon
	zdpsrf(offset+i)=zdpsrf_omp(i)
      enddo 
      

      klon=klon_mpi
500   CONTINUE
c$OMP BARRIER

c$OMP MASTER
      call stop_timer(timer_physic)
c$OMP END MASTER

      IF (using_mpi) THEN
            
      if (MPI_rank>0) then

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
       DO l=1,llm      
        du_send(1:iim,l)=zdufi(1:iim,l)
        dv_send(1:iim,l)=zdvfi(1:iim,l)
       ENDDO
c$OMP END DO NOWAIT       

c$OMP BARRIER
#ifdef CPP_MPI 
c$OMP MASTER
!$OMP CRITICAL (MPI)
        call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
     &                   COMM_LMDZ,Req(1),ierr)
        call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
     &                  COMM_LMDZ,Req(2),ierr)
!$OMP END CRITICAL (MPI)
c$OMP END MASTER
#endif
c$OMP BARRIER
     
      endif
   
      if (MPI_rank<MPI_Size-1) then
c$OMP BARRIER
#ifdef CPP_MPI 
c$OMP MASTER      
!$OMP CRITICAL (MPI)
        call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
     &                 COMM_LMDZ,Req(3),ierr)
        call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
     &                 COMM_LMDZ,Req(4),ierr)
!$OMP END CRITICAL (MPI)
c$OMP END MASTER
#endif
      endif

c$OMP BARRIER


#ifdef CPP_MPI 
c$OMP MASTER    
!$OMP CRITICAL (MPI)
      if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
        call MPI_WAITALL(4,Req(1),Status,ierr)
      else if (MPI_rank>0) then
        call MPI_WAITALL(2,Req(1),Status,ierr)
      else if (MPI_rank <MPI_Size-1) then
        call MPI_WAITALL(2,Req(3),Status,ierr)
      endif
!$OMP END CRITICAL (MPI)
c$OMP END MASTER
#endif

c$OMP BARRIER     

      ENDIF ! using_mpi
      
      
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l=1,llm
            
        zdufi2(1:klon,l)=zdufi(1:klon,l)
        zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
            
        zdvfi2(1:klon,l)=zdvfi(1:klon,l)
        zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l) 
  
         pdhfi(:,jj_begin,l)=0
         pdqfi(:,jj_begin,l,:)=0
         pdufi(:,jj_begin,l)=0
         pdvfi(:,jj_begin,l)=0
         
         if (.not. is_south_pole) then
           pdhfi(:,jj_end,l)=0
           pdqfi(:,jj_end,l,:)=0
           pdufi(:,jj_end,l)=0
           pdvfi(:,jj_end,l)=0
         endif
      
       ENDDO 
c$OMP END DO NOWAIT

c$OMP MASTER
       pdpsfi(:,jj_begin)=0    
       if (.not. is_south_pole) then
	 pdpsfi(:,jj_end)=0
       endif
c$OMP END MASTER
c-----------------------------------------------------------------------
c   transformation des tendances physiques en tendances dynamiques:
c   ---------------------------------------------------------------

c  tendance sur la pression :
c  -----------------------------------
      CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
c
c   62. enthalpie potentielle
c   ---------------------
      
      kstart=1
      kend=klon

      if (is_north_pole) kstart=2
      if (is_south_pole)  kend=klon-1

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l=1,llm

!CDIR ON_ADB(index_i)
!CDIR ON_ADB(index_j) 
!cdir NODEP
        do ig0=kstart,kend
          i=index_i(ig0)
          j=index_j(ig0)
          pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
          if (i==1) pdhfi(iip1,j,l) =  cpp * zdtfi(ig0,l) / ppk(i,j,l)
         enddo          

        if (is_north_pole) then
            DO i=1,iip1
              pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
            enddo
        endif
        
        if (is_south_pole) then
            DO i=1,iip1
              pdhfi(i,jjp1,l) = cpp *  zdtfi(klon,l)/ ppk(i,jjp1,l)
            ENDDO
        endif
      ENDDO
c$OMP END DO NOWAIT
      
c   62. humidite specifique
c   ---------------------
! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
!      DO iq=1,nqtot
!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
!         DO l=1,llm
!!!cdir NODEP 
!           do ig0=kstart,kend
!             i=index_i(ig0)
!             j=index_j(ig0)
!             pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq) 
!             if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq) 
!           enddo
!           
!           if (is_north_pole) then
!             do i=1,iip1
!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)             
!             enddo
!           endif
!           
!           if (is_south_pole) then
!             do i=1,iip1
!               pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq) 
!             enddo
!           endif
!         ENDDO
!c$OMP END DO NOWAIT
!      ENDDO

c   63. traceurs
c   ------------
C     initialisation des tendances

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l=1,llm
        pdqfi(:,:,l,:)=0.
      ENDDO
c$OMP END DO NOWAIT	 

C
!cdir NODEP
      DO iq=1,nqtot
         iiq=niadv(iq)
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
         DO l=1,llm
!CDIR ON_ADB(index_i)
!CDIR ON_ADB(index_j) 
!cdir NODEP           
	     DO ig0=kstart,kend
              i=index_i(ig0)
              j=index_j(ig0)
              pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
              if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
            ENDDO
	    
	    IF (is_north_pole) then
	      DO i=1,iip1
                pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
	      ENDDO
	    ENDIF
	    
	    IF (is_south_pole) then
	      DO i=1,iip1
                pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
	      ENDDO
	    ENDIF
	    
         ENDDO
c$OMP END DO NOWAIT	 
      ENDDO
      
c   65. champ u:
c   ------------
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l=1,llm
!CDIR ON_ADB(index_i)
!CDIR ON_ADB(index_j) 
!cdir NODEP
         do ig0=kstart,kend
           i=index_i(ig0)
           j=index_j(ig0)
           
           if (i/=iim) then
             pdufi(i,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
           endif
           
           if (i==1) then
              pdufi(iim,j,l)=0.5*(  zdufi2(ig0,l)
     $                            + zdufi2(ig0+iim-1,l))*cu(iim,j)
             pdufi(iip1,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
           endif
         
         enddo
         
         if (is_north_pole) then
           DO i=1,iip1
            pdufi(i,1,l)    = 0.
           ENDDO
         endif
         
         if (is_south_pole) then
           DO i=1,iip1
            pdufi(i,jjp1,l) = 0.
           ENDDO
         endif
         
      ENDDO
c$OMP END DO NOWAIT

c   67. champ v:
c   ------------

      kstart=1
      kend=klon

      if (is_north_pole) kstart=2
      if (is_south_pole)  kend=klon-1-iim
      
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
      DO l=1,llm
!CDIR ON_ADB(index_i)
!CDIR ON_ADB(index_j) 
!cdir NODEP
        do ig0=kstart,kend
           i=index_i(ig0)
           j=index_j(ig0)
           pdvfi(i,j,l)=0.5*(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))*cv(i,j)
           if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
     $	                                    zdvfi2(ig0+iim,l))
     $				          *cv(i,j)
        enddo
         
      ENDDO
c$OMP END DO NOWAIT


c   68. champ v pres des poles:
c   ---------------------------
c      v = U * cos(long) + V * SIN(long)

      if (is_north_pole) then

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
        DO l=1,llm

          DO i=1,iim
            pdvfi(i,1,l)=
     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
       
            pdvfi(i,1,l)=
     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
          ENDDO

          pdvfi(iip1,1,l)  = pdvfi(1,1,l)

        ENDDO
c$OMP END DO NOWAIT

      endif    
      
      if (is_south_pole) then

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
         DO l=1,llm
  
           DO i=1,iim
              pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
     $        +zdvfi(klon,l)*SIN(rlonv(i))

              pdvfi(i,jjm,l)=
     $        0.5*(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))*cv(i,jjm)
           ENDDO

           pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)

        ENDDO
c$OMP END DO NOWAIT
     
      endif
c-----------------------------------------------------------------------

700   CONTINUE
 
      firstcal = .FALSE.

#else 
      write(lunout,*) 
     & "calfis_p: for now can only work with parallel physics" 
      stop 
#endif 
! of #ifdef CPP_PHYS
      RETURN
      END