!
! $Id: calfis.F 1225 2009-08-07 16:29:23Z oboucher $
!
C
C
      SUBROUTINE calfis(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)
c
c    Auteur :  P. Le Van, F. Hourdin 
c   .........
      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       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"

c    Arguments :
c    -----------
      LOGICAL  lafin


      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 zpsrf(ngridmx)
      REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
      REAL zphi(ngridmx,llm),zphis(ngridmx)
c
      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
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 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(ngridmx,ntetaSTD)
c
      REAL flxw(iip1,jjp1,llm)  ! Flux de masse verticale sur la grille dynamique
      REAL flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
c
      
      REAL SSUM

      LOGICAL firstcal, debut
      DATA firstcal/.true./
      SAVE firstcal,debut
!      REAL rdayvrai
      REAL, intent(in):: jD_cur, jH_cur
c
c-----------------------------------------------------------------------
c
c    1. Initialisations :
c    --------------------
c
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
      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 :
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
        zplev( 1,l ) = pp(1,1,l)
        ig0 = 2
          DO j = 2, jjm
             DO i =1, iim
              zplev( ig0,l ) = pp(i,j,l)
              ig0 = ig0 +1
             ENDDO
          ENDDO
        zplev( ngridmx,l ) = pp(1,jjp1,l)
      ENDDO
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
      DO iq=1,nqtot
          iiq=niadv(iq) 
         DO l=1,llm
            zqfi(1,l,iq) = pq(1,1,l,iiq)
            ig0          = 2
            DO j=2,jjm
               DO i = 1, iim
                  zqfi(ig0,l,iq)  = pq(i,j,l,iiq)
                  ig0             = ig0 + 1
               ENDDO
            ENDDO
            zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
         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   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
         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

      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

      ENDDO
c
      if (planet_type=="earth") then
#ifdef CPP_EARTH
cIM calcul PV a teta=350, 380, 405K
      CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
     $           ztfi,zplay,zplev,
     $           ntetaSTD,rtetaSTD,PVteta)
#endif
      endif
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   ---------------------


      if (planet_type=="earth") then
#ifdef CPP_EARTH
      CALL physiq (ngridmx,
     .             llm,
     .             debut,
     .             lafin,
     .             jD_cur,
     .             jH_cur,
     .             dtphys,
     .             zplev,
     .             zplay,
     .             zphi,
     .             zphis,
     .             presnivs,
     .             clesphy0,
     .             zufi,
     .             zvfi,
     .             ztfi,
     .             zqfi,
     .             flxwfi,
     .             zdufi,
     .             zdvfi,
     .             zdtfi,
     .             zdqfi,
     .             zdpsrf,
cIM diagnostique PVteta, Amip2          
     .             pducov,
     .             PVteta)
#endif
      endif !of if (planet_type=="earth")

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
      DO iq=1,nqtot
         iiq=niadv(iq)
         DO l=1,llm
            DO i=1,iip1
               pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
               pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
            ENDDO
            DO j=2,jjm
               ig0=1+(j-2)*iim
               DO i=1,iim
                  pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
               ENDDO
               pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
            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