!
! $Id: calfis_loc.F 4368 2022-12-05 23:01:16Z jyg $
!
C
C
      SUBROUTINE calfis_loc(lafin,
     $                  jD_cur, jH_cur,
     $                  pucov,
     $                  pvcov,
     $                  pteta,
     $                  pq,
     $                  pmasse,
     $                  pps,
     $                  pp,
     $                  ppk,
     $                  pphis,
     $                  pphi,
     $                  pducov,
     $                  pdvcov,
     $                  pdteta,
     $                  pdq,
     $                  flxw,
     $                  pdufi,
     $                  pdvfi,
     $                  pdhfi,
     $                  pdqfi,
     $                  pdpsfi)
#ifdef CPP_PHYS
! If using physics
c
c    Auteur :  P. Le Van, F. Hourdin 
c   .........
      USE dimphy
      USE mod_phys_lmdz_mpi_data, mpi_root_xx=>mpi_master
      USE mod_phys_lmdz_omp_data, ONLY: klon_omp, klon_omp_begin
      USE mod_const_mpi, ONLY: COMM_LMDZ
      USE mod_interface_dyn_phys
      USE IOPHY
#endif
#ifdef CPP_PARA
      USE parallel_lmdz,ONLY:omp_chunk,using_mpi,jjb_u,jje_u,jjb_v,jje_v
     $                        ,jj_begin_dyn=>jj_begin,jj_end_dyn=>jj_end
      USE Write_Field
      Use Write_field_p
      USE Times
#endif
      USE infotrac, ONLY: nqtot, tracers
      USE control_mod, ONLY: planet_type, nsplit_phys
#ifdef CPP_PHYS
      USE callphysiq_mod, ONLY: call_physiq
#endif 
      USE comvert_mod, ONLY: preff, presnivs
      USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi

#ifdef CPP_PARA
      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"

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

      include "comgeom2.h"
      include "iniprint.h"
#ifdef CPP_MPI
      include 'mpif.h'
#endif
c    Arguments :
c    -----------
      LOGICAL,INTENT(IN) ::  lafin ! .true. for the very last call to physics
      REAL,INTENT(IN):: jD_cur, jH_cur
      REAL,INTENT(IN):: pvcov(iip1,jjb_v:jje_v,llm) ! covariant meridional velocity
      REAL,INTENT(IN):: pucov(iip1,jjb_u:jje_u,llm) ! covariant zonal velocity
      REAL,INTENT(IN):: pteta(iip1,jjb_u:jje_u,llm) ! potential temperature
      REAL,INTENT(IN):: pmasse(iip1,jjb_u:jje_u,llm) ! mass in each cell ! not used
      REAL,INTENT(IN):: pq(iip1,jjb_u:jje_u,llm,nqtot) ! tracers
      REAL,INTENT(IN):: pphis(iip1,jjb_u:jje_u) ! surface geopotential
      REAL,INTENT(IN):: pphi(iip1,jjb_u:jje_u,llm) ! geopotential

      REAL,INTENT(IN) :: pdvcov(iip1,jjb_v:jje_v,llm) ! dynamical tendency on vcov ! not used
      REAL,INTENT(IN) :: pducov(iip1,jjb_u:jje_u,llm) ! dynamical tendency on ucov
      REAL,INTENT(IN) :: pdteta(iip1,jjb_u:jje_u,llm) ! dynamical tendency on teta ! not used
      REAL,INTENT(IN) :: pdq(iip1,jjb_u:jje_u,llm,nqtot) ! dynamical tendency on tracers ! not used

      REAL,INTENT(IN) :: pps(iip1,jjb_u:jje_u) ! surface pressure (Pa)
      REAL,INTENT(IN) :: pp(iip1,jjb_u:jje_u,llmp1) ! pressure at mesh interfaces (Pa)
      REAL,INTENT(IN) :: ppk(iip1,jjb_u:jje_u,llm) ! Exner at mid-layer
      REAL,INTENT(IN) :: flxw(iip1,jjb_u:jje_u,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0)

      ! tendencies (in */s) from the physics
      REAL,INTENT(OUT) :: pdvfi(iip1,jjb_v:jje_v,llm) ! tendency on covariant meridional wind
      REAL,INTENT(OUT) :: pdufi(iip1,jjb_u:jje_u,llm) ! tendency on covariant zonal wind
      REAL,INTENT(OUT) :: pdhfi(iip1,jjb_u:jje_u,llm) ! tendency on potential temperature (K/s)
      REAL,INTENT(OUT) :: pdqfi(iip1,jjb_u:jje_u,llm,nqtot) ! tendency on tracers
      REAL,INTENT(OUT) :: pdpsfi(iip1,jjb_u:jje_u) ! tendency on surface pressure (Pa/s)

#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,itr
      REAL,ALLOCATABLE,SAVE :: zpsrf(:)
      REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
      REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
c
      REAL zrot(iip1,jjb_v:jje_v,llm) ! AdlC May 2014
      REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:), zrfi(:,:)
      REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
      REAL,ALLOCATABLE,SAVE ::  zpk(:,:)
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 :: zpk_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 :: zrfi_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,zpk_omp,zphi_omp,zphis_omp,
c$OMP+                 presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
c$OMP+                 zrfi_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
      REAL SSUM

      LOGICAL,SAVE :: firstcal=.true., debut=.true.
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
      INTEGER :: jjb,jje

c
c-----------------------------------------------------------------------
c
c    1. Initialisations :
c    --------------------
c

      klon=klon_mpi
      
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),zrfi(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))
      ALLOCATE(zpk(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$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l=1,llm
        do ig0=1,klon
          i=index_i(ig0)
          j=index_j(ig0)
          zpk(ig0,l)=ppk(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

      itr = 0
      DO iq=1,nqtot
         IF(.NOT.tracers(iq)%isAdvected) CYCLE
         itr = itr + 1
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,itr)  = pq(i,j,l,iq)
           enddo
         ENDDO
c$OMP END DO NOWAIT         
      ENDDO


c   Geopotentiel calcule par rapport a la surface locale:
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)
             zphi(ig0,l)  = pphi(i,j,l)
           enddo
         ENDDO
c$OMP END DO NOWAIT         

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

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)
             zphis(ig0)  = pphis(i,j)
           enddo
c$OMP END MASTER


c      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_dyn) kstart=2
      if (is_south_pole_dyn) 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
C  Alvaro de la Camara (May 2014)
C  46.1 Calcul de la vorticite et passage sur la grille physique
C  --------------------------------------------------------------

      jjb=jj_begin_dyn-1
      jje=jj_end_dyn+1
      if (is_north_pole_dyn) jjb=1
      if (is_south_pole_dyn) jje=jjm

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)

      DO l=1,llm
        do i=1,iim
          do j=jjb,jje
            zrot(i,j,l) = (pvcov(i+1,j,l) - pvcov(i,j,l)
     $                   + pucov(i,j+1,l) - pucov(i,j,l))
     $                   / (cu(i,j)+cu(i,j+1))
     $                   / (cv(i+1,j)+cv(i,j)) *4
          enddo
        enddo
      ENDDO


c   46.2champ 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) )
          if (j==1 .OR. j==jjp1) then !  AdlC MAY 2014
            zrfi(ig0,l) = 0 !  AdlC MAY 2014
          else
            if(i==1)then
            zrfi(ig0,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l)
     $                   +zrot(1,j-1,l)+zrot(1,j,l))   !  AdlC MAY 2014
            else
            zrfi(ig0,l)= 0.25 *(zrot(i-1,j-1,l)+zrot(i-1,j,l)
     $                   +zrot(i,j-1,l)+zrot(i,j,l))   !  AdlC MAY 2014
            endif
          endif

    
         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_dyn) 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
           zrfi(1,l)  = 0.
  
        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_dyn) 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
           zrfi(klon,l)  = 0.
        ENDDO
c$OMP END DO NOWAIT       
      endif

c On change de grille, dynamique vers physiq, pour le flux de masse verticale
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)
             flxwfi(ig0,l)  = flxw(i,j,l)
           enddo
         ENDDO
c$OMP END DO NOWAIT

c      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(zpk_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(zrfi_omp(klon,llm))  ! LG Ari 2014
        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  
	  zpk_omp(i,l)=zpk(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
          zrfi_omp(i,l)=zrfi(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

        CALL call_physiq(klon,llm,nqtot,tracers(:)%name,
     &                   debut_split,lafin_split,
     &                   jD_cur,jH_cur_split,zdt_split,
     &                   zplev_omp,zplay_omp,
     &                   zpk_omp,zphi_omp,zphis_omp,
     &                   presnivs_omp,
     &                   zufi_omp,zvfi_omp,zrfi_omp,ztfi_omp,zqfi_omp,
     &                   flxwfi_omp,pducov,
     &                   zdufi_omp,zdvfi_omp,zdtfi_omp,zdqfi_omp,
     &                   zdpsrf_omp)


         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_dyn) then
          pdhfi(:,jj_end:jj_end+1,l)=0
          pdqfi(:,jj_end:jj_end+1,l,:)=0
          pdufi(:,jj_end:jj_end+1,l)=0
          pdvfi(:,jj_end:jj_end+1,l)=0
        endif
      
       ENDDO 
c$OMP END DO NOWAIT

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

c  tendance sur la pression :
c  -----------------------------------
c      CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)

c$OMP MASTER
      kstart=1
      kend=klon

      if (is_north_pole_dyn) kstart=2
      if (is_south_pole_dyn)  kend=klon-1

!CDIR ON_ADB(index_i)
!CDIR ON_ADB(index_j) 
!cdir NODEP
        do ig0=kstart,kend
          i=index_i(ig0)
          j=index_j(ig0)
          pdpsfi(i,j) = zdpsrf(ig0)
          if (i==1) pdpsfi(iip1,j) =  zdpsrf(ig0)
         enddo          

        if (is_north_pole_dyn) then
            DO i=1,iip1
              pdpsfi(i,1)    = zdpsrf(1)
            enddo
        endif
        
        if (is_south_pole_dyn) then
            DO i=1,iip1
              pdpsfi(i,jjp1) = zdpsrf(klon)
            ENDDO
        endif
c$OMP END MASTER
cc$OMP BARRIER

c
c   62. enthalpie potentielle
c   ---------------------
      
      kstart=1
      kend=klon

      if (is_north_pole_dyn) kstart=2
      if (is_south_pole_dyn)  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_dyn) then
            DO i=1,iip1
              pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
            enddo
        endif
        
        if (is_south_pole_dyn) 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_dyn) then
!             do i=1,iip1
!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)             
!             enddo
!           endif
!           
!           if (is_south_pole_dyn) 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(:,jj_begin:jj_end,l,:)=0.
      ENDDO
c$OMP END DO NOWAIT         

C
!cdir NODEP
      itr = 0
      DO iq=1,nqtot
         IF(.NOT.tracers(iq)%isAdvected) CYCLE
         itr = itr + 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)
              pdqfi(i,j,l,iq) = zdqfi(ig0,l,itr)
              if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,itr)
            ENDDO
            
            IF (is_north_pole_dyn) then
              DO i=1,iip1
                pdqfi(i,1,l,iq)    = zdqfi(1,l,itr)
              ENDDO
            ENDIF
            
            IF (is_south_pole_dyn) then
              DO i=1,iip1
                pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,itr)
              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_dyn) then
           DO i=1,iip1
            pdufi(i,1,l)    = 0.
           ENDDO
         endif
         
         if (is_south_pole_dyn) 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_dyn) kstart=2
      if (is_south_pole_dyn)  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_dyn) 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_dyn) 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
#endif
! of #ifdef CPP_PARA
      END