Index: LMDZ6/trunk/libf/dynphy_lonlat/calfis.F
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/calfis.F	(revision 5204)
+++ 	(revision )
@@ -1,689 +1,0 @@
-!
-! $Id$
-!
-C
-C
-      SUBROUTINE calfis(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)
-c
-c    Auteur :  P. Le Van, F. Hourdin 
-c   .........
-      USE infotrac, ONLY: nqtot, tracers
-      USE control_mod, ONLY: planet_type, nsplit_phys
-#ifdef CPP_PHYS
-      USE callphysiq_mod, ONLY: call_physiq
-#endif 
-      USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
-      USE comvert_mod, ONLY: preff, presnivs
-      
-      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       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"
-
-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,jjm,llm) ! covariant meridional velocity
-      REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
-      REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
-      REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
-      REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
-      REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
-      REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
-
-      REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
-      REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
-      REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
-      ! NB: pdteta is used only to compute pcvgt which is in fact not used...
-      REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
-      ! NB: pdq is only used to compute pcvgq which is in fact not used...
-
-      REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
-      REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
-      REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
-      REAL,INTENT(IN) :: flxw(iip1,jjp1,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,jjm,llm) ! tendency on covariant meridional wind
-      REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
-      REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
-      REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
-      REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
-
-
-c    Local variables :
-c    -----------------
-
-      INTEGER i,j,l,ig0,ig,iq,itr
-      REAL zpsrf(ngridmx)
-      REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
-      REAL zphi(ngridmx,llm),zphis(ngridmx)
-c
-      REAL zrot(iip1,jjm,llm) ! AdlC May 2014
-      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
-      REAL zrfi(ngridmx,llm) ! relative wind vorticity
-      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
-      REAL zpk(ngridmx,llm)
-c
-      REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
-      REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
-c
-      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
-      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
-      REAL zdpsrf(ngridmx)
-c
-      REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
-      REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
-      REAL jH_cur_split,zdt_split
-      LOGICAL debut_split,lafin_split
-      INTEGER isplit
-
-      REAL zsin(iim),zcos(iim),z1(iim)
-      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
-      REAL unskap, pksurcp
-c
-      REAL flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
-c
-      
-      REAL SSUM
-
-      LOGICAL,SAVE :: firstcal=.true., debut=.true.
-!      REAL rdayvrai
-
-c
-c-----------------------------------------------------------------------
-c
-c    1. Initialisations :
-c    --------------------
-c
-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
-         call abort_gcm("calfis", "", 1)
-        ENDIF
-      ELSE
-        debut = .FALSE.
-      ENDIF ! of IF (firstcal)
-
-c
-c
-c-----------------------------------------------------------------------
-c   40. transformation des variables dynamiques en variables physiques:
-c   ---------------------------------------------------------------
-
-c   41. pressions au sol (en Pascals)
-c   ----------------------------------
-
-       
-      zpsrf(1) = pps(1,1)
-
-      ig0  = 2
-      DO j = 2,jjm
-         CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
-         ig0 = ig0+iim
-      ENDDO
-
-      zpsrf(ngridmx) = pps(1,jjp1)
-
-
-c   42. pression intercouches et fonction d'Exner:
-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
-      DO l = 1, llm
-        zpk(   1,l ) = ppk(1,1,l)
-        zplev( 1,l ) = pp(1,1,l)
-        ig0 = 2
-          DO j = 2, jjm
-             DO i =1, iim
-              zpk(   ig0,l ) = ppk(i,j,l)
-              zplev( ig0,l ) = pp(i,j,l)
-              ig0 = ig0 +1
-             ENDDO
-          ENDDO
-        zpk(   ngridmx,l ) = ppk(1,jjp1,l)
-        zplev( ngridmx,l ) = pp(1,jjp1,l)
-      ENDDO
-        zplev( 1,llmp1 ) = pp(1,1,llmp1)
-        ig0 = 2
-          DO j = 2, jjm
-             DO i =1, iim
-              zplev( ig0,llmp1 ) = pp(i,j,llmp1)
-              ig0 = ig0 +1
-             ENDDO
-          ENDDO
-        zplev( ngridmx,llmp1 ) = pp(1,jjp1,llmp1)
-c
-c
-
-c   43. temperature naturelle (en K) et pressions milieux couches .
-c   ---------------------------------------------------------------
-
-      DO l=1,llm
-
-         pksurcp     =  ppk(1,1,l) / cpp
-         zplay(1,l)  =  preff * pksurcp ** unskap
-         ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
-         pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
-         ig0         = 2
-
-         DO j = 2, jjm
-            DO i = 1, iim
-              pksurcp        = ppk(i,j,l) / cpp
-              zplay(ig0,l)   = preff * pksurcp ** unskap
-              ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
-              pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
-              ig0            = ig0 + 1
-            ENDDO
-         ENDDO
-
-         pksurcp       = ppk(1,jjp1,l) / cpp
-         zplay(ig0,l)  = preff * pksurcp ** unskap
-         ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
-         pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
-
-      ENDDO
-
-c   43.bis traceurs
-c   ---------------
-c
-      itr=0
-      DO iq=1,nqtot
-         IF(.NOT.tracers(iq)%isAdvected) CYCLE
-         itr = itr + 1
-         DO l=1,llm
-            zqfi(1,l,itr) = pq(1,1,l,iq)
-            ig0           = 2
-            DO j=2,jjm
-               DO i = 1, iim
-                  zqfi(ig0,l,itr) = pq(i,j,l,iq)
-                  ig0             = ig0 + 1
-               ENDDO
-            ENDDO
-            zqfi(ig0,l,itr) = pq(1,jjp1,l,iq)
-         ENDDO
-      ENDDO
-
-c   convergence dynamique pour les traceurs "EAU"
-! Earth-specific treatment of first 2 tracers (water)
-       if (planet_type=="earth") then
-        DO iq=1,2
-         DO l=1,llm
-            pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
-            ig0          = 2
-            DO j=2,jjm
-               DO i = 1, iim
-                  pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
-                  ig0             = ig0 + 1
-               ENDDO
-            ENDDO
-            pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
-         ENDDO
-        ENDDO
-       endif ! of if (planet_type=="earth")
-
-
-c   Geopotentiel calcule par rapport a la surface locale:
-c   -----------------------------------------------------
-
-      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
-      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
-      DO l=1,llm
-         DO ig=1,ngridmx
-           zphi(ig,l)=zphi(ig,l)-zphis(ig)
-         ENDDO
-      ENDDO
-
-c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
-c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux 
-c    de masse est calclue dans advtrac.F  
-c      DO l=1,llm
-c        pvervel(1,l)=pw(1,1,l) * g /apoln
-c        ig0=2
-c       DO j=2,jjm
-c           DO i = 1, iim
-c              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
-c              ig0 = ig0 + 1
-c           ENDDO
-c       ENDDO
-c        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
-c      ENDDO
-
-c
-c   45. champ u:
-c   ------------
-
-      DO 50 l=1,llm
-
-         DO 25 j=2,jjm
-            ig0 = 1+(j-2)*iim
-            zufi(ig0+1,l)= 0.5 * 
-     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
-            pcvgu(ig0+1,l)= 0.5 * 
-     $      ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
-            DO 10 i=2,iim
-               zufi(ig0+i,l)= 0.5 *
-     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
-               pcvgu(ig0+i,l)= 0.5 *
-     $         ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
-10         CONTINUE
-25      CONTINUE
-
-50    CONTINUE
-
-
-C  Alvaro de la Camara (May 2014)
-C  46.1 Calcul de la vorticite et passage sur la grille physique
-C  --------------------------------------------------------------
-      DO l=1,llm
-        do i=1,iim
-          do j=1,jjm
-            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.champ v:
-c   -----------
-
-      DO l=1,llm
-         DO j=2,jjm
-            ig0=1+(j-2)*iim
-            DO i=1,iim
-               zvfi(ig0+i,l)= 0.5 *
-     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
-               pcvgv(ig0+i,l)= 0.5 *
-     $         ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
-            ENDDO
-               zrfi(ig0 + 1,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l)
-     &                                +zrot(1,j-1,l)+zrot(1,j,l))
-            DO i=2,iim
-               zrfi(ig0 + i,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
-            ENDDO
-         ENDDO
-      ENDDO
-
-
-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 ]
-
-      DO l=1,llm
-
-         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
-         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
-         DO i=2,iim
-            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
-            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
-         ENDDO
-
-         DO i=1,iim
-            zcos(i)   = COS(rlonv(i))*z1(i)
-            zcosbis(i)= COS(rlonv(i))*z1bis(i)
-            zsin(i)   = SIN(rlonv(i))*z1(i)
-            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
-         ENDDO
-
-         zufi(1,l)  = SSUM(iim,zcos,1)/pi
-         pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
-         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
-         pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
-         zrfi(1, l) = 0.
-      ENDDO
-
-
-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 ]
-
-      DO l=1,llm
-
-         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
-         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
-         DO i=2,iim
-            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
-            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
-         ENDDO
-
-         DO i=1,iim
-            zcos(i)    = COS(rlonv(i))*z1(i)
-            zcosbis(i) = COS(rlonv(i))*z1bis(i)
-            zsin(i)    = SIN(rlonv(i))*z1(i)
-            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
-         ENDDO
-
-         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
-         pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
-         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
-         pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
-         zrfi(ngridmx, l) = 0.
-      ENDDO
-c
-c On change de grille, dynamique vers physiq, pour le flux de masse verticale
-      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
-
-c-----------------------------------------------------------------------
-c   Appel de la physique:
-c   ---------------------
-
-
-
-!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
-      zdt_split=dtphys/nsplit_phys
-      zdufic(:,:)=0.
-      zdvfic(:,:)=0.
-      zdtfic(:,:)=0.
-      zdqfic(:,:,:)=0.
-
-#ifdef CPP_PHYS
-
-       do isplit=1,nsplit_phys
-
-         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
-         debut_split=debut.and.isplit==1
-         lafin_split=lafin.and.isplit==nsplit_phys
-
-!      if (planet_type=="earth") then
-        CALL call_physiq(ngridmx,llm,nqtot,tracers(:)%name,
-     &                   debut_split,lafin_split,
-     &                   jD_cur,jH_cur_split,zdt_split,
-     &                   zplev,zplay,
-     &                   zpk,zphi,zphis,
-     &                   presnivs,
-     &                   zufi,zvfi,zrfi,ztfi,zqfi,
-     &                   flxwfi,pducov,
-     &                   zdufi,zdvfi,zdtfi,zdqfi,zdpsrf)
-!
-!      else if ( planet_type=="generic" ) then
-!
-!         CALL physiq (ngridmx,     !! ngrid
-!     .             llm,            !! nlayer
-!     .             nqtot,          !! nq
-!     .             tracers(:)%name,!! tracer names from dynamical core (given in infotrac)
-!     .             debut_split,    !! firstcall 
-!     .             lafin_split,    !! lastcall
-!     .             jD_cur,         !! pday. see leapfrog
-!     .             jH_cur_split,   !! ptime "fraction of day"
-!     .             zdt_split,      !! ptimestep
-!     .             zplev,          !! pplev
-!     .             zplay,          !! pplay
-!     .             zphi,           !! pphi
-!     .             zufi,           !! pu
-!     .             zvfi,           !! pv
-!     .             ztfi,           !! pt
-!     .             zqfi,           !! pq
-!     .             flxwfi,         !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
-!     .             zdufi,          !! pdu
-!     .             zdvfi,          !! pdv
-!     .             zdtfi,          !! pdt
-!     .             zdqfi,          !! pdq
-!     .             zdpsrf,         !! pdpsrf
-!     .             tracerdyn)      !! tracerdyn <-- utilite ???
-!
-!      endif ! of if (planet_type=="earth")
-
-         zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
-         zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
-         ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
-         zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
-
-         zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
-         zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
-         zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
-         zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
-
-       enddo ! of do isplit=1,nsplit_phys
-
-#endif
-! of #ifdef CPP_PHYS
-
-      zdufi(:,:)=zdufic(:,:)/nsplit_phys
-      zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
-      zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
-      zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
-
-
-500   CONTINUE
-
-c-----------------------------------------------------------------------
-c   transformation des tendances physiques en tendances dynamiques:
-c   ---------------------------------------------------------------
-
-c  tendance sur la pression :
-c  -----------------------------------
-
-      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
-c
-c   62. enthalpie potentielle
-c   ---------------------
-
-      DO l=1,llm
-
-         DO i=1,iip1
-          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
-          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
-         ENDDO
-
-         DO j=2,jjm
-            ig0=1+(j-2)*iim
-            DO i=1,iim
-               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
-            ENDDO
-               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
-         ENDDO
-
-      ENDDO
-
-
-c   62. humidite specifique
-c   ---------------------
-! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
-!      DO iq=1,nqtot
-!         DO l=1,llm
-!            DO i=1,iip1
-!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
-!               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
-!            ENDDO
-!            DO j=2,jjm
-!               ig0=1+(j-2)*iim
-!               DO i=1,iim
-!                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
-!               ENDDO
-!               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
-!            ENDDO
-!         ENDDO
-!      ENDDO
-
-c   63. traceurs
-c   ------------
-C     initialisation des tendances
-      pdqfi(:,:,:,:)=0.
-C
-      itr = 0
-      DO iq=1,nqtot
-         IF(.NOT.tracers(iq)%isAdvected) CYCLE
-         itr = itr + 1
-         DO l=1,llm
-            DO i=1,iip1
-               pdqfi(i,1,l,iq)    = zdqfi(1,l,itr)
-               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,itr)
-            ENDDO
-            DO j=2,jjm
-               ig0=1+(j-2)*iim
-               DO i=1,iim
-                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,itr)
-               ENDDO
-               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,itr)
-            ENDDO
-         ENDDO
-      ENDDO
-
-c   65. champ u:
-c   ------------
-
-      DO l=1,llm
-
-         DO i=1,iip1
-            pdufi(i,1,l)    = 0.
-            pdufi(i,jjp1,l) = 0.
-         ENDDO
-
-         DO j=2,jjm
-            ig0=1+(j-2)*iim
-            DO i=1,iim-1
-               pdufi(i,j,l)=
-     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
-            ENDDO
-            pdufi(iim,j,l)=
-     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
-            pdufi(iip1,j,l)=pdufi(1,j,l)
-         ENDDO
-
-      ENDDO
-
-
-c   67. champ v:
-c   ------------
-
-      DO l=1,llm
-
-         DO j=2,jjm-1
-            ig0=1+(j-2)*iim
-            DO i=1,iim
-               pdvfi(i,j,l)=
-     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
-            ENDDO
-            pdvfi(iip1,j,l) = pdvfi(1,j,l)
-         ENDDO
-      ENDDO
-
-
-c   68. champ v pres des poles:
-c   ---------------------------
-c      v = U * cos(long) + V * SIN(long)
-
-      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,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
-     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
-            pdvfi(i,1,l)=
-     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
-            pdvfi(i,jjm,l)=
-     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
-          ENDDO
-
-         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
-         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
-
-      ENDDO
-
-c-----------------------------------------------------------------------
-
-700   CONTINUE
- 
-      firstcal = .FALSE.
-
-      RETURN
-      END
Index: LMDZ6/trunk/libf/dynphy_lonlat/calfis.F90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/calfis.F90	(revision 5246)
+++ LMDZ6/trunk/libf/dynphy_lonlat/calfis.F90	(revision 5246)
@@ -0,0 +1,689 @@
+!
+! $Id$
+!
+!
+!
+SUBROUTINE calfis(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 infotrac, ONLY: nqtot, tracers
+  USE control_mod, ONLY: planet_type, nsplit_phys
+#ifdef CPP_PHYS
+  USE callphysiq_mod, ONLY: call_physiq
+#endif
+  USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
+  USE comvert_mod, ONLY: preff, presnivs
+
+  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 :
+  ! -------
+  !   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,jjm,llm) ! covariant meridional velocity
+  REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
+  REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
+  REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
+  REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
+  REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
+  REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
+
+  REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
+  REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
+  REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
+  ! ! NB: pdteta is used only to compute pcvgt which is in fact not used...
+  REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
+  ! ! NB: pdq is only used to compute pcvgq which is in fact not used...
+
+  REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
+  REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
+  REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
+  REAL,INTENT(IN) :: flxw(iip1,jjp1,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,jjm,llm) ! tendency on covariant meridional wind
+  REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
+  REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
+  REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
+  REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
+
+
+  !    Local variables :
+  !    -----------------
+
+  INTEGER :: i,j,l,ig0,ig,iq,itr
+  REAL :: zpsrf(ngridmx)
+  REAL :: zplev(ngridmx,llm+1),zplay(ngridmx,llm)
+  REAL :: zphi(ngridmx,llm),zphis(ngridmx)
+  !
+  REAL :: zrot(iip1,jjm,llm) ! AdlC May 2014
+  REAL :: zufi(ngridmx,llm), zvfi(ngridmx,llm)
+  REAL :: zrfi(ngridmx,llm) ! relative wind vorticity
+  REAL :: ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
+  REAL :: zpk(ngridmx,llm)
+  !
+  REAL :: pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
+  REAL :: pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
+  !
+  REAL :: zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
+  REAL :: zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
+  REAL :: zdpsrf(ngridmx)
+  !
+  REAL :: zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
+  REAL :: zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
+  REAL :: jH_cur_split,zdt_split
+  LOGICAL :: debut_split,lafin_split
+  INTEGER :: isplit
+
+  REAL :: zsin(iim),zcos(iim),z1(iim)
+  REAL :: zsinbis(iim),zcosbis(iim),z1bis(iim)
+  REAL :: unskap, pksurcp
+  !
+  REAL :: flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
+  !
+
+  REAL :: SSUM
+
+  LOGICAL,SAVE :: firstcal=.true., debut=.true.
+   ! REAL rdayvrai
+
+  !
+  !-----------------------------------------------------------------------
+  !
+  !    1. Initialisations :
+  !    --------------------
+  !
+  !
+  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", "", 1)
+    ENDIF
+  ELSE
+    debut = .FALSE.
+  ENDIF ! of IF (firstcal)
+
+  !
+  !
+  !-----------------------------------------------------------------------
+  !   40. transformation des variables dynamiques en variables physiques:
+  !   ---------------------------------------------------------------
+
+  !   41. pressions au sol (en Pascals)
+  !   ----------------------------------
+
+
+  zpsrf(1) = pps(1,1)
+
+  ig0  = 2
+  DO j = 2,jjm
+     CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
+     ig0 = ig0+iim
+  ENDDO
+
+  zpsrf(ngridmx) = pps(1,jjp1)
+
+
+  !   42. pression intercouches et fonction d'Exner:
+  !
+  !   -----------------------------------------------------------------
+  ! .... zplev  definis aux (llm +1) interfaces des couches  ....
+  ! .... zplay  definis aux (  llm )    milieux des couches  ....
+  !   -----------------------------------------------------------------
+
+  !    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
+  !
+   unskap   = 1./ kappa
+  !
+  DO l = 1, llm
+    zpk(   1,l ) = ppk(1,1,l)
+    zplev( 1,l ) = pp(1,1,l)
+    ig0 = 2
+      DO j = 2, jjm
+         DO i =1, iim
+          zpk(   ig0,l ) = ppk(i,j,l)
+          zplev( ig0,l ) = pp(i,j,l)
+          ig0 = ig0 +1
+         ENDDO
+      ENDDO
+    zpk(   ngridmx,l ) = ppk(1,jjp1,l)
+    zplev( ngridmx,l ) = pp(1,jjp1,l)
+  ENDDO
+    zplev( 1,llmp1 ) = pp(1,1,llmp1)
+    ig0 = 2
+      DO j = 2, jjm
+         DO i =1, iim
+          zplev( ig0,llmp1 ) = pp(i,j,llmp1)
+          ig0 = ig0 +1
+         ENDDO
+      ENDDO
+    zplev( ngridmx,llmp1 ) = pp(1,jjp1,llmp1)
+  !
+  !
+
+  !   43. temperature naturelle (en K) et pressions milieux couches .
+  !   ---------------------------------------------------------------
+
+  DO l=1,llm
+
+     pksurcp     =  ppk(1,1,l) / cpp
+     zplay(1,l)  =  preff * pksurcp ** unskap
+     ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
+     pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
+     ig0         = 2
+
+     DO j = 2, jjm
+        DO i = 1, iim
+          pksurcp        = ppk(i,j,l) / cpp
+          zplay(ig0,l)   = preff * pksurcp ** unskap
+          ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
+          pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
+          ig0            = ig0 + 1
+        ENDDO
+     ENDDO
+
+     pksurcp       = ppk(1,jjp1,l) / cpp
+     zplay(ig0,l)  = preff * pksurcp ** unskap
+     ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
+     pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
+
+  ENDDO
+
+  !   43.bis traceurs
+  !   ---------------
+  !
+  itr=0
+  DO iq=1,nqtot
+     IF(.NOT.tracers(iq)%isAdvected) CYCLE
+     itr = itr + 1
+     DO l=1,llm
+        zqfi(1,l,itr) = pq(1,1,l,iq)
+        ig0           = 2
+        DO j=2,jjm
+           DO i = 1, iim
+              zqfi(ig0,l,itr) = pq(i,j,l,iq)
+              ig0             = ig0 + 1
+           ENDDO
+        ENDDO
+        zqfi(ig0,l,itr) = pq(1,jjp1,l,iq)
+     ENDDO
+  ENDDO
+
+  !   convergence dynamique pour les traceurs "EAU"
+  ! Earth-specific treatment of first 2 tracers (water)
+   if (planet_type=="earth") then
+    DO iq=1,2
+     DO l=1,llm
+        pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
+        ig0          = 2
+        DO j=2,jjm
+           DO i = 1, iim
+              pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
+              ig0             = ig0 + 1
+           ENDDO
+        ENDDO
+        pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
+     ENDDO
+    ENDDO
+   endif ! of if (planet_type=="earth")
+
+
+  !   Geopotentiel calcule par rapport a la surface locale:
+  !   -----------------------------------------------------
+
+  CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
+  CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
+  DO l=1,llm
+     DO ig=1,ngridmx
+       zphi(ig,l)=zphi(ig,l)-zphis(ig)
+     ENDDO
+  ENDDO
+
+  !   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
+  ! JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
+  !    de masse est calclue dans advtrac.F
+   ! DO l=1,llm
+   !   pvervel(1,l)=pw(1,1,l) * g /apoln
+   !   ig0=2
+   !  DO j=2,jjm
+   !      DO i = 1, iim
+   !         pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
+   !         ig0 = ig0 + 1
+   !      ENDDO
+  !       ENDDO
+   !   pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
+   ! ENDDO
+
+  !
+  !   45. champ u:
+  !   ------------
+
+  DO l=1,llm
+
+     DO j=2,jjm
+        ig0 = 1+(j-2)*iim
+        zufi(ig0+1,l)= 0.5 * &
+              ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
+        pcvgu(ig0+1,l)= 0.5 * &
+              ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
+        DO i=2,iim
+           zufi(ig0+i,l)= 0.5 * &
+                 ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
+           pcvgu(ig0+i,l)= 0.5 * &
+                 ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
+        END DO
+     END DO
+
+  END DO
+
+
+  !  Alvaro de la Camara (May 2014)
+  !  46.1 Calcul de la vorticite et passage sur la grille physique
+  !  --------------------------------------------------------------
+  DO l=1,llm
+    do i=1,iim
+      do j=1,jjm
+        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.champ v:
+  !   -----------
+
+  DO l=1,llm
+     DO j=2,jjm
+        ig0=1+(j-2)*iim
+        DO i=1,iim
+           zvfi(ig0+i,l)= 0.5 * &
+                 ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
+           pcvgv(ig0+i,l)= 0.5 * &
+                 ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
+        ENDDO
+           zrfi(ig0 + 1,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l) &
+                 +zrot(1,j-1,l)+zrot(1,j,l))
+        DO i=2,iim
+           zrfi(ig0 + i,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
+        ENDDO
+     ENDDO
+  ENDDO
+
+
+  !   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 ]
+
+  DO l=1,llm
+
+     z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
+     z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
+     DO i=2,iim
+        z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
+        z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
+     ENDDO
+
+     DO i=1,iim
+        zcos(i)   = COS(rlonv(i))*z1(i)
+        zcosbis(i)= COS(rlonv(i))*z1bis(i)
+        zsin(i)   = SIN(rlonv(i))*z1(i)
+        zsinbis(i)= SIN(rlonv(i))*z1bis(i)
+     ENDDO
+
+     zufi(1,l)  = SSUM(iim,zcos,1)/pi
+     pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
+     zvfi(1,l)  = SSUM(iim,zsin,1)/pi
+     pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
+     zrfi(1, l) = 0.
+  ENDDO
+
+
+  !   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 ]
+
+  DO l=1,llm
+
+     z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
+     z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
+     DO i=2,iim
+        z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
+        z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
+     ENDDO
+
+     DO i=1,iim
+        zcos(i)    = COS(rlonv(i))*z1(i)
+        zcosbis(i) = COS(rlonv(i))*z1bis(i)
+        zsin(i)    = SIN(rlonv(i))*z1(i)
+        zsinbis(i) = SIN(rlonv(i))*z1bis(i)
+     ENDDO
+
+     zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
+     pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
+     zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
+     pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
+     zrfi(ngridmx, l) = 0.
+  ENDDO
+  !
+  ! On change de grille, dynamique vers physiq, pour le flux de masse verticale
+  CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
+
+  !-----------------------------------------------------------------------
+  !   Appel de la physique:
+  !   ---------------------
+
+
+
+   ! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
+  zdt_split=dtphys/nsplit_phys
+  zdufic(:,:)=0.
+  zdvfic(:,:)=0.
+  zdtfic(:,:)=0.
+  zdqfic(:,:,:)=0.
+
+#ifdef CPP_PHYS
+
+   do isplit=1,nsplit_phys
+
+     jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
+     debut_split=debut.and.isplit==1
+     lafin_split=lafin.and.isplit==nsplit_phys
+
+   ! if (planet_type=="earth") then
+    CALL call_physiq(ngridmx,llm,nqtot,tracers(:)%name, &
+          debut_split,lafin_split, &
+          jD_cur,jH_cur_split,zdt_split, &
+          zplev,zplay, &
+          zpk,zphi,zphis, &
+          presnivs, &
+          zufi,zvfi,zrfi,ztfi,zqfi, &
+          flxwfi,pducov, &
+          zdufi,zdvfi,zdtfi,zdqfi,zdpsrf)
+  !
+  !  else if ( planet_type=="generic" ) then
+  !
+  !     CALL physiq (ngridmx,     !! ngrid
+  ! .             llm,            !! nlayer
+  ! .             nqtot,          !! nq
+  ! .             tracers(:)%name,!! tracer names from dynamical core (given in infotrac)
+  ! .             debut_split,    !! firstcall
+  ! .             lafin_split,    !! lastcall
+  ! .             jD_cur,         !! pday. see leapfrog
+  ! .             jH_cur_split,   !! ptime "fraction of day"
+  ! .             zdt_split,      !! ptimestep
+  ! .             zplev,          !! pplev
+  ! .             zplay,          !! pplay
+  ! .             zphi,           !! pphi
+  ! .             zufi,           !! pu
+  ! .             zvfi,           !! pv
+  ! .             ztfi,           !! pt
+  ! .             zqfi,           !! pq
+  ! .             flxwfi,         !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
+  ! .             zdufi,          !! pdu
+  ! .             zdvfi,          !! pdv
+  ! .             zdtfi,          !! pdt
+  ! .             zdqfi,          !! pdq
+  ! .             zdpsrf,         !! pdpsrf
+  ! .             tracerdyn)      !! tracerdyn <-- utilite ???
+  !
+  !  endif ! of if (planet_type=="earth")
+
+     zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
+     zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
+     ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
+     zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
+
+     zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
+     zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
+     zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
+     zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
+
+   enddo ! of do isplit=1,nsplit_phys
+
+#endif
+  ! of #ifdef CPP_PHYS
+
+  zdufi(:,:)=zdufic(:,:)/nsplit_phys
+  zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
+  zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
+  zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
+
+
+500   CONTINUE
+
+  !-----------------------------------------------------------------------
+  !   transformation des tendances physiques en tendances dynamiques:
+  !   ---------------------------------------------------------------
+
+  !  tendance sur la pression :
+  !  -----------------------------------
+
+  CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
+  !
+  !   62. enthalpie potentielle
+  !   ---------------------
+
+  DO l=1,llm
+
+     DO i=1,iip1
+      pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
+      pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
+     ENDDO
+
+     DO j=2,jjm
+        ig0=1+(j-2)*iim
+        DO i=1,iim
+           pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
+        ENDDO
+           pdhfi(iip1,j,l) =  pdhfi(1,j,l)
+     ENDDO
+
+  ENDDO
+
+
+  !   62. humidite specifique
+  !   ---------------------
+  ! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
+   ! DO iq=1,nqtot
+   !    DO l=1,llm
+   !       DO i=1,iip1
+   !          pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
+   !          pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
+   !       ENDDO
+   !       DO j=2,jjm
+   !          ig0=1+(j-2)*iim
+   !          DO i=1,iim
+   !             pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
+   !          ENDDO
+   !          pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
+   !       ENDDO
+   !    ENDDO
+   ! ENDDO
+
+  !   63. traceurs
+  !   ------------
+  ! initialisation des tendances
+  pdqfi(:,:,:,:)=0.
+  !
+  itr = 0
+  DO iq=1,nqtot
+     IF(.NOT.tracers(iq)%isAdvected) CYCLE
+     itr = itr + 1
+     DO l=1,llm
+        DO i=1,iip1
+           pdqfi(i,1,l,iq)    = zdqfi(1,l,itr)
+           pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,itr)
+        ENDDO
+        DO j=2,jjm
+           ig0=1+(j-2)*iim
+           DO i=1,iim
+              pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,itr)
+           ENDDO
+           pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,itr)
+        ENDDO
+     ENDDO
+  ENDDO
+
+  !   65. champ u:
+  !   ------------
+
+  DO l=1,llm
+
+     DO i=1,iip1
+        pdufi(i,1,l)    = 0.
+        pdufi(i,jjp1,l) = 0.
+     ENDDO
+
+     DO j=2,jjm
+        ig0=1+(j-2)*iim
+        DO i=1,iim-1
+           pdufi(i,j,l)= &
+                 0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
+        ENDDO
+        pdufi(iim,j,l)= &
+              0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
+        pdufi(iip1,j,l)=pdufi(1,j,l)
+     ENDDO
+
+  ENDDO
+
+
+  !   67. champ v:
+  !   ------------
+
+  DO l=1,llm
+
+     DO j=2,jjm-1
+        ig0=1+(j-2)*iim
+        DO i=1,iim
+           pdvfi(i,j,l)= &
+                 0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
+        ENDDO
+        pdvfi(iip1,j,l) = pdvfi(1,j,l)
+     ENDDO
+  ENDDO
+
+
+  !   68. champ v pres des poles:
+  !   ---------------------------
+   ! v = U * cos(long) + V * SIN(long)
+
+  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,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i)) &
+              +zdvfi(ngridmx,l)*SIN(rlonv(i))
+        pdvfi(i,1,l)= &
+              0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
+        pdvfi(i,jjm,l)= &
+              0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
+      ENDDO
+
+     pdvfi(iip1,1,l)  = pdvfi(1,1,l)
+     pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
+
+  ENDDO
+
+  !-----------------------------------------------------------------------
+
+700   CONTINUE
+
+  firstcal = .FALSE.
+
+  RETURN
+END SUBROUTINE calfis
Index: LMDZ6/trunk/libf/dynphy_lonlat/calfis_loc.F
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/calfis_loc.F	(revision 5204)
+++ 	(revision )
@@ -1,1219 +1,0 @@
-!
-! $Id$
-!
-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
-      USE lmdz_mpi
-
-#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"
-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
-      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
-
-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
-          call abort_gcm("calfis_loc", "", 1)
-        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
-
-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
-
-c$OMP BARRIER
-     
-      endif
-   
-      if (MPI_rank<MPI_Size-1) then
-c$OMP BARRIER
-
-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
-
-c$OMP BARRIER
-
-
-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
-
-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 
-      call abort_gcm("calfis_loc",
-     & "calfis_p: for now can only work with parallel physics", 1)
-#endif 
-! of #ifdef CPP_PHYS
-#endif
-! of #ifdef CPP_PARA
-      END
Index: LMDZ6/trunk/libf/dynphy_lonlat/calfis_loc.F90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/calfis_loc.F90	(revision 5246)
+++ LMDZ6/trunk/libf/dynphy_lonlat/calfis_loc.F90	(revision 5246)
@@ -0,0 +1,1219 @@
+!
+! $Id$
+!
+!
+!
+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
+  !
+  !    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
+#endif
+  USE lmdz_mpi
+
+#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
+  !=======================================================================
+  !
+  !   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)
+
+#ifdef CPP_PHYS
+  ! 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
+
+  !
+  !-----------------------------------------------------------------------
+  !
+  !    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.
+
+#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
+
+!$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)
+#endif
+  ! of #ifdef CPP_PHYS
+#endif
+  ! of #ifdef CPP_PARA
+END SUBROUTINE calfis_loc
Index: LMDZ6/trunk/libf/dynphy_lonlat/gr_dyn_fi.F
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/gr_dyn_fi.F	(revision 5204)
+++ 	(revision )
@@ -1,40 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
-      IMPLICIT NONE
-c=======================================================================
-c   passage d'un champ de la grille scalaire a la grille physique
-c=======================================================================
-
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-      INTEGER im,jm,ngrid,nfield
-      REAL pdyn(im,jm,nfield)
-      REAL pfi(ngrid,nfield)
-
-      INTEGER j,ifield,ig
-
-c-----------------------------------------------------------------------
-c   calcul:
-c   -------
-
-      IF (ngrid.NE.2+(jm-2)*(im-1)) then
-         call abort_gcm("gr_dyn_fi", 'probleme de dim', 1)
-      end if
-c   traitement des poles
-      CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
-      CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
-
-c   traitement des point normaux
-      DO ifield=1,nfield
-         DO j=2,jm-1
-	    ig=2+(j-2)*(im-1)
-            CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
-         ENDDO
-      ENDDO
-
-      RETURN
-      END
Index: LMDZ6/trunk/libf/dynphy_lonlat/gr_dyn_fi.f90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/gr_dyn_fi.f90	(revision 5246)
+++ LMDZ6/trunk/libf/dynphy_lonlat/gr_dyn_fi.f90	(revision 5246)
@@ -0,0 +1,40 @@
+!
+! $Header$
+!
+SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
+  IMPLICIT NONE
+  !=======================================================================
+  !   passage d'un champ de la grille scalaire a la grille physique
+  !=======================================================================
+
+  !-----------------------------------------------------------------------
+  !   declarations:
+  !   -------------
+
+  INTEGER :: im,jm,ngrid,nfield
+  REAL :: pdyn(im,jm,nfield)
+  REAL :: pfi(ngrid,nfield)
+
+  INTEGER :: j,ifield,ig
+
+  !-----------------------------------------------------------------------
+  !   calcul:
+  !   -------
+
+  IF (ngrid.NE.2+(jm-2)*(im-1)) then
+     call abort_gcm("gr_dyn_fi", 'probleme de dim', 1)
+  end if
+  !   traitement des poles
+  CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
+  CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
+
+  !   traitement des point normaux
+  DO ifield=1,nfield
+     DO j=2,jm-1
+        ig=2+(j-2)*(im-1)
+        CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
+     ENDDO
+  ENDDO
+
+  RETURN
+END SUBROUTINE gr_dyn_fi
Index: LMDZ6/trunk/libf/dynphy_lonlat/gr_dyn_fi_p.F
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/gr_dyn_fi_p.F	(revision 5204)
+++ 	(revision )
@@ -1,44 +1,0 @@
-!
-! $Id$
-!
-      SUBROUTINE gr_dyn_fi_p(nfield,im,jm,ngrid,pdyn,pfi)
-#ifdef CPP_PARA
-! Interface with parallel physics,
-      USE mod_interface_dyn_phys
-      USE dimphy
-      USE parallel_lmdz
-      IMPLICIT NONE
-c=======================================================================
-c   passage d'un champ de la grille scalaire a la grille physique
-c=======================================================================
-
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-      INTEGER im,jm,ngrid,nfield
-      REAL pdyn(im,jm,nfield)
-      REAL pfi(ngrid,nfield)
-
-      INTEGER i,j,ig,l
-
-c-----------------------------------------------------------------------
-c   calcul:
-c   -------
-
-c      IF(ngrid.NE.2+(jm-2)*(im-1)) STOP 'probleme de dim'
-c   traitement des poles
-c   traitement des point normaux
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l=1,nfield    
-       DO ig=1,klon
-         i=index_i(ig)
-         j=index_j(ig)
-         pfi(ig,l)=pdyn(i,j,l)
-       ENDDO
-      ENDDO
-c$OMP END DO NOWAIT
-#endif
-! of #ifdef CPP_PARA
-      RETURN
-      END
Index: LMDZ6/trunk/libf/dynphy_lonlat/gr_dyn_fi_p.F90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/gr_dyn_fi_p.F90	(revision 5246)
+++ LMDZ6/trunk/libf/dynphy_lonlat/gr_dyn_fi_p.F90	(revision 5246)
@@ -0,0 +1,44 @@
+!
+! $Id$
+!
+SUBROUTINE gr_dyn_fi_p(nfield,im,jm,ngrid,pdyn,pfi)
+#ifdef CPP_PARA
+  ! Interface with parallel physics,
+  USE mod_interface_dyn_phys
+  USE dimphy
+  USE parallel_lmdz
+  IMPLICIT NONE
+  !=======================================================================
+  !   passage d'un champ de la grille scalaire a la grille physique
+  !=======================================================================
+
+  !-----------------------------------------------------------------------
+  !   declarations:
+  !   -------------
+
+  INTEGER :: im,jm,ngrid,nfield
+  REAL :: pdyn(im,jm,nfield)
+  REAL :: pfi(ngrid,nfield)
+
+  INTEGER :: i,j,ig,l
+
+  !-----------------------------------------------------------------------
+  !   calcul:
+  !   -------
+
+   ! IF(ngrid.NE.2+(jm-2)*(im-1)) STOP 'probleme de dim'
+  !   traitement des poles
+  !   traitement des point normaux
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+  DO l=1,nfield
+   DO ig=1,klon
+     i=index_i(ig)
+     j=index_j(ig)
+     pfi(ig,l)=pdyn(i,j,l)
+   ENDDO
+  ENDDO
+!$OMP END DO NOWAIT
+#endif
+  ! of #ifdef CPP_PARA
+  RETURN
+END SUBROUTINE gr_dyn_fi_p
Index: LMDZ6/trunk/libf/dynphy_lonlat/gr_fi_dyn.F
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/gr_fi_dyn.F	(revision 5204)
+++ 	(revision )
@@ -1,40 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
-      IMPLICIT NONE
-c=======================================================================
-c   passage d'un champ de la grille scalaire a la grille physique
-c=======================================================================
-
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-      INTEGER im,jm,ngrid,nfield
-      REAL pdyn(im,jm,nfield)
-      REAL pfi(ngrid,nfield)
-
-      INTEGER i,j,ifield,ig
-
-c-----------------------------------------------------------------------
-c   calcul:
-c   -------
-
-      DO ifield=1,nfield
-c   traitement des poles
-         DO i=1,im
-            pdyn(i,1,ifield)=pfi(1,ifield)
-            pdyn(i,jm,ifield)=pfi(ngrid,ifield)
-         ENDDO
-
-c   traitement des point normaux
-         DO j=2,jm-1
-	    ig=2+(j-2)*(im-1)
-            CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
-	    pdyn(im,j,ifield)=pdyn(1,j,ifield)
-         ENDDO
-      ENDDO
-
-      RETURN
-      END
Index: LMDZ6/trunk/libf/dynphy_lonlat/gr_fi_dyn.f90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/gr_fi_dyn.f90	(revision 5246)
+++ LMDZ6/trunk/libf/dynphy_lonlat/gr_fi_dyn.f90	(revision 5246)
@@ -0,0 +1,40 @@
+!
+! $Header$
+!
+SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
+  IMPLICIT NONE
+  !=======================================================================
+  !   passage d'un champ de la grille scalaire a la grille physique
+  !=======================================================================
+
+  !-----------------------------------------------------------------------
+  !   declarations:
+  !   -------------
+
+  INTEGER :: im,jm,ngrid,nfield
+  REAL :: pdyn(im,jm,nfield)
+  REAL :: pfi(ngrid,nfield)
+
+  INTEGER :: i,j,ifield,ig
+
+  !-----------------------------------------------------------------------
+  !   calcul:
+  !   -------
+
+  DO ifield=1,nfield
+  !   traitement des poles
+     DO i=1,im
+        pdyn(i,1,ifield)=pfi(1,ifield)
+        pdyn(i,jm,ifield)=pfi(ngrid,ifield)
+     ENDDO
+
+  !   traitement des point normaux
+     DO j=2,jm-1
+        ig=2+(j-2)*(im-1)
+        CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
+        pdyn(im,j,ifield)=pdyn(1,j,ifield)
+     ENDDO
+  ENDDO
+
+  RETURN
+END SUBROUTINE gr_fi_dyn
Index: LMDZ6/trunk/libf/dynphy_lonlat/gr_fi_dyn_p.F
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/gr_fi_dyn_p.F	(revision 5204)
+++ 	(revision )
@@ -1,56 +1,0 @@
-!
-! $Id$
-!
-      SUBROUTINE gr_fi_dyn_p(nfield,ngrid,im,jm,pfi,pdyn)
-#ifdef CPP_PARA
-! Interface with parallel physics,
-      USE mod_interface_dyn_phys
-      USE dimphy
-      USE parallel_lmdz
-      IMPLICIT NONE
-c=======================================================================
-c   passage d'un champ de la grille scalaire a la grille physique
-c=======================================================================
-
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-      INTEGER im,jm,ngrid,nfield
-      REAL pdyn(im,jm,nfield)
-      REAL pfi(ngrid,nfield)
-
-      INTEGER i,j,ifield,ig
-
-c-----------------------------------------------------------------------
-c   calcul:
-c   -------
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO ifield=1,nfield
-
-        do ig=1,klon
-          i=index_i(ig)
-          j=index_j(ig)
-          pdyn(i,j,ifield)=pfi(ig,ifield)
-          if (i==1) pdyn(im,j,ifield)=pdyn(i,j,ifield)
-	enddo
-
-c   traitement des poles
-      if (pole_nord) then
-        do i=1,im
-	  pdyn(i,1,ifield)=pdyn(1,1,ifield)
-	enddo
-      endif
-       
-      if (pole_sud) then
-        do i=1,im
-	  pdyn(i,jm,ifield)=pdyn(1,jm,ifield)
-	enddo
-      endif
-      
-      ENDDO
-c$OMP END DO NOWAIT
-#endif
-! of #ifdef CPP_PARA
-      RETURN
-      END
Index: LMDZ6/trunk/libf/dynphy_lonlat/gr_fi_dyn_p.F90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/gr_fi_dyn_p.F90	(revision 5246)
+++ LMDZ6/trunk/libf/dynphy_lonlat/gr_fi_dyn_p.F90	(revision 5246)
@@ -0,0 +1,56 @@
+!
+! $Id$
+!
+SUBROUTINE gr_fi_dyn_p(nfield,ngrid,im,jm,pfi,pdyn)
+#ifdef CPP_PARA
+  ! Interface with parallel physics,
+  USE mod_interface_dyn_phys
+  USE dimphy
+  USE parallel_lmdz
+  IMPLICIT NONE
+  !=======================================================================
+  !   passage d'un champ de la grille scalaire a la grille physique
+  !=======================================================================
+
+  !-----------------------------------------------------------------------
+  !   declarations:
+  !   -------------
+
+  INTEGER :: im,jm,ngrid,nfield
+  REAL :: pdyn(im,jm,nfield)
+  REAL :: pfi(ngrid,nfield)
+
+  INTEGER :: i,j,ifield,ig
+
+  !-----------------------------------------------------------------------
+  !   calcul:
+  !   -------
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+  DO ifield=1,nfield
+
+    do ig=1,klon
+      i=index_i(ig)
+      j=index_j(ig)
+      pdyn(i,j,ifield)=pfi(ig,ifield)
+      if (i==1) pdyn(im,j,ifield)=pdyn(i,j,ifield)
+    enddo
+
+  !   traitement des poles
+  if (pole_nord) then
+    do i=1,im
+      pdyn(i,1,ifield)=pdyn(1,1,ifield)
+    enddo
+  endif
+
+  if (pole_sud) then
+    do i=1,im
+      pdyn(i,jm,ifield)=pdyn(1,jm,ifield)
+    enddo
+  endif
+
+  ENDDO
+!$OMP END DO NOWAIT
+#endif
+  ! of #ifdef CPP_PARA
+  RETURN
+END SUBROUTINE gr_fi_dyn_p
