source: LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/lmdz_calfis_loc.F90 @ 5117

#ifdef CPP_PARA
MODULE lmdz_calfis_loc
  IMPLICIT NONE
  PRIVATE
  PUBLIC calfis_loc
CONTAINS

  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)

    USE dimphy
    USE lmdz_phys_mpi_data, mpi_root_xx=>mpi_master
    USE lmdz_phys_omp_data, ONLY: klon_omp, klon_omp_begin
    USE mod_const_mpi, ONLY: COMM_LMDZ
    USE mod_interface_dyn_phys
    USE IOPHY

    USE lmdz_mpi

    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 lmdz_write_field
    Use Write_field_p
    USE Times
    USE infotrac, ONLY: nqtot, tracers
    USE control_mod, ONLY: planet_type, nsplit_phys
    USE callphysiq_mod, ONLY: call_physiq
    USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS

    USE comvert_mod, ONLY: preff, presnivs
    USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi

    !=======================================================================

    !   1. rearrangement des tableaux et transformation
    !  variables dynamiques  >  variables physiques
    !   2. calcul des termes physiques
    !   3. retransformation des tendances physiques en tendances dynamiques

    !   remarques:
    !   ----------

    !    - les vents sont donnes dans la physique par leurs composantes
    !  naturelles.
    !    - la variable thermodynamique de la physique est une variable
    !  intensive :   T
    !  pour la dynamique on prend    T * ( preff / p(l) ) **kappa
    !    - les deux seules variables dependant de la geometrie necessaires
    !  pour la physique sont la latitude pour le rayonnement et
    !  l'aire de la maille quand on veut integrer une grandeur
    !  horizontalement.
    !    - les points de la physique sont les points scalaires de la
    !  la dynamique; numerotation:
    !      1 pour le pole nord
    !      (jjm-1)*iim pour l'interieur du domaine
    !      ngridmx pour le pole sud
    !  ---> ngridmx=2+(jjm-1)*iim

    ! Input :
    ! -------
    !   ecritphy        frequence d'ecriture (en jours)de histphy
    !   pucov           covariant zonal velocity
    !   pvcov           covariant meridional velocity
    !   pteta           potential temperature
    !   pps             surface pressure
    !   pmasse          masse d'air dans chaque maille
    !   pts             surface temperature  (K)
    !   callrad         clef d'appel au rayonnement

    !    Output :
    !    --------
    !    pdufi          tendency for the natural zonal velocity (ms-1)
    !    pdvfi          tendency for the natural meridional velocity
    !    pdhfi          tendency for the potential temperature
    !    pdtsfi         tendency for the surface temperature

    !    pdtrad         radiative tendencies  \  both input
    !    pfluxrad       radiative fluxes      /  and output

    !=======================================================================

    !-----------------------------------------------------------------------

    !    0.  Declarations :
    !    ------------------

    include "dimensions.h"
    include "paramet.h"

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

    include "comgeom2.h"
    include "iniprint.h"
    !    Arguments :
    !    -----------
    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)


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

  INTEGER :: i,j,l,ig0,ig,iq,itr
  REAL,ALLOCATABLE,SAVE :: zpsrf(:)
  REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
  REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)

  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(:,:)

  REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
  REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)

  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

  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
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

  LOGICAL,SAVE :: first_omp=.TRUE.
!$OMP THREADPRIVATE(first_omp)

  REAL :: zsin(iim),zcos(iim),z1(iim)
  REAL :: zsinbis(iim),zcosbis(iim),z1bis(iim)
  REAL :: unskap, pksurcp

  REAL :: SSUM

  LOGICAL,SAVE :: firstcal=.TRUE., debut=.TRUE.
!$OMP THREADPRIVATE(firstcal,debut)

  REAL,SAVE,DIMENSION(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
  INTEGER :: ierr
  INTEGER,DIMENSION(MPI_STATUS_SIZE,4) :: Status
  INTEGER, DIMENSION(4) :: Req
  REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
  INTEGER :: k,kstart,kend
  INTEGER :: offset
  INTEGER :: jjb,jje

    IF (CPPKEY_PHYS) THEN

  !-----------------------------------------------------------------------

  !    1. Initialisations :
  !    --------------------

  klon=klon_mpi

  IF ( firstcal )  THEN
    debut = .TRUE.
    IF (ngridmx/=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
      CALL abort_gcm("calfis_loc", "", 1)
    ENDIF
!$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))
!$OMP END MASTER
!$OMP BARRIER
  ELSE
      debut = .FALSE.
  ENDIF


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

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

!$OMP MASTER
  CALL start_timer(timer_physic)
!$OMP END MASTER

!$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
!$OMP END MASTER


  !   42. pression intercouches :

  !   -----------------------------------------------------------------
  ! .... zplev  definis aux (llm +1) interfaces des couches  ....
  ! .... zplay  definis aux (  llm )    milieux des couches  ....
  !   -----------------------------------------------------------------

  !    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....

   unskap   = 1./ kappa

  !  print *,omp_rank,'klon--->',klon
!$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
!$OMP END DO NOWAIT

!$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
!$OMP END DO NOWAIT


  !   43. temperature naturelle (en K) et pressions milieux couches .
  !   ---------------------------------------------------------------
!$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
!$OMP END DO NOWAIT

  !   43.bis traceurs
  !   ---------------

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


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

!$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
!$OMP END DO NOWAIT

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

!$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
!$OMP END MASTER


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

!$OMP BARRIER

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

  !   45. champ u:
  !   ------------

  kstart=1
  kend=klon

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

!$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
!$OMP END DO NOWAIT

  !  Alvaro de la Camara (May 2014)
  !  46.1 Calcul de la vorticite et passage sur la grille physique
  !  --------------------------------------------------------------

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

!$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


  !   46.2champ v:
  !   -----------

!$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
!$OMP END DO NOWAIT

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

  IF (is_north_pole_dyn) THEN
!$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
!$OMP END DO NOWAIT
  ENDIF


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

  IF (is_south_pole_dyn) THEN
!$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
!$OMP END DO NOWAIT
  ENDIF

  ! On change de grille, dynamique vers physiq, pour le flux de masse verticale
!$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
!$OMP END DO NOWAIT

   ! CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)

  !-----------------------------------------------------------------------
  !   Appel de la physique:
  !   ---------------------


!$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

!$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.

  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


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

!$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
!$OMP BARRIER

!$OMP MASTER
  CALL stop_timer(timer_physic)
!$OMP END MASTER

  IF (using_mpi) THEN

  IF (MPI_rank>0) THEN
!$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
!$OMP END DO NOWAIT

!$OMP BARRIER

!$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)
!$OMP END MASTER

!$OMP BARRIER

  ENDIF

  IF (MPI_rank<MPI_Size-1) THEN
!$OMP BARRIER

!$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)
!$OMP END MASTER

  ENDIF

!$OMP BARRIER


!$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)
!$OMP END MASTER

!$OMP BARRIER

  ENDIF ! using_mpi


!$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
!$OMP END DO NOWAIT

!$OMP MASTER
    pdpsfi(:,jj_begin)=0

   IF (.NOT. is_south_pole_dyn) THEN
     pdpsfi(:,jj_end:jj_end+1)=0
   endif
!$OMP END MASTER
  !-----------------------------------------------------------------------
  !   transformation des tendances physiques en tendances dynamiques:
  !   ---------------------------------------------------------------

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

!$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
!$OMP END MASTER
  !c$OMP BARRIER

  !   62. enthalpie potentielle
  !   ---------------------

  kstart=1
  kend=klon

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

!$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
!$OMP END DO NOWAIT

  !   62. humidite specifique
  !   ---------------------
  ! 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

  !   63. traceurs
  !   ------------
  ! initialisation des tendances

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

  !cdir NODEP
  itr = 0
  DO iq=1,nqtot
     IF(.NOT.tracers(iq)%isAdvected) CYCLE
     itr = itr + 1
!$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
!$OMP END DO NOWAIT
  ENDDO

  !   65. champ u:
  !   ------------
!$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
!$OMP END DO NOWAIT

  !   67. champ v:
  !   ------------

  kstart=1
  kend=klon

  IF (is_north_pole_dyn) kstart=2
  IF (is_south_pole_dyn)  kend=klon-1-iim

!$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
!$OMP END DO NOWAIT


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

  IF (is_north_pole_dyn) THEN
!$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
!$OMP END DO NOWAIT

  ENDIF

  IF (is_south_pole_dyn) THEN
!$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
!$OMP END DO NOWAIT

  ENDIF
  !-----------------------------------------------------------------------

  firstcal = .FALSE.
END IF
  END SUBROUTINE calfis_loc

END MODULE lmdz_calfis_loc
#endif