! $Id: leapfrog.F90 5103 2024-07-23 13:29:36Z abarral $

!
!
SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)


  !IM : pour sortir les param. du modele dans un fis. netcdf 110106
  use IOIPSL
  USE infotrac, ONLY: nqtot, isoCheck
  USE guide_mod, ONLY: guide_main
  USE write_field, ONLY: writefield
  USE control_mod, ONLY: nday, day_step, planet_type, offline, &
          iconser, iphysiq, iperiod, dissip_period, &
          iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins, &
          periodav, ok_dyn_ave, output_grads_dyn
  use exner_hyb_m, only: exner_hyb
  use exner_milieu_m, only: exner_milieu
  USE comvert_mod, ONLY: ap, bp, pressure_exner, presnivs
  USE comconst_mod, ONLY: cpp, dtphys, dtvr, pi, ihf
  USE logic_mod, ONLY: iflag_phys, ok_guide, forward, leapf, apphys, &
          statcl, conser, apdiss, purmats, ok_strato
  USE temps_mod, ONLY: jD_ref, jH_ref, itaufin, day_ini, day_ref, &
          start_time, dt
  USE strings_mod, ONLY: msg
  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS

  IMPLICIT NONE

  ! ......   Version  du 10/01/98    ..........

  !        avec  coordonnees  verticales hybrides
  !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )

  !=======================================================================
  !
  !   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
  !   -------
  !
  !   Objet:
  !   ------
  !
  !   GCM LMD nouvelle grille
  !
  !=======================================================================
  !
  !  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
  !  et possibilite d'appeler une fonction f(y)  a derivee tangente
  !  hyperbolique a la  place de la fonction a derivee sinusoidale.

  !  ... Possibilite de choisir le shema pour l'advection de
  !    q  , en modifiant iadv dans traceur.def  (10/02) .
  !
  !  Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
  !  Pour Van-Leer iadv=10
  !
  !-----------------------------------------------------------------------
  !   Declarations:
  !   -------------

  include "dimensions.h"
  include "paramet.h"
  include "comdissnew.h"
  include "comgeom.h"
  include "description.h"
  include "iniprint.h"
  include "academic.h"

  REAL, INTENT(IN) :: time_0 ! not used

  !   dynamical variables:
  REAL, INTENT(INOUT) :: ucov(ip1jmp1, llm)    ! zonal covariant wind
  REAL, INTENT(INOUT) :: vcov(ip1jm, llm)      ! meridional covariant wind
  REAL, INTENT(INOUT) :: teta(ip1jmp1, llm)    ! potential temperature
  REAL, INTENT(INOUT) :: ps(ip1jmp1)          ! surface pressure (Pa)
  REAL, INTENT(INOUT) :: masse(ip1jmp1, llm)   ! air mass
  REAL, INTENT(INOUT) :: phis(ip1jmp1)        ! geopotentiat at the surface
  REAL, INTENT(INOUT) :: q(ip1jmp1, llm, nqtot) ! advected tracers

  REAL :: p (ip1jmp1, llmp1)               ! interlayer pressure
  REAL :: pks(ip1jmp1)                      ! exner at the surface
  REAL :: pk(ip1jmp1, llm)                   ! exner at mid-layer
  REAL :: pkf(ip1jmp1, llm)                  ! filtered exner at mid-layer
  REAL :: phi(ip1jmp1, llm)                  ! geopotential
  REAL :: w(ip1jmp1, llm)                    ! vertical velocity

  real :: zqmin, zqmax

  ! variables dynamiques intermediaire pour le transport
  REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse

  !   variables dynamiques au pas -1
  REAL :: vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)
  REAL :: tetam1(ip1jmp1, llm), psm1(ip1jmp1)
  REAL :: massem1(ip1jmp1, llm)

  !   tendances dynamiques
  REAL :: dv(ip1jm, llm), du(ip1jmp1, llm)
  REAL :: dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqtot), dp(ip1jmp1)

  !   tendances de la dissipation
  REAL :: dvdis(ip1jm, llm), dudis(ip1jmp1, llm)
  REAL :: dtetadis(ip1jmp1, llm)

  !   tendances physiques
  REAL :: dvfi(ip1jm, llm), dufi(ip1jmp1, llm)
  REAL :: dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqtot), dpfi(ip1jmp1)

  !   variables pour le fichier histoire
  REAL :: dtav      ! intervalle de temps elementaire

  REAL :: tppn(iim), tpps(iim), tpn, tps
  !
  INTEGER :: itau, itaufinp1, iav
  ! INTEGER  iday ! jour julien
  REAL :: time

  REAL :: SSUM
  ! REAL finvmaold(ip1jmp1,llm)

  !ym      LOGICAL  lafin
  LOGICAL :: lafin = .FALSE.
  INTEGER :: ij, iq, l
  INTEGER :: ik

  real :: time_step, t_wrt, t_ops

  ! REAL rdayvrai,rdaym_ini
  ! jD_cur: jour julien courant
  ! jH_cur: heure julienne courante
  REAL :: jD_cur, jH_cur
  INTEGER :: an, mois, jour
  REAL :: secondes

  LOGICAL :: first, callinigrads
  !IM : pour sortir les param. du modele dans un fis. netcdf 110106
  save first
  data first/.TRUE./
  real :: dt_cum
  character(len = 10) :: infile
  integer :: zan, tau0, thoriid
  integer :: nid_ctesGCM
  save nid_ctesGCM
  real :: degres
  real :: rlong(iip1), rlatg(jjp1)
  real :: zx_tmp_2d(iip1, jjp1)
  integer :: ndex2d(iip1 * jjp1)
  logical :: ok_sync
  parameter (ok_sync = .TRUE.)
  logical :: physic

  data callinigrads/.TRUE./
  character(len = 10) :: string10

  REAL :: flxw(ip1jmp1, llm)  ! flux de masse verticale

  !+jld variables test conservation energie
  REAL :: ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)
  ! Tendance de la temp. potentiel d (theta)/ d t due a la
  ! tansformation d'energie cinetique en energie thermique
  ! cree par la dissipation
  REAL :: dtetaecdt(ip1jmp1, llm)
  REAL :: vcont(ip1jm, llm), ucont(ip1jmp1, llm)
  REAL :: vnat(ip1jm, llm), unat(ip1jmp1, llm)
  REAL :: d_h_vcol, d_qt, d_qw, d_ql, d_ec
  CHARACTER(len = 15) :: ztit
  !IM   INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
  !IM   SAVE      ip_ebil_dyn
  !IM   DATA      ip_ebil_dyn/0/
  !-jld

  character(len = 80) :: dynhist_file, dynhistave_file
  character(len = *), parameter :: modname = "leapfrog"
  character(len = 80) :: abort_message

  logical :: dissip_conservative
  save dissip_conservative
  data dissip_conservative/.TRUE./

  LOGICAL :: prem
  save prem
  DATA prem/.TRUE./
  INTEGER :: testita
  PARAMETER (testita = 9)

  logical, parameter :: flag_verif = .FALSE.

  integer :: itau_w   ! pas de temps ecriture = itap + itau_phy

  if (nday>=0) then
    itaufin = nday * day_step
  else
    itaufin = -nday
  endif
  itaufinp1 = itaufin + 1
  itau = 0
  physic = .TRUE.
  if (iflag_phys==0.or.iflag_phys==2) physic = .FALSE.

  ! iday = day_ini+itau/day_step
  ! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
  !    IF(time.GT.1.) THEN
  !     time = time-1.
  !     iday = iday+1
  !    ENDIF


  !-----------------------------------------------------------------------
  !   On initialise la pression et la fonction d'Exner :
  !   --------------------------------------------------

  dq(:, :, :) = 0.
  CALL pression (ip1jmp1, ap, bp, ps, p)
  if (pressure_exner) then
    CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
  else
    CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
  endif

  !-----------------------------------------------------------------------
  !   Debut de l'integration temporelle:
  !   ----------------------------------

  1   CONTINUE ! Matsuno Forward step begins here

  !   date: (NB: date remains unchanged for Backward step)
  !   -----

  jD_cur = jD_ref + day_ini - day_ref +                             & &
  (itau+1)/day_step
  jH_cur = jH_ref + start_time +                                    & &
          mod(itau+1, day_step)/float(day_step)
  jD_cur = jD_cur + int(jH_cur)
  jH_cur = jH_cur - int(jH_cur)

  CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 321')

  if (ok_guide) then
    CALL guide_main(itau,ucov,vcov,teta,q,masse,ps)
  endif


  !
  ! IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
  !   CALL  test_period ( ucov,vcov,teta,q,p,phis )
  !   PRINT *,' ----   Test_period apres continue   OK ! -----', itau
  ! ENDIF
  !

  ! Save fields obtained at previous time step as '...m1'
  CALL SCOPY(ijmllm, vcov, 1, vcovm1, 1)
  CALL SCOPY(ijp1llm, ucov, 1, ucovm1, 1)
  CALL SCOPY(ijp1llm, teta, 1, tetam1, 1)
  CALL SCOPY(ijp1llm, masse, 1, massem1, 1)
  CALL SCOPY(ip1jmp1, ps, 1, psm1, 1)

  forward = .TRUE.
  leapf = .FALSE.
  dt = dtvr

  !   ...    P.Le Van .26/04/94  ....
  ! Ehouarn: finvmaold is actually not used
  ! CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
  ! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )

  CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 400')

  2   CONTINUE ! Matsuno backward or leapfrog step begins here

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

  !   date: (NB: only leapfrog step requires recomputing date)
  !   -----

  IF (leapf) THEN
    jD_cur = jD_ref + day_ini - day_ref + &
            (itau + 1) / day_step
    jH_cur = jH_ref + start_time + &
            mod(itau + 1, day_step) / float(day_step)
    jD_cur = jD_cur + int(jH_cur)
    jH_cur = jH_cur - int(jH_cur)
  ENDIF


  !   gestion des appels de la physique et des dissipations:
  !   ------------------------------------------------------
  !
  !   ...    P.Le Van  ( 6/02/95 )  ....

  apphys = .FALSE.
  statcl = .FALSE.
  conser = .FALSE.
  apdiss = .FALSE.

  IF(purmats) THEN
    ! ! Purely Matsuno time stepping
    IF(MOD(itau, iconser) ==0.AND.  forward) conser = .TRUE.
    IF(MOD(itau, dissip_period)==0.AND..NOT.forward) &
            apdiss = .TRUE.
    IF(MOD(itau, iphysiq)==0.AND..NOT.forward &
            .and. physic) apphys = .TRUE.
  ELSE
    ! ! Leapfrog/Matsuno time stepping
    IF(MOD(itau, iconser) == 0) conser = .TRUE.
    IF(MOD(itau + 1, dissip_period)==0 .AND. .NOT. forward) &
            apdiss = .TRUE.
    IF(MOD(itau + 1, iphysiq)==0.AND.physic) apphys = .TRUE.
  END IF

  ! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
  ! supress dissipation step
  if (llm==1) then
    apdiss = .FALSE.
  endif

  CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 589')

  !-----------------------------------------------------------------------
  !   calcul des tendances dynamiques:
  !   --------------------------------

  ! ! compute geopotential phi()
  CALL geopot  (ip1jmp1, teta, pk, pks, phis, phi)

  time = jD_cur + jH_cur
  CALL caldyn &
          (itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
          phi, conser, du, dv, dteta, dp, w, pbaru, pbarv, time)


  !-----------------------------------------------------------------------
  !   calcul des tendances advection des traceurs (dont l'humidite)
  !   -------------------------------------------------------------

  CALL check_isotopes_seq(q, ip1jmp1, &
          'leapfrog 686: avant caladvtrac')

  IF(forward .OR.  leapf)  THEN
    ! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
    CALL caladvtrac(q, pbaru, pbarv, &
            p, masse, dq, teta, &
            flxw, pk)
    ! !write(*,*) 'caladvtrac 346'

    IF (offline) THEN
      !maf stokage du flux de masse pour traceurs OFF-LINE

       CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis, &
             dtvr, itau)


    ENDIF ! of IF (offline)
    !
  ENDIF ! of IF( forward .OR.  leapf )


  !-----------------------------------------------------------------------
  !   integrations dynamique et traceurs:
  !   ----------------------------------

  CALL msg('720', modname, isoCheck)
  CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 756')

  CALL integrd (nqtot, vcovm1, ucovm1, tetam1, psm1, massem1, &
          dv, du, dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis)
  ! $              finvmaold                                    )

  CALL msg('724', modname, isoCheck)
  CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 762')

  ! .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
  !
  !-----------------------------------------------------------------------
  !   calcul des tendances physiques:
  !   -------------------------------
  !    ########   P.Le Van ( Modif le  6/02/95 )   ###########
  !
  IF(purmats)  THEN
    IF(itau==itaufin.AND..NOT.forward) lafin = .TRUE.
  ELSE
    IF(itau + 1 == itaufin)              lafin = .TRUE.
  ENDIF
  !
  !
  IF(apphys)  THEN
    !
    ! .......   Ajout   P.Le Van ( 17/04/96 )   ...........
    !

    CALL pression (ip1jmp1, ap, bp, ps, p)
    if (pressure_exner) then
      CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
    else
      CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
    endif

    ! Appel a geopot ajoute le 2014/05/08 pour garantir la convergence numerique
    ! avec dyn3dmem
    CALL geopot  (ip1jmp1, teta, pk, pks, phis, phi)

    ! rdaym_ini  = itau * dtvr / daysec
    ! rdayvrai   = rdaym_ini  + day_ini
    ! jD_cur = jD_ref + day_ini - day_ref
    ! $        + int (itau * dtvr / daysec)
    !       jH_cur = jH_ref +                                            &
    ! &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
    jD_cur = jD_ref + day_ini - day_ref +                        & &
    (itau+1)/day_step

    IF (planet_type =="generic") THEN
      ! ! AS: we make jD_cur to be pday
      jD_cur = int(day_ini + itau / day_step)
    ENDIF

    jH_cur = jH_ref + start_time +                               & &
            mod(itau+1, day_step)/float(day_step)
    jD_cur = jD_cur + int(jH_cur)
    jH_cur = jH_cur - int(jH_cur)
    ! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
    ! CALL ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
    ! write(lunout,*)'current date = ',an, mois, jour, secondes

    ! rajout debug
    ! lafin = .TRUE.


    !   Inbterface avec les routines de phylmd (phymars ... )
    !   -----------------------------------------------------

    !+jld

    !  Diagnostique de conservation de l'energie : initialisation
    IF (ip_ebil_dyn>=1) THEN
      ztit = 'bil dyn'
      ! Ehouarn: be careful, diagedyn is Earth-specific!
      IF (planet_type=="earth") THEN
        CALL diagedyn(ztit, 2, 1, 1, dtphys &
                , ucov, vcov, ps, p, pk, teta, q(:, :, 1), q(:, :, 2))
      ENDIF
    ENDIF ! of IF (ip_ebil_dyn.ge.1 )
    IF (CPPKEY_PHYS) THEN
      CALL calfis(lafin, jD_cur, jH_cur, &
              ucov, vcov, teta, q, masse, ps, p, pk, phis, phi, &
              du, dv, dteta, dq, &
              flxw, dufi, dvfi, dtetafi, dqfi, dpfi)
    END IF
    ! ajout des tendances physiques:
    ! ------------------------------
    CALL addfi(dtphys, leapf, forward, &
            ucov, vcov, teta, q, ps, &
            dufi, dvfi, dtetafi, dqfi, dpfi)
    ! ! since addfi updates ps(), also update p(), masse() and pk()
    CALL pression (ip1jmp1, ap, bp, ps, p)
    CALL massdair(p, masse)
    if (pressure_exner) then
      CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
    else
      CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
    endif

    IF (ok_strato) THEN
      CALL top_bound(vcov, ucov, teta, masse, dtphys)
    ENDIF

    !
    !  Diagnostique de conservation de l'energie : difference
    IF (ip_ebil_dyn>=1) THEN
      ztit = 'bil phys'
      IF (planet_type=="earth") THEN
        CALL diagedyn(ztit, 2, 1, 1, dtphys &
                , ucov, vcov, ps, p, pk, teta, q(:, :, 1), q(:, :, 2))
      ENDIF
    ENDIF ! of IF (ip_ebil_dyn.ge.1 )

  ENDIF ! of IF( apphys )

  IF(iflag_phys==2) THEN ! "Newtonian" case
    !   Academic case : Simple friction and Newtonan relaxation
    !   -------------------------------------------------------
    DO l = 1, llm
      DO ij = 1, ip1jmp1
        teta(ij, l) = teta(ij, l) - dtvr * &
                (teta(ij, l) - tetarappel(ij, l)) * (knewt_g + knewt_t(l) * clat4(ij))
      ENDDO
    ENDDO ! of DO l=1,llm

    if (planet_type=="giant") then
      ! ! add an intrinsic heat flux at the base of the atmosphere
      teta(:, 1) = teta(:, 1) + dtvr * aire(:) * ihf / cpp / masse(:, 1)
    endif

    CALL friction(ucov, vcov, dtvr)

    ! ! Sponge layer (if any)
    IF (ok_strato) THEN
      ! dufi(:,:)=0.
      ! dvfi(:,:)=0.
      ! dtetafi(:,:)=0.
      ! dqfi(:,:,:)=0.
      !          dpfi(:)=0.
      ! CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
      CALL top_bound(vcov, ucov, teta, masse, dtvr)
      ! CALL addfi( dtvr, leapf, forward   ,
      ! $                  ucov, vcov, teta , q   ,ps ,
      ! $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
    ENDIF ! of IF (ok_strato)
  ENDIF ! of IF (iflag_phys.EQ.2)


  !-jld

  CALL pression (ip1jmp1, ap, bp, ps, p)
  if (pressure_exner) then
    CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
  else
    CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
  endif
  CALL massdair(p, masse)

  CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1196')

  !-----------------------------------------------------------------------
  !   dissipation horizontale et verticale  des petites echelles:
  !   ----------------------------------------------------------

  IF(apdiss) THEN


    !   calcul de l'energie cinetique avant dissipation
    CALL covcont(llm, ucov, vcov, ucont, vcont)
    CALL enercin(vcov, ucov, vcont, ucont, ecin0)

    !   dissipation
    CALL dissip(vcov, ucov, teta, p, dvdis, dudis, dtetadis)
    ucov = ucov + dudis
    vcov = vcov + dvdis
    ! teta=teta+dtetadis


    !------------------------------------------------------------------------
    if (dissip_conservative) then
      ! On rajoute la tendance due a la transform. Ec -> E therm. cree
      ! lors de la dissipation
      CALL covcont(llm, ucov, vcov, ucont, vcont)
      CALL enercin(vcov, ucov, vcont, ucont, ecin)
      dtetaecdt = (ecin0 - ecin) / pk
      ! teta=teta+dtetaecdt
      dtetadis = dtetadis + dtetaecdt
    endif
    teta = teta + dtetadis
    !------------------------------------------------------------------------


    !    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
    !   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
    !

    DO l = 1, llm
      DO ij = 1, iim
        tppn(ij) = aire(ij) * teta(ij, l)
        tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)
      ENDDO
      tpn = SSUM(iim, tppn, 1) / apoln
      tps = SSUM(iim, tpps, 1) / apols

      DO ij = 1, iip1
        teta(ij, l) = tpn
        teta(ij + ip1jm, l) = tps
      ENDDO
    ENDDO

    if (1 == 0) then
      !!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
      !!!                     2) should probably not be here anyway
      !!! but are kept for those who would want to revert to previous behaviour
      DO ij = 1, iim
        tppn(ij) = aire(ij) * ps (ij)
        tpps(ij) = aire(ij + ip1jm) * ps (ij + ip1jm)
      ENDDO
      tpn = SSUM(iim, tppn, 1) / apoln
      tps = SSUM(iim, tpps, 1) / apols

      DO ij = 1, iip1
        ps(ij) = tpn
        ps(ij + ip1jm) = tps
      ENDDO
    endif ! of if (1 == 0)

  END IF ! of IF(apdiss)

  ! ajout debug
  ! IF( lafin ) then
  !   abort_message = 'Simulation finished'
  !   CALL abort_gcm(modname,abort_message,0)
  ! ENDIF

  !   ********************************************************************
  !   ********************************************************************
  !   .... fin de l'integration dynamique  et physique pour le pas itau ..
  !   ********************************************************************
  !   ********************************************************************

  !   preparation du pas d'integration suivant  ......

  CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1509')

  IF (.NOT.purmats) THEN
    ! ........................................................
    ! ..............  schema matsuno + leapfrog  ..............
    ! ........................................................

    IF(forward .OR. leapf) THEN
      itau = itau + 1
      ! iday= day_ini+itau/day_step
      ! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
      !   IF(time.GT.1.) THEN
      !     time = time-1.
      !     iday = iday+1
      !   ENDIF
    ENDIF

    IF(itau == itaufinp1) then
      if (flag_verif) then
        write(79, *) 'ucov', ucov
        write(80, *) 'vcov', vcov
        write(81, *) 'teta', teta
        write(82, *) 'ps', ps
        write(83, *) 'q', q
        WRITE(85, *) 'q1 = ', q(:, :, 1)
        WRITE(86, *) 'q3 = ', q(:, :, 3)
      endif

      abort_message = 'Simulation finished'

      CALL abort_gcm(modname, abort_message, 0)
    ENDIF
    !-----------------------------------------------------------------------
    !   ecriture du fichier histoire moyenne:
    !   -------------------------------------

    IF(MOD(itau, iperiod)==0 .OR. itau==itaufin) THEN
      IF(itau==itaufin) THEN
        iav = 1
      ELSE
        iav = 0
      ENDIF

      ! ! Ehouarn: re-compute geopotential for outputs
      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)

      IF (ok_dynzon) THEN
             CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav, &
                   ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
      END IF
      IF (ok_dyn_ave) THEN
             CALL writedynav(itau,vcov, &
                   ucov,teta,pk,phi,q,masse,ps,phis)
      ENDIF

    ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))

    CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1584')

    !-----------------------------------------------------------------------
    !   ecriture de la bande histoire:
    !   ------------------------------

    IF(MOD(itau, iecri)==0) THEN
      ! ! Ehouarn: output only during LF or Backward Matsuno
      if (leapf.or.(.not.leapf.and.(.not.forward))) then
        CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
        unat = 0.
        do l = 1, llm
          unat(iip2:ip1jm, l) = ucov(iip2:ip1jm, l) / cu(iip2:ip1jm)
          vnat(:, l) = vcov(:, l) / cv(:)
        enddo
          if (ok_dyn_ins) then
            ! write(lunout,*) "leapfrog: CALL writehist, itau=",itau
           CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
            ! CALL WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
            ! CALL WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
           ! CALL WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
           !  CALL WriteField('ps',reshape(ps,(/iip1,jmp1/)))
           !  CALL WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
          endif ! of if (ok_dyn_ins)
        ! For some Grads outputs of fields
        if (output_grads_dyn) then
          include "write_grads_dyn.h"
        endif
      endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
    ENDIF ! of IF(MOD(itau,iecri).EQ.0)

    IF(itau==itaufin) THEN


      ! if (planet_type.eq."earth") then
      ! Write an Earth-format restart file
      CALL dynredem1("restart.nc", start_time, &
              vcov, ucov, teta, q, masse, ps)
      ! endif ! of if (planet_type.eq."earth")

      CLOSE(99)
      if (ok_guide) then
        ! ! set ok_guide to false to avoid extra output
        ! ! in following forward step
        ok_guide = .FALSE.
      endif
      ! !!! Ehouarn: Why not stop here and now?
    ENDIF ! of IF (itau.EQ.itaufin)

    !-----------------------------------------------------------------------
    !   gestion de l'integration temporelle:
    !   ------------------------------------

    IF(MOD(itau, iperiod)==0)    THEN
      GO TO 1
    ELSE IF (MOD(itau - 1, iperiod) == 0) THEN

      IF(forward)  THEN
        ! fin du pas forward et debut du pas backward

        forward = .FALSE.
        leapf = .FALSE.
        GO TO 2

      ELSE
        ! fin du pas backward et debut du premier pas leapfrog

        leapf = .TRUE.
        dt = 2. * dtvr
        GO TO 2
      END IF ! of IF (forward)
    ELSE

      ! ......   pas leapfrog  .....

      leapf = .TRUE.
      dt = 2. * dtvr
      GO TO 2
    END IF ! of IF (MOD(itau,iperiod).EQ.0)
    ! !    ELSEIF (MOD(itau-1,iperiod).EQ.0)

  ELSE ! of IF (.not.purmats)

    CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1664')

    ! ........................................................
    ! ..............       schema  matsuno        ...............
    ! ........................................................
    IF(forward)  THEN

      itau = itau + 1
      ! iday = day_ini+itau/day_step
      ! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
      !
      !              IF(time.GT.1.) THEN
      !               time = time-1.
      !               iday = iday+1
      !              ENDIF

      forward = .FALSE.
      IF(itau == itaufinp1) then
        abort_message = 'Simulation finished'
        CALL abort_gcm(modname, abort_message, 0)
      ENDIF
      GO TO 2

    ELSE ! of IF(forward) i.e. backward step

      CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1698')

      IF(MOD(itau, iperiod)==0 .OR. itau==itaufin) THEN
        IF(itau==itaufin) THEN
          iav = 1
        ELSE
          iav = 0
        ENDIF

        ! ! Ehouarn: re-compute geopotential for outputs
        CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)

        IF (ok_dynzon) THEN
             CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav, &
                   ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
        ENDIF
        IF (ok_dyn_ave) THEN
             CALL writedynav(itau,vcov, &
                   ucov,teta,pk,phi,q,masse,ps,phis)
        ENDIF

      ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)

      IF(MOD(itau, iecri)==0) THEN
        ! IF(MOD(itau,iecri*day_step).EQ.0) THEN
        CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
        unat = 0.
        do l = 1, llm
          unat(iip2:ip1jm, l) = ucov(iip2:ip1jm, l) / cu(iip2:ip1jm)
          vnat(:, l) = vcov(:, l) / cv(:)
        enddo
          if (ok_dyn_ins) then
             ! write(lunout,*) "leapfrog: CALL writehist (b)",
  ! &                        itau,iecri
            CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
          endif ! of if (ok_dyn_ins)
        ! For some Grads outputs
        if (output_grads_dyn) then
          include "write_grads_dyn.h"
        endif

      ENDIF ! of IF(MOD(itau,iecri         ).EQ.0)

      IF(itau==itaufin) THEN
        ! if (planet_type.eq."earth") then
        CALL dynredem1("restart.nc", start_time, &
                vcov, ucov, teta, q, masse, ps)
        ! endif ! of if (planet_type.eq."earth")
        if (ok_guide) then
          ! ! set ok_guide to false to avoid extra output
          ! ! in following forward step
          ok_guide = .FALSE.
        endif
      ENDIF ! of IF(itau.EQ.itaufin)

      forward = .TRUE.
      GO TO  1

    ENDIF ! of IF (forward)

  END IF ! of IF(.not.purmats)

END SUBROUTINE leapfrog
