source: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/calfis_p.F @ 1227

!
! $Id: calfis_p.F 1222 2009-08-07 11:48:33Z jghattas $
!
C
C
      SUBROUTINE calfis_p(lafin,
     $                  jD_cur, jH_cur,
     $                  pucov,
     $                  pvcov,
     $                  pteta,
     $                  pq,
     $                  pmasse,
     $                  pps,
     $                  pp,
     $                  ppk,
     $                  pphis,
     $                  pphi,
     $                  pducov,
     $                  pdvcov,
     $                  pdteta,
     $                  pdq,
     $                  flxw,
     $                  clesphy0,
     $                  pdufi,
     $                  pdvfi,
     $                  pdhfi,
     $                  pdqfi,
     $                  pdpsfi)
#ifdef CPP_EARTH
! Ehouarn: For now, calfis_p needs Earth physics
c
c    Auteur :  P. Le Van, F. Hourdin 
c   .........
      USE dimphy
      USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root 
      USE parallel, ONLY : omp_chunk, using_mpi
      USE mod_interface_dyn_phys
      USE Write_Field
      Use Write_field_p
      USE Times
      USE IOPHY
      USE infotrac

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

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

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

#include "comconst.h"
#include "comvert.h"
#include "comgeom2.h"
#include "control.h"
#ifdef CPP_MPI
      include 'mpif.h'
#endif
c    Arguments :
c    -----------
      LOGICAL  lafin
      REAL heure

      REAL pvcov(iip1,jjm,llm)
      REAL pucov(iip1,jjp1,llm)
      REAL pteta(iip1,jjp1,llm)
      REAL pmasse(iip1,jjp1,llm)
      REAL pq(iip1,jjp1,llm,nqtot)
      REAL pphis(iip1,jjp1)
      REAL pphi(iip1,jjp1,llm)
c
      REAL pdvcov(iip1,jjm,llm)
      REAL pducov(iip1,jjp1,llm)
      REAL pdteta(iip1,jjp1,llm)
      REAL pdq(iip1,jjp1,llm,nqtot)
c
      REAL pps(iip1,jjp1)
      REAL pp(iip1,jjp1,llmp1)
      REAL ppk(iip1,jjp1,llm)
c
      REAL pdvfi(iip1,jjm,llm)
      REAL pdufi(iip1,jjp1,llm)
      REAL pdhfi(iip1,jjp1,llm)
      REAL pdqfi(iip1,jjp1,llm,nqtot)
      REAL pdpsfi(iip1,jjp1)

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


c    Local variables :
c    -----------------

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

c
      REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
      REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
      REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:) 
      REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
      REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
c      REAL,ALLOCATABLE,SAVE :: pvervel_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

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

      LOGICAL,SAVE :: first_omp=.true.
c$OMP THREADPRIVATE(first_omp)
      
      REAL zsin(iim),zcos(iim),z1(iim)
      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
      REAL unskap, pksurcp
c
cIM diagnostique PVteta, Amip2
      INTEGER ntetaSTD
      PARAMETER(ntetaSTD=3)
      REAL rtetaSTD(ntetaSTD)
      DATA rtetaSTD/350., 380., 405./
      REAL PVteta(klon,ntetaSTD)
      
      REAL flxw(iip1,jjp1,llm)  ! Flux de masse verticale sur la grille dynamique
      
      REAL SSUM

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

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

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

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

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

c$OMP MASTER             
      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
        do ig0=1,klon
          i=index_i(ig0)
          j=index_j(ig0)
          zplev( ig0,l ) = pp(i,j,l)
        enddo
      ENDDO
c$OMP END DO NOWAIT
c
c

c   43. temperature naturelle (en K) et pressions milieux couches .
c   ---------------------------------------------------------------
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l=1,llm

        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
c          pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
        enddo

      ENDDO
c$OMP END DO NOWAIT

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

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

c   convergence dynamique pour les traceurs "EAU"
! Earth-specific treatment of first 2 tracers (water)
      if (planet_type=="earth") then
       DO iq=1,2
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
         DO l=1,llm
           do ig0=1,klon
             i=index_i(ig0)
             j=index_j(ig0)
c             pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
           enddo
         ENDDO
c$OMP END DO NOWAIT      
       ENDDO
      endif ! of if (planet_type=="earth")


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

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

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

c$OMP BARRIER

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l=1,llm
         DO ig=1,klon
           zphi(ig,l)=zphi(ig,l)-zphis(ig)
         ENDDO
      ENDDO
c$OMP END DO NOWAIT
      
c   ....  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_p.F  
c
cc$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
c      DO l=1,llm
c        do ig0=1,klon
c           i=index_i(ig0)
c           j=index_j(ig0)
c           pvervel(ig0,l) = pw(i,j,l)*g* unsaire(i,j)
c        enddo
c       if (is_north_pole) pvervel(1,l)=pw(1,1,l)*g /apoln
c       if (is_south_pole) pvervel(klon,l)=pw(1,jjp1,l)*g/apols
c      ENDDO
cc$OMP END DO NOWAIT

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

      kstart=1
      kend=klon
      
      if (is_north_pole) kstart=2
      if (is_south_pole) kend=klon-1
      
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
      DO l=1,llm
        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) )
c            pcvgu(ig0,l)= 0.5*(  pducov(iim,j,l)/cu(iim,j) 
c     $                         + pducov(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) )
c            pcvgu(ig0,l)= 0.5*(  pducov(i-1,j,l)/cu(i-1,j) 
c     $                        + pducov(i,j,l)/cu(i,j) )
          endif
        enddo
      ENDDO
c$OMP END DO NOWAIT
c   46.champ v:
c   -----------
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l=1,llm
        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) )
    
c          pcvgv(ig0+i,l)= 0.5 * (  pdvcov(i,j-1,l)/cv(i,j-1) 
c     $                           + pdvcov(i,j,l)/cv(i,j) )
         ENDDO
      ENDDO
c$OMP END DO NOWAIT

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

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

           z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
c           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)
c              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)
c              zcosbis(i)= COS(rlonv(i))*z1bis(i)
              zsin(i)   = SIN(rlonv(i))*z1(i)
c              zsinbis(i)= SIN(rlonv(i))*z1bis(i)
           ENDDO
  
           zufi(1,l)  = SSUM(iim,zcos,1)/pi
c           pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
           zvfi(1,l)  = SSUM(iim,zsin,1)/pi
c           pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
  
        ENDDO
c$OMP END DO NOWAIT      
      endif


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

      if (is_south_pole) then
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
        DO l=1,llm
  
         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
c         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)
c           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)
c              zcosbis(i) = COS(rlonv(i))*z1bis(i)
              zsin(i)    = SIN(rlonv(i))*z1(i)
c              zsinbis(i) = SIN(rlonv(i))*z1bis(i)
           ENDDO
  
           zufi(klon,l)  = SSUM(iim,zcos,1)/pi
c           pcvgu(klon,l) = SSUM(iim,zcosbis,1)/pi
           zvfi(klon,l)  = SSUM(iim,zsin,1)/pi
c           pcvgv(klon,l) = SSUM(iim,zsinbis,1)/pi

        ENDDO
c$OMP END DO NOWAIT       
      endif


      IF (is_sequential) THEN
        if (planet_type=="earth") then
#ifdef CPP_EARTH
c
cIM calcul PV a teta=350, 380, 405K
        CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
     $           ztfi,zplay,zplev,
     $           ntetaSTD,rtetaSTD,PVteta)
c
#endif
        endif 
      ENDIF

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

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

cc$OMP  PARALLEL DEFAULT(NONE) 
cc$OMP+ PRIVATE(i,l,offset,iq)
cc$OMP+ SHARED(klon_omp_nb,nqtot,klon_omp_begin,
cc$OMP+        debut,lafin,rdayvrai,heure,dtphys,zplev,zplay,
cc$OMP+        zphi,zphis,presnivs,clesphy0,zufi,zvfi,ztfi,
cc$OMP+        zqfi,pvervel,zdufi,zdvfi,zdtfi,zdqfi,zdpsrf)

c PRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
c c$OMP+                 presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
c c$OMP+                 zqfi_omp,pvervel_omp,zdufi_omp,zdvfi_omp,
c c$OMP+                 zdtfi_omp,zdqfi_omp,zdpsrf_omp)

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

        allocate(zplev_omp(klon,llm+1))
        allocate(zplay_omp(klon,llm))
        allocate(zphi_omp(klon,llm))
        allocate(zphis_omp(klon))
        allocate(presnivs_omp(llm))
        allocate(zufi_omp(klon,llm))
        allocate(zvfi_omp(klon,llm))
        allocate(ztfi_omp(klon,llm))
        allocate(zqfi_omp(klon,llm,nqtot))
c        allocate(pvervel_omp(klon,llm))
        allocate(zdufi_omp(klon,llm))
        allocate(zdvfi_omp(klon,llm))
        allocate(zdtfi_omp(klon,llm))
        allocate(zdqfi_omp(klon,llm,nqtot))
        allocate(zdpsrf_omp(klon))
        allocate(flxwfi_omp(klon,llm))
        first_omp=.false.
      endif
       
           
      klon=klon_omp
      offset=klon_omp_begin-1
      
      do l=1,llm+1
        do i=1,klon
          zplev_omp(i,l)=zplev(offset+i,l)
        enddo 
      enddo
          
       do l=1,llm
        do i=1,klon  
          zplay_omp(i,l)=zplay(offset+i,l)
        enddo 
      enddo
        
      do l=1,llm
        do i=1,klon
          zphi_omp(i,l)=zphi(offset+i,l)
        enddo 
      enddo
        
      do i=1,klon
        zphis_omp(i)=zphis(offset+i)
      enddo 
     
        
      do l=1,llm
        presnivs_omp(l)=presnivs(l)
      enddo 
        
      do l=1,llm
        do i=1,klon
          zufi_omp(i,l)=zufi(offset+i,l)
        enddo 
      enddo
        
      do l=1,llm
        do i=1,klon
          zvfi_omp(i,l)=zvfi(offset+i,l)
        enddo 
      enddo
        
      do l=1,llm
        do i=1,klon
          ztfi_omp(i,l)=ztfi(offset+i,l)
        enddo 
      enddo
        
      do iq=1,nqtot
        do l=1,llm
          do i=1,klon
            zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
          enddo
        enddo 
      enddo
        
c      do l=1,llm
c        do i=1,klon
c         pvervel_omp(i,l)=pvervel(offset+i,l)
c       enddo 
c      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
cym      call WriteField_phy_p('zdtfi_omp',zdtfi_omp(:,:),llm)
      
      if (planet_type=="earth") then
#ifdef CPP_EARTH
      CALL physiq (klon,
     .             llm,
     .             debut,
     .             lafin,
     .             jD_cur,
     .             jH_cur,
     .             dtphys,
     .             zplev_omp,
     .             zplay_omp,
     .             zphi_omp,
     .             zphis_omp,
     .             presnivs_omp,
     .             clesphy0,
     .             zufi_omp,
     .             zvfi_omp,
     .             ztfi_omp,
     .             zqfi_omp,
c     .             pvervel_omp,
c#ifdef INCA
     .             flxwfi_omp,
c#endif
     .             zdufi_omp,
     .             zdvfi_omp,
     .             zdtfi_omp,
     .             zdqfi_omp,
     .             zdpsrf_omp,
cIM diagnostique PVteta, Amip2          
     .             pducov,
     .             PVteta)
#endif
      endif !of if (planet_type=="earth")

cym      call WriteField_phy_p('zdtfi_omp',zdtfi_omp(:,:),llm)

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
        
c      do l=1,llm
c        do i=1,klon
c         pvervel(offset+i,l)=pvervel_omp(i,l)
c       enddo 
c      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 
      

cc$OMP END PARALLEL
      klon=klon_mpi
500   CONTINUE
c$OMP BARRIER

c$OMP MASTER
cym      call WriteField_phy('zdtfi',zdtfi(:,:),llm)
      call stop_timer(timer_physic)
c$OMP END MASTER

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

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

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

c$OMP BARRIER


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

c$OMP BARRIER     

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

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

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

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

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

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

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

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

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

C

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

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

      kstart=1
      kend=klon

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


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

      if (is_north_pole) then

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

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

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

        ENDDO
c$OMP END DO NOWAIT

      endif    
      
      if (is_south_pole) then

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

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

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

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

700   CONTINUE
 
      firstcal = .FALSE.

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