source: LMDZ6/branches/Amaury_dev/libf/dyn3d/lmdz_leapfrog.f90 @ 5214

MODULE lmdz_leapfrog

  IMPLICIT NONE; PRIVATE
  PUBLIC leapfrog

CONTAINS

  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 lmdz_infotrac, ONLY: nqtot, isoCheck
    USE guide_mod, ONLY: guide_main
    USE lmdz_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 lmdz_strings, ONLY: msg
    USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS
    USE lmdz_iniprint, ONLY: lunout, prt_level
    USE lmdz_ssum_scopy, ONLY: scopy, ssum
    USE lmdz_academic, ONLY: tetarappel, knewt_t, kfrict, knewt_g, clat4
    USE lmdz_comdissnew, ONLY: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
            tetagrot, tetatemp, coefdis, vert_prof_dissip
    USE lmdz_comgeom

    USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
    USE lmdz_paramet
    USE lmdz_addfi, ONLY: addfi
    USE lmdz_caldyn, ONLY: caldyn
    USE lmdz_caladvtrac, ONLY: caladvtrac
    USE lmdz_bilan_dyn, ONLY: bilan_dyn
    USE lmdz_check_isotopes, ONLY: check_isotopes_seq
    USE lmdz_writedynav, ONLY: writedynav
    USE lmdz_writehist, ONLY: writehist
    USE lmdz_dissip, ONLY: dissip
    USE lmdz_dynredem, ONLY: dynredem1

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

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

END MODULE lmdz_leapfrog