source: LMDZ6/trunk/libf/dyn3d/leapfrog.f90 @ 5431

!
! $Id: leapfrog.f90 5312 2024-11-01 12:22:08Z aborella $
!
!
!
SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
  !IM : pour sortir les param. du modele dans un fis. netcdf 110106
  USE iniprint_mod_h
  USE comgeom_mod_h
  USE comdissnew_mod_h
  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
  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
  USE paramet_mod_h
  USE academic_mod_h, ONLY: tetarappel, knewt_t, knewt_g, clat4
  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 :: 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) .EQ.0.AND.  forward    ) conser = .TRUE.
     IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward ) &
           apdiss = .TRUE.
     IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward &
           .and. physic                        ) apphys = .TRUE.
  ELSE
  ! ! Leapfrog/Matsuno time stepping
     IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
     IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward ) &
           apdiss = .TRUE.
     IF( MOD(itau+1,iphysiq).EQ.0.AND.physic       ) apphys=.TRUE.
  END IF

  ! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
       ! supress dissipation step
  if (llm.eq.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.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
   ELSE
      IF( itau+1.EQ. 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 .eq."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.ge.1 ) THEN
      ztit='bil dyn'
  ! Ehouarn: be careful, diagedyn is Earth-specific!
       IF (planet_type.eq."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 )
  !-jld

  !IM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
  !IM uncomment next 6 lines to get some parameters for LMDZ dynamics
     ! IF (first) THEN
     !  first=.false.
  !INCLUDE "ini_paramLMDZ_dyn.h"
     ! ENDIF
  !
  !INCLUDE "write_paramLMDZ_dyn.h"

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.ge.1 ) THEN
      ztit='bil phys'
      IF (planet_type.eq."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.EQ.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.eq."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.EQ. 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).EQ.0 .OR. itau.EQ.itaufin) THEN
           IF(itau.EQ.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).EQ.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.EQ.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).EQ.0 )    THEN
                GO TO 1
        ELSE IF ( MOD(itau-1,iperiod).EQ. 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.EQ. 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).EQ.0 .OR. itau.EQ.itaufin) THEN
           IF(itau.EQ.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         ).EQ.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.EQ.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