!
! $Header$
!
C
C
      SUBROUTINE calfis_p(nq,
     $                  lafin,
     $                  rdayvrai,
     $                  heure,
     $                  pucov,
     $                  pvcov,
     $                  pteta,
     $                  pq,
     $                  pmasse,
     $                  pps,
     $                  pp,
     $                  ppk,
     $                  pphis,
     $                  pphi,
     $                  pducov,
     $                  pdvcov,
     $                  pdteta,
     $                  pdq,
     $                  pw,
#ifdef INCA_CH4
     $                  flxw,
#endif
     $                  clesphy0,
     $                  pdufi,
     $                  pdvfi,
     $                  pdhfi,
     $                  pdqfi,
     $                  pdpsfi)
c
c    Auteur :  P. Le Van, F. Hourdin 
c   .........
      USE dimphy
      USE parallel
      USE Write_Field
      Use Write_field_p
      USE Times
      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"
#include "advtrac.h"

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

#include "comconst.h"
#include "comvert.h"
#include "comgeom2.h"
#include "control.h"
      include 'mpif.h'

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,nqmx)
      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,nqmx)
c
      REAL pw(iip1,jjp1,llm)

      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,nqmx)
      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 zpsrf(klon)
      REAL zplev(klon,llm+1),zplay(klon,llm)
      REAL zphi(klon,llm),zphis(klon)
c
      REAL zufi(klon,llm), zvfi(klon,llm)
      REAL ztfi(klon,llm),zqfi(klon,llm,nqmx)
c
      REAL pcvgu(klon,llm), pcvgv(klon,llm)
      REAL pcvgt(klon,llm), pcvgq(klon,llm,2)
c
      REAL pvervel(klon,llm)
c
      REAL zdufi(klon,llm),zdvfi(klon,llm)
      REAL zdtfi(klon,llm),zdqfi(klon,llm,nqmx)
      REAL zdpsrf(klon)
c
      REAL zsin(iim),zcos(iim),z1(iim)
      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
      REAL unskap, pksurcp

#ifdef INCA_CH4
      REAL flxw(iip1,jjp1,llm)
      REAL flxwfi(klon,llm)
#endif
c
      
      REAL SSUM

      LOGICAL firstcal, debut
      DATA firstcal/.true./
      SAVE firstcal,debut
      REAL rdayvrai
      
      REAL,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 zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm)
      integer :: k,kstart,kend     
c
c-----------------------------------------------------------------------
c
c    1. Initialisations :
c    --------------------
c

      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-----------------------------------------------------------------------
c   latitude, longitude et aires des mailles pour la physique:
c   ----------------------------------------------------------

c
      IF ( firstcal )  THEN
          debut = .TRUE.
      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   ----------------------------------

      call start_timer(timer_physic)
             
      do ig0=1,klon
        i=Liste_i(ig0)
        j=Liste_j(ig0)
        zpsrf(ig0)=pps(i,j)
      enddo



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
      DO l = 1, llmp1
        do ig0=1,klon
          i=Liste_i(ig0)
          j=Liste_j(ig0)
          zplev( ig0,l ) = pp(i,j,l)
        enddo
      ENDDO
c
c

c   43. temperature naturelle (en K) et pressions milieux couches .
c   ---------------------------------------------------------------

      DO l=1,llm

        do ig0=1,klon
          i=Liste_i(ig0)
          j=Liste_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   43.bis traceurs
c   ---------------
c

      DO iq=1,nq
         iiq=niadv(iq)
         DO l=1,llm
           do ig0=1,klon
             i=Liste_i(ig0)
             j=Liste_j(ig0)
             zqfi(ig0,l,iq)  = pq(i,j,l,iiq)
           enddo
         ENDDO
      ENDDO

c   convergence dynamique pour les traceurs "EAU"

      DO iq=1,2
         DO l=1,llm
           do ig0=1,klon
             i=Liste_i(ig0)
             j=Liste_j(ig0)
c             pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
           enddo
         ENDDO
      ENDDO



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)

      DO l=1,llm
	 DO ig=1,klon
	   zphi(ig,l)=zphi(ig,l)-zphis(ig)
	 ENDDO
      ENDDO
      
c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
c
       
      DO l=1,llm
        do ig0=1,klon
           i=Liste_i(ig0)
           j=Liste_j(ig0)
           pvervel(ig0,l) = pw(i,j,l)*g* unsaire(i,j)
        enddo
	if (pole_nord) pvervel(1,l)=pw(1,1,l)*g /apoln
	if (pole_sud) pvervel(klon,l)=pw(1,jjp1,l)*g/apols
      ENDDO


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

      kstart=1
      kend=klon
      
      if (pole_nord) kstart=2
      if (pole_sud) kend=klon-1
      
      DO l=1,llm
        do ig0=kstart,kend
          i=Liste_i(ig0)
          j=Liste_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   46.champ v:
c   -----------

      DO l=1,llm
        DO ig0=kstart,kend
          i=Liste_i(ig0)
          j=Liste_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   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 (pole_nord) then
      
        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
      
      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 (pole_sud) then
      
        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
      
      endif


#ifdef INCA_CH4
      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
#endif


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


      CALL physiq (klon,
     .             llm,
     .             nq,
     .             debut,
     .             lafin,
     .             rdayvrai,
     .             heure,
     .             dtphys,
     .             zplev,
     .             zplay,
     .             zphi,
     .             zphis,
     .             presnivs,
     .             clesphy0,
     .             zufi,
     .             zvfi,
     .             ztfi,
     .             zqfi,
     .             pvervel,
#ifdef INCA_CH4
     .             flxwfi,
#endif
     .             zdufi,
     .             zdvfi,
     .             zdtfi,
     .             zdqfi,
     .             zdpsrf)

500   CONTINUE

      call stop_timer(timer_physic)
      
      if (MPI_rank>0) then
      
        du_send(1:iim,1:llm)=zdufi(1:iim,1:llm)
        dv_send(1:iim,1:llm)=zdvfi(1:iim,1:llm)
        
        call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
     &                   MPI_COMM_WORLD,Req(1),ierr)
        call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
     &                  MPI_COMM_WORLD,Req(2),ierr)
     
      endif
   
      if (MPI_rank<MPI_Size-1) then
      
        call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
     &                 MPI_COMM_WORLD,Req(3),ierr)
        call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
     &                 MPI_COMM_WORLD,Req(4),ierr)
     
      endif
   
      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
      
      zdufi2(1:klon,:)=zdufi(1:klon,:)
      zdufi2(klon+1:klon+iim,:)=du_recv(1:iim,:)
          
      zdvfi2(1:klon,:)=zdvfi(1:klon,:)
      zdvfi2(klon+1:klon+iim,:)=dv_recv(1:iim,:)

       pdhfi(:,jjphy_begin,:)=0
       pdqfi(:,jjphy_begin,:,:)=0
       pdufi(:,jjphy_begin,:)=0
       pdvfi(:,jjphy_begin,:)=0
       pdpsfi(:,jjphy_begin)=0

       if (.not. pole_sud) then
         pdhfi(:,jjphy_end,:)=0
         pdqfi(:,jjphy_end,:,:)=0
         pdufi(:,jjphy_end,:)=0
         pdvfi(:,jjphy_end,:)=0
	 pdpsfi(:,jjphy_end)=0
       endif

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 (pole_nord) kstart=2
      if (pole_sud)  kend=klon-1

      DO l=1,llm

!cdir NODEP
        do ig0=kstart,kend
          i=Liste_i(ig0)
          j=Liste_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 (pole_nord) then
            DO i=1,iip1
              pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
            enddo
        endif
        
        if (pole_sud) then
            DO i=1,iip1
              pdhfi(i,jjp1,l) = cpp *  zdtfi(klon,l)/ ppk(i,jjp1,l)
            ENDDO
        endif
      ENDDO
      
c   62. humidite specifique
c   ---------------------

      DO iq=1,nqmx
         DO l=1,llm
!cdir NODEP 
           do ig0=kstart,kend
             i=Liste_i(ig0)
             j=Liste_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 (pole_nord) then
             do i=1,iip1
               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)             
             enddo
           endif
           
           if (pole_sud) then
             do i=1,iip1
               pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq) 
             enddo
           endif
           
         ENDDO
      ENDDO

c   63. traceurs
c   ------------
C     initialisation des tendances
      pdqfi=0.
C

      DO iq=1,nq
         iiq=niadv(iq)
         DO l=1,llm

!cdir NODEP           
	     DO ig0=kstart,kend
              i=Liste_i(ig0)
              j=Liste_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 (pole_nord) then
	      DO i=1,iip1
                pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
	      ENDDO
	    ENDIF
	    
	    IF (pole_sud) then
	      DO i=1,iip1
                pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
	      ENDDO
	    ENDIF
	    
         ENDDO
      ENDDO
      
c   65. champ u:
c   ------------

      DO l=1,llm
!cdir NODEP
         do ig0=kstart,kend
           i=Liste_i(ig0)
           j=Liste_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 (Pole_nord) then
           DO i=1,iip1
            pdufi(i,1,l)    = 0.
           ENDDO
         endif
         
         if (Pole_sud) then
           DO i=1,iip1
            pdufi(i,jjp1,l) = 0.
           ENDDO
         endif
         
      ENDDO


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

      kstart=1
      kend=klon

      if (pole_nord) kstart=2
      if (pole_sud)  kend=klon-1-iim
      
      DO l=1,llm
!cdir NODEP
        do ig0=kstart,kend
           i=Liste_i(ig0)
           j=Liste_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   68. champ v pres des poles:
c   ---------------------------
c      v = U * cos(long) + V * SIN(long)

      if (pole_nord) then
        
        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

      endif    
      
      if (pole_sud) then
      
        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
     
      endif
c-----------------------------------------------------------------------

700   CONTINUE
 
      firstcal = .FALSE.

      RETURN
      END
