!
! $Id: leapfrog.F 1446 2010-10-22 09:27:25Z emillour $
!
c
c
      SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,
     &                    time_0)


cIM : pour sortir les param. du modele dans un fis. netcdf 110106
#ifdef CPP_IOIPSL
      use IOIPSL
#endif
      USE infotrac
      USE guide_mod, ONLY : guide_main
      USE write_field
      USE control_mod
      IMPLICIT NONE

c      ......   Version  du 10/01/98    ..........

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

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

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

#include "dimensions.h"
#include "paramet.h"
#include "comconst.h"
#include "comdissnew.h"
#include "comvert.h"
#include "comgeom.h"
#include "logic.h"
#include "temps.h"
#include "ener.h"
#include "description.h"
#include "serre.h"
!#include "com_io_dyn.h"
#include "iniprint.h"
#include "academic.h"

! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
! #include "clesphys.h"

      real zqmin,zqmax
      INTEGER nbetatmoy, nbetatdem,nbetat

c   variables dynamiques
      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
      REAL teta(ip1jmp1,llm)                 ! temperature potentielle 
      REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
      REAL ps(ip1jmp1)                       ! pression  au sol
      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
      REAL pks(ip1jmp1)                      ! exner au  sol
      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
      REAL masse(ip1jmp1,llm)                ! masse d'air
      REAL phis(ip1jmp1)                     ! geopotentiel au sol
      REAL phi(ip1jmp1,llm)                  ! geopotentiel
      REAL w(ip1jmp1,llm)                    ! vitesse verticale
! ADAPTATION GCM POUR CP(T)
      REAL temp(ip1jmp1,llm)                 ! temperature  
      REAL tsurpk(ip1jmp1,llm)               ! cpp*T/pk  

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

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

c   tendances dynamiques en */s
      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)

c   tendances de la dissipation en */s
      REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
      REAL dtetadis(ip1jmp1,llm)

c   tendances de la couche superieure */s
      REAL dvtop(ip1jm,llm),dutop(ip1jmp1,llm)
      REAL dtetatop(ip1jmp1,llm)
      REAL dqtop(ip1jmp1,llm,nqtot),dptop(ip1jmp1)

c   tendances physiques */s
      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
      REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)

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

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

      REAL  SSUM
      REAL time_0 , finvmaold(ip1jmp1,llm)

cym      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
cIM : pour sortir les param. du modele dans un fis. netcdf 110106
      save first
      data first/.true./
      real dt_cum
      character*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.) 

      data callinigrads/.true./
      character*10 string10

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

c+jld variables test conservation energie
      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
C     Tendance de la temp. potentiel d (theta)/ d t due a la 
C     tansformation d'energie cinetique en energie thermique
C     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*15 ztit
!IM   INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
!IM   SAVE      ip_ebil_dyn
!IM   DATA      ip_ebil_dyn/0/
c-jld 

      integer :: itau_w ! for write_paramLMDZ_dyn.h

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

      logical dissip_conservative
      save dissip_conservative
      data dissip_conservative/.true./

      INTEGER testita
      PARAMETER (testita = 9)

      logical , parameter :: flag_verif = .false.
      
      ! for CP(T)
      real :: dtec
      real,external :: cpdet
      real :: ztetaec(ip1jmp1,llm)

c dummy: sinon cette routine n'est jamais compilee...
      if(1.eq.0) then
        CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
      endif

      itaufin   = nday*day_step
      if (less1day) then
c MODIF VENUS: to run less than one day:
        itaufin   = int(fractday*day_step)
      endif
      itaufinp1 = itaufin +1
      
c INITIALISATIONS
        dudis(:,:)   =0.
        dvdis(:,:)   =0.
        dtetadis(:,:)=0.
        dutop(:,:)   =0.
        dvtop(:,:)   =0.
        dtetatop(:,:)=0.
        dqtop(:,:,:) =0.
        dptop(:)     =0.
        dufi(:,:)   =0.
        dvfi(:,:)   =0.
        dtetafi(:,:)=0.
        dqfi(:,:,:) =0.
        dpfi(:)     =0.

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


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

      dq(:,:,:)=0.
      CALL pression ( ip1jmp1, ap, bp, ps, p       )
      if (disvert_type==1) then
        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
      else ! we assume that we are in the disvert_type==2 case
        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
      endif

c------------------
c TEST PK MONOTONE
c------------------
      write(*,*) "Test PK"
      do ij=1,ip1jmp1
        do l=2,llm
          if(pk(ij,l).gt.pk(ij,l-1)) then
c           write(*,*) ij,l,pk(ij,l)
            abort_message = 'PK non strictement decroissante'
            call abort_gcm(modname,abort_message,1)
c           write(*,*) "ATTENTION, Test PK deconnecté..."
          endif
        enddo
      enddo
      write(*,*) "Fin Test PK"
c     stop
c------------------

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

   1  CONTINUE

      jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec) 
      jH_cur = jH_ref +                                                 &
     &          (itau * dtvr / daysec - int(itau * dtvr / daysec)) 


#ifdef CPP_IOIPSL
      if (ok_guide) then
        call guide_main(itau,ucov,vcov,teta,q,masse,ps)
      endif
#endif


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

! 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

c   ...    P.Le Van .26/04/94  ....

      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )

   2  CONTINUE

c-----------------------------------------------------------------------

c   date:
c   -----


c   gestion des appels de la physique et des dissipations:
c   ------------------------------------------------------
c
c   ...    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,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward 
     s          .and. iflag_phys.EQ.1                 ) apphys = .TRUE.
      ELSE
      ! Leapfrog/Matsuno time stepping 
         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
      END IF

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

c-----------------------------------------------------------------------
c   calcul des tendances dynamiques:
c   --------------------------------

! ADAPTATION GCM POUR CP(T)
      call tpot2t(ijp1llm,teta,temp,pk)
      tsurpk = cpp*temp/pk
      CALL geopot  ( ip1jmp1, tsurpk  , pk , pks,  phis  , phi   )

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


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

      IF( forward. OR . leapf )  THEN

         CALL caladvtrac(q,pbaru,pbarv,
     *        p, masse, dq,  teta,
     .        flxw, pk)
         
         IF (offline) THEN
Cmaf stokage du flux de masse pour traceurs OFF-LINE

#ifdef CPP_IOIPSL
           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
     .   dtvr, itau)
#endif


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


c-----------------------------------------------------------------------
c   integrations dynamique et traceurs:
c   ----------------------------------


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


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

         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
         if (disvert_type==1) then
           CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
         else ! we assume that we are in the disvert_type==2 case
           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
         endif

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

c rajout debug
c       lafin = .true.


c   Interface avec les routines de phylmd (phymars ... )
c   -----------------------------------------------------

c+jld

c  Diagnostique de conservation de l'énergie : initialisation
         IF (ip_ebil_dyn.ge.1 ) THEN 
          ztit='bil dyn'
! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
           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 )
c-jld
#ifdef CPP_IOIPSL
cIM : pour sortir les param. du modele dans un fis. netcdf 110106
         IF (first) THEN
#include "ini_paramLMDZ_dyn.h"
          first=.false.
         ENDIF
c
#include "write_paramLMDZ_dyn.h"
c
#endif
! #endif of #ifdef CPP_IOIPSL

         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  )

c      ajout des tendances physiques:
c      ------------------------------
          CALL addfi( dtphys, leapf, forward   ,
     $                  ucov, vcov, teta , q   ,ps ,
     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )

c      Couche superieure :
c      -------------------
         IF (ok_strato) THEN
           CALL top_bound( vcov,ucov,teta,phi,masse,
     $                     dutop,dvtop,dtetatop)
c dqtop=0, dptop=0
           CALL addfi( dtphys, leapf, forward   ,
     $                  ucov, vcov, teta , q   ,ps ,
     $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
         ENDIF

c  Diagnostique de conservation de l'énergie : 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 
          call friction(ucov,vcov,dtvr)
    
       if (planet_type.eq."giant") then
          ! Intrinsic heat flux
          ! Aymeric -- for giant planets
          if (ihf .gt. 1.e-6) then
          !print *, '**** INTRINSIC HEAT FLUX ****', ihf 
            DO ij=1,ip1jmp1
              teta(ij,1) = teta(ij,1)
     &        + dtvr * aire(ij) * ihf / cpp / masse(ij,1)
            ENDDO
          !print *, '**** d teta ' 
          !print *, dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
          endif
       endif

    
        ! Sponge layer (if any)
        IF (ok_strato) THEN
          CALL top_bound(vcov,ucov,teta,phi,
     $                   masse,dutop,dvtop,dtetatop)
c dqtop=0, dptop=0
          CALL addfi( dtvr, leapf, forward   ,
     $                  ucov, vcov, teta , q   ,ps ,
     $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
        ENDIF ! of IF (ok_strato) 
      ENDIF ! of IF (iflag_phys.EQ.2)


c-jld

        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
        if (disvert_type==1) then
          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
        else ! we assume that we are in the disvert_type==2 case
          CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
        endif


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

      IF(apdiss) THEN


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

c   dissipation
! ADAPTATION GCM POUR CP(T)
        call tpot2t(ijp1llm,teta,temp,pk)

        CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
        ucov=ucov+dudis
        vcov=vcov+dvdis
        dudis=dudis/dtdiss   ! passage en (m/s)/s
        dvdis=dvdis/dtdiss   ! passage en (m/s)/s

c------------------------------------------------------------------------
        if (dissip_conservative) then
C       On rajoute la tendance due a la transform. Ec -> E therm. cree
C       lors de la dissipation
            call covcont(llm,ucov,vcov,ucont,vcont)
            call enercin(vcov,ucov,vcont,ucont,ecin)
! ADAPTATION GCM POUR CP(T)
            do ij=1,ip1jmp1
              do l=1,llm
                dtec = (ecin0(ij,l)-ecin(ij,l))/cpdet(temp(ij,l))
                temp(ij,l) = temp(ij,l) + dtec
              enddo
            enddo
            call t2tpot(ijp1llm,temp,ztetaec,pk)
            dtetaecdt=ztetaec-teta
            dtetadis=dtetadis+dtetaecdt
        endif
        teta=teta+dtetadis
        dtetadis=dtetadis/dtdiss   ! passage en K/s
c------------------------------------------------------------------------


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

        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

        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


      END IF ! of IF(apdiss)

c ajout debug
c              IF( lafin ) then  
c                abort_message = 'Simulation finished'
c                call abort_gcm(modname,abort_message,0)
c              ENDIF
        
c   ********************************************************************
c   ********************************************************************
c   .... fin de l'integration dynamique  et physique pour le pas itau ..
c   ********************************************************************
c   ********************************************************************

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

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

            IF(forward. OR. leapf) THEN
              itau= itau + 1
c              iday= day_ini+itau/day_step
c              time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
c                IF(time.GT.1.) THEN
c                  time = time-1.
c                  iday = iday+1
c                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
c-----------------------------------------------------------------------
c   ecriture du fichier histoire moyenne:
c   -------------------------------------

            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
               IF(itau.EQ.itaufin) THEN
                  iav=1
               ELSE
                  iav=0
               ENDIF
               
               IF (ok_dynzon) THEN
#ifdef CPP_IOIPSL
c les traceurs ne sont pas sortis, trop lourd. 
c Peut changer eventuellement si besoin.
                 CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav,
     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
     &                 du,dudis,dutop,dufi)
#endif
               END IF
               IF (ok_dyn_ave) THEN
#ifdef CPP_IOIPSL
                 CALL writedynav(itau,vcov,
     &                 ucov,teta,pk,phi,q,masse,ps,phis)
#endif
               ENDIF

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

c-----------------------------------------------------------------------
c   ecriture de la bande histoire:
c   ------------------------------

            IF( MOD(itau,iecri).EQ.0) THEN
             ! Ehouarn: output only during LF or Backward Matsuno
	     if (leapf.or.(.not.leapf.and.(.not.forward))) then
              nbetat = nbetatdem
! ADAPTATION GCM POUR CP(T)
              call tpot2t(ijp1llm,teta,temp,pk)
              tsurpk = cpp*temp/pk
              CALL geopot(ip1jmp1,tsurpk,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
#ifdef CPP_IOIPSL
              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)
#endif
! 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."mars") then
! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem1_mars
                abort_message = 'dynredem1_mars A FAIRE'
                call abort_gcm(modname,abort_message,0)
              else
                CALL dynredem1("restart.nc",0.0,
     &                         vcov,ucov,teta,q,masse,ps)
              endif ! of if (planet_type.eq."mars")

              CLOSE(99)
            ENDIF ! of IF (itau.EQ.itaufin)

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

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

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

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

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

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

c      ......   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)

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

             itau =  itau + 1
c             iday = day_ini+itau/day_step
c             time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
c
c                  IF(time.GT.1.) THEN
c                   time = time-1.
c                   iday = iday+1
c                  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

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

               IF (ok_dynzon) THEN 
#ifdef CPP_IOIPSL
c les traceurs ne sont pas sortis, trop lourd. 
c Peut changer eventuellement si besoin.
                 CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav,
     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
     &                 du,dudis,dutop,dufi)
#endif
               ENDIF
               IF (ok_dyn_ave) THEN
#ifdef CPP_IOIPSL
                 CALL writedynav(itau,vcov,
     &                 ucov,teta,pk,phi,q,masse,ps,phis)
#endif
               ENDIF

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

              IF(MOD(itau,iecri         ).EQ.0) THEN
c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
                nbetat = nbetatdem
! ADAPTATION GCM POUR CP(T)
                call tpot2t(ijp1llm,teta,temp,pk)
                tsurpk = cpp*temp/pk
                CALL geopot(ip1jmp1,tsurpk,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
#ifdef CPP_IOIPSL
              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)
#endif
! 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."mars") then
! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem1_mars
                  abort_message = 'dynredem1_mars A FAIRE'
                  call abort_gcm(modname,abort_message,0)
                else
                  CALL dynredem1("restart.nc",0.0,
     &                         vcov,ucov,teta,q,masse,ps)
                endif ! of if (planet_type.eq."mars")
              ENDIF ! of IF(itau.EQ.itaufin)

              forward = .TRUE.
              GO TO  1

            ENDIF ! of IF (forward)

      END IF ! of IF(.not.purmats)

      STOP
      END