source: LMDZ6/trunk/libf/dynphy_lonlat/calfis.f90 @ 5272

!
! $Id: calfis.f90 5272 2024-10-24 15:53:15Z abarral $
!
!
!
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
  USE callphysiq_mod, ONLY: call_physiq
  USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
  USE comvert_mod, ONLY: preff, presnivs
  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS
  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
  USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
          ip1jmi1, ijp1llm, ijmllm, mvar, jcfil, jcfllm
  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 :
  !    ------------------
  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.

IF (CPPKEY_PHYS) THEN

   do isplit=1,nsplit_phys

     jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
     debut_split=debut.and.isplit==1
     lafin_split=lafin.and.isplit==nsplit_phys

   ! 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

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