!
! $Header$
!
c
c
      SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0,
     &                    time_0)


cIM : pour sortir les param. du modele dans un fis. netcdf 110106
      USE IOIPSL 

      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 "control.h"
#include "ener.h"
#include "description.h"
#include "serre.h"
#include "com_io_dyn.h"
#include "iniprint.h"
#include "advtrac.h"
c#include "tracstoke.h"

#include "academic.h"

#include "clesphys.h"

      integer nq

      INTEGER         longcles
      PARAMETER     ( longcles = 20 )
      REAL  clesphy0( longcles )

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

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
      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)

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

c   tendances physiques
      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
      REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),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*4  iday ! jour julien
      REAL       time ! Heure de la journee en fraction d'1 jour

      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
      LOGICAL first,callinigrads
cIM : pour sortir les param. du modele dans un fis. netcdf 110106
      save first
      data first/.true./
      real dt_cum, zjulian
      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)
#ifdef INCA
      REAL :: flxw(ip1jmp1,llm)
#endif

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 

      character*80 dynhist_file, dynhistave_file
      character*20 modname
      character*80 abort_message

C Calendrier
      LOGICAL true_calendar
      PARAMETER (true_calendar = .false.)

      logical dissip_conservative
      save dissip_conservative
      data dissip_conservative/.true./

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

      integer itau_w   ! pas de temps ecriture = itap + itau_phy


      itaufin   = nday*day_step
      itaufinp1 = itaufin +1


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


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

      dq=0.
      CALL pression ( ip1jmp1, ap, bp, ps, p       )
      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )

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

   1  CONTINUE


#ifdef CPP_IOIPSL
      if (ok_guide.and.(itaufin-itau-1)*dtvr.gt.21600) then
        call guide(itau,ucov,vcov,teta,q,masse,ps)
      else
        IF(prt_level>9)WRITE(*,*)'attention on ne guide pas les ',
     .    '6 dernieres heures'
      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
      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 )

      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)

   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
         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.NE.0                 ) apphys = .TRUE.
      ELSE
         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.NE.0) apphys=.TRUE.
      END IF

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

      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )

      CALL caldyn 
     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )

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

      IF( forward. OR . leapf )  THEN

c
#ifdef INCA
             CALL caladvtrac(q,pbaru,pbarv,
     *                      p, masse, dq,  teta,
     .             flxw,
     .             pk)
#else
             CALL caladvtrac(q,pbaru,pbarv,
     *                      p, masse, dq,  teta,
     .             pk)
#endif

         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
c
      ENDIF


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      )
         CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )

           rdaym_ini  = itau * dtvr / daysec
           rdayvrai   = rdaym_ini  + day_ini


c rajout debug
c       lafin = .true.


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

#ifdef CPP_PHYS
c+jld

c  Diagnostique de conservation de l'énergie : initialisation
      IF (ip_ebil_dyn.ge.1 ) THEN 
          ztit='bil dyn'
          CALL diagedyn(ztit,2,1,1,dtphys
     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
      ENDIF 
c-jld
cIM : pour sortir les param. du modele dans un fis. netcdf 110106
      IF (first) THEN
       first=.false.
#include "ini_paramLMDZ_dyn.h"
      ENDIF
c
#include "write_paramLMDZ_dyn.h"
c

        CALL calfis( nq, lafin ,rdayvrai,time  ,
     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
     $               du,dv,dteta,dq,w,
#ifdef INCA
     $               flxw,
#endif
     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )

c      ajout des tendances physiques:
c      ------------------------------
          CALL addfi( nqmx, dtphys, leapf, forward   ,
     $                  ucov, vcov, teta , q   ,ps ,
     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
c
c  Diagnostique de conservation de l'énergie : difference
      IF (ip_ebil_dyn.ge.1 ) THEN 
          ztit='bil phys'
          CALL diagedyn(ztit,2,1,1,dtphys
     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
      ENDIF 
#else

c   Calcul academique de la physique = Rappel Newtonien + fritcion 
c   --------------------------------------------------------------
       teta(:,:)=teta(:,:)
     s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
       call friction(ucov,vcov,iphysiq*dtvr)

#endif

c-jld
       ENDIF

        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )


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
        CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
        ucov=ucov+dudis
        vcov=vcov+dvdis
c       teta=teta+dtetadis


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)
            dtetaecdt= (ecin0-ecin)/ pk
c           teta=teta+dtetaecdt
            dtetadis=dtetadis+dtetaecdt
        endif
        teta=teta+dtetadis
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

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
              iday= day_ini+itau/day_step
              time= FLOAT(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  
c              write(79,*) 'ucov',ucov
c              write(80,*) 'vcov',vcov
c              write(81,*) 'teta',teta
c              write(82,*) 'ps',ps
c              write(83,*) 'q',q
c              WRITE(85,*) 'q1 = ',q(:,:,1)
c              WRITE(86,*) 'q3 = ',q(:,:,3)

              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
#ifdef CPP_IOIPSL
              CALL writedynav(histaveid, nqmx, itau,vcov ,
     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
               call bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
#endif

            ENDIF

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

            IF( MOD(itau,iecri         ).EQ.0) THEN
c           IF( MOD(itau,iecri*day_step).EQ.0) THEN

               nbetat = nbetatdem
       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
#ifdef CPP_IOIPSL
c        CALL writehist(histid,histvid, nqmx,itau,vcov, 
c     s                       ucov,teta,phi,q,masse,ps,phis)
#else
#include "write_grads_dyn.h"
#endif


            ENDIF

            IF(itau.EQ.itaufin) THEN


#ifdef CPP_IOIPSL
       CALL dynredem1("restart.nc",0.0,
     ,                     vcov,ucov,teta,q,nqmx,masse,ps)
#endif

              CLOSE(99)
            ENDIF

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
            ELSE

c      ......   pas leapfrog  .....

                 leapf = .TRUE.
                 dt  = 2.*dtvr
                 GO TO 2
            END IF

      ELSE

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

             itau =  itau + 1
             iday = day_ini+itau/day_step
             time = FLOAT(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

            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
               IF(itau.EQ.itaufin) THEN
                  iav=1
               ELSE
                  iav=0
               ENDIF
#ifdef CPP_IOIPSL
              CALL writedynav(histaveid, nqmx, itau,vcov ,
     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
               call bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
#endif

            ENDIF

               IF(MOD(itau,iecri         ).EQ.0) THEN
c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
                  nbetat = nbetatdem
       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
#ifdef CPP_IOIPSL
c       CALL writehist( histid, histvid, nqmx, itau,vcov , 
c    ,                           ucov,teta,phi,q,masse,ps,phis)
#else
#include "write_grads_dyn.h"
#endif


               ENDIF

#ifdef CPP_IOIPSL
                 IF(itau.EQ.itaufin)
     . CALL dynredem1("restart.nc",0.0,
     .                     vcov,ucov,teta,q,nqmx,masse,ps)
#endif

                 forward = .TRUE.
                 GO TO  1

            ENDIF

      END IF

      STOP
      END
