source: LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/calfis.f90 @ 5182

! $Id: calfis.f90 5182 2024-09-10 14:25:29Z 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 lmdz_infotrac, ONLY: nqtot, tracers
  USE control_mod, ONLY: planet_type, nsplit_phys
  USE callphysiq_mod, ONLY: call_physiq
  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS
  USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
  USE comvert_mod, ONLY: preff, presnivs
  USE lmdz_iniprint, ONLY: lunout, prt_level
  USE lmdz_ssum_scopy, ONLY: scopy, ssum
  USE lmdz_comgeom2
  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
  USE lmdz_paramet

  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)

  !    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

  LOGICAL, SAVE :: firstcal = .TRUE., debut = .TRUE.
  ! REAL rdayvrai


  !-----------------------------------------------------------------------

  !    1. Initialisations :
  !    --------------------


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

  zdufi(:, :) = zdufic(:, :) / nsplit_phys
  zdvfi(:, :) = zdvfic(:, :) / nsplit_phys
  zdtfi(:, :) = zdtfic(:, :) / nsplit_phys
  zdqfi(:, :, :) = zdqfic(:, :, :) / nsplit_phys

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

  !-----------------------------------------------------------------------
  firstcal = .FALSE.

  RETURN
END SUBROUTINE calfis