source: LMDZ6/trunk/libf/dynphy_lonlat/calfis_loc.F90 @ 5421

!
! $Id: calfis_loc.F90 5367 2024-12-03 09:40:13Z fhourdin $
!
!
!
#ifdef CPP_PARA

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)

  !    Auteur :  P. Le Van, F. Hourdin
  !   .........
  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
  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 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 comvert_mod, ONLY: preff, presnivs
  USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS
  USE iniprint_mod_h
  USE comgeom2_mod_h
  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
  USE paramet_mod_h
  IMPLICIT NONE
  !=======================================================================
  !
  !   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 :
  !    ------------------
  INTEGER :: ngridmx
  PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )

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

IF (CPPKEY_PHYS) THEN
  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

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

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

700   CONTINUE

  firstcal = .FALSE.

ELSE
  call abort_gcm("calfis_loc", &
        "calfis_p: for now can only work with parallel physics", 1)
END IF
END SUBROUTINE calfis_loc
#endif