Index: LMDZ6/trunk/libf/dyn3d/leapfrog.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3d/leapfrog.F	(revision 4143)
+++ LMDZ6/trunk/libf/dyn3d/leapfrog.F90	(revision 5246)
@@ -2,867 +2,867 @@
 ! $Id$
 !
-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, 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
-
-      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 "comdissnew.h"
-      include "comgeom.h"
-      include "description.h"
-      include "iniprint.h"
-      include "academic.h"
-
-      REAL,INTENT(IN) :: time_0 ! not used
-
-c   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
-
-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,nqtot),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,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 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.) 
-      logical physic
-
-      data callinigrads/.true./
-      character*10 string10
-
-      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 
-
-      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./
-
-      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
+!
+!
+SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
+
+
+  !IM : pour sortir les param. du modele dans un fis. netcdf 110106
+#ifdef CPP_IOIPSL
+  use IOIPSL
+#endif
+  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
+
+  IMPLICIT NONE
+
+   ! ......   Version  du 10/01/98    ..........
+
+   !        avec  coordonnees  verticales hybrides
+  !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
+
+  !=======================================================================
+  !
+  !   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
+  !   -------
+  !
+  !   Objet:
+  !   ------
+  !
+  !   GCM LMD nouvelle grille
+  !
+  !=======================================================================
+  !
+  !  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
+  !  et possibilite d'appeler une fonction f(y)  a derivee tangente
+  !  hyperbolique a la  place de la fonction a derivee sinusoidale.
+
+  !  ... Possibilite de choisir le shema pour l'advection de
+  !    q  , en modifiant iadv dans traceur.def  (10/02) .
+  !
+  !  Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
+  !  Pour Van-Leer iadv=10
+  !
+  !-----------------------------------------------------------------------
+  !   Declarations:
+  !   -------------
+
+  include "dimensions.h"
+  include "paramet.h"
+  include "comdissnew.h"
+  include "comgeom.h"
+  include "description.h"
+  include "iniprint.h"
+  include "academic.h"
+
+  REAL,INTENT(IN) :: time_0 ! not used
+
+  !   dynamical variables:
+  REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm)    ! zonal covariant wind
+  REAL,INTENT(INOUT) :: vcov(ip1jm,llm)      ! meridional covariant wind
+  REAL,INTENT(INOUT) :: teta(ip1jmp1,llm)    ! potential temperature
+  REAL,INTENT(INOUT) :: ps(ip1jmp1)          ! surface pressure (Pa)
+  REAL,INTENT(INOUT) :: masse(ip1jmp1,llm)   ! air mass
+  REAL,INTENT(INOUT) :: phis(ip1jmp1)        ! geopotentiat at the surface
+  REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
+
+  REAL :: p (ip1jmp1,llmp1  )               ! interlayer pressure
+  REAL :: pks(ip1jmp1)                      ! exner at the surface
+  REAL :: pk(ip1jmp1,llm)                   ! exner at mid-layer
+  REAL :: pkf(ip1jmp1,llm)                  ! filtered exner at mid-layer
+  REAL :: phi(ip1jmp1,llm)                  ! geopotential
+  REAL :: w(ip1jmp1,llm)                    ! vertical velocity
+
+  real :: zqmin,zqmax
+
+  ! variables dynamiques intermediaire pour le transport
+  REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
+
+  !   variables dynamiques au pas -1
+  REAL :: vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
+  REAL :: tetam1(ip1jmp1,llm),psm1(ip1jmp1)
+  REAL :: massem1(ip1jmp1,llm)
+
+  !   tendances dynamiques
+  REAL :: dv(ip1jm,llm),du(ip1jmp1,llm)
+  REAL :: dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
+
+  !   tendances de la dissipation
+  REAL :: dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
+  REAL :: dtetadis(ip1jmp1,llm)
+
+  !   tendances physiques
+  REAL :: dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
+  REAL :: dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
+
+  !   variables pour le fichier histoire
+  REAL :: dtav      ! intervalle de temps elementaire
+
+  REAL :: tppn(iim),tpps(iim),tpn,tps
+  !
+  INTEGER :: itau,itaufinp1,iav
+   ! INTEGER  iday ! jour julien
+  REAL :: time
+
+  REAL :: SSUM
+  ! REAL finvmaold(ip1jmp1,llm)
+
+  !ym      LOGICAL  lafin
+  LOGICAL :: lafin=.false.
+  INTEGER :: ij,iq,l
+  INTEGER :: ik
+
+  real :: time_step, t_wrt, t_ops
+
+   ! REAL rdayvrai,rdaym_ini
+  ! jD_cur: jour julien courant
+  ! jH_cur: heure julienne courante
+  REAL :: jD_cur, jH_cur
+  INTEGER :: an, mois, jour
+  REAL :: secondes
+
+  LOGICAL :: first,callinigrads
+  !IM : pour sortir les param. du modele dans un fis. netcdf 110106
+  save first
+  data first/.true./
+  real :: dt_cum
+  character(len=10) :: infile
+  integer :: zan, tau0, thoriid
+  integer :: nid_ctesGCM
+  save nid_ctesGCM
+  real :: degres
+  real :: rlong(iip1), rlatg(jjp1)
+  real :: zx_tmp_2d(iip1,jjp1)
+  integer :: ndex2d(iip1*jjp1)
+  logical :: ok_sync
+  parameter (ok_sync = .true.)
+  logical :: physic
+
+  data callinigrads/.true./
+  character(len=10) :: string10
+
+  REAL :: flxw(ip1jmp1,llm)  ! flux de masse verticale
+
+  !+jld variables test conservation energie
+  REAL :: ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
+  ! Tendance de la temp. potentiel d (theta)/ d t due a la
+  ! tansformation d'energie cinetique en energie thermique
+  ! cree par la dissipation
+  REAL :: dtetaecdt(ip1jmp1,llm)
+  REAL :: vcont(ip1jm,llm),ucont(ip1jmp1,llm)
+  REAL :: vnat(ip1jm,llm),unat(ip1jmp1,llm)
+  REAL :: d_h_vcol, d_qt, d_qw, d_ql, d_ec
+  CHARACTER(len=15) :: ztit
+  !IM   INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
+  !IM   SAVE      ip_ebil_dyn
+  !IM   DATA      ip_ebil_dyn/0/
+  !-jld
+
+  character(len=80) :: dynhist_file, dynhistave_file
+  character(len=*),parameter :: modname="leapfrog"
+  character(len=80) :: abort_message
+
+  logical :: dissip_conservative
+  save dissip_conservative
+  data dissip_conservative/.true./
+
+  LOGICAL :: prem
+  save prem
+  DATA prem/.true./
+  INTEGER :: testita
+  PARAMETER (testita = 9)
+
+  logical , parameter :: flag_verif = .false.
+
+
+  integer :: itau_w   ! pas de temps ecriture = itap + itau_phy
+
+
+  if (nday>=0) then
+     itaufin   = nday*day_step
+  else
+     itaufin   = -nday
+  endif
+  itaufinp1 = itaufin +1
+  itau = 0
+  physic=.true.
+  if (iflag_phys==0.or.iflag_phys==2) physic=.false.
+
+   ! iday = day_ini+itau/day_step
+   ! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
+   !    IF(time.GT.1.) THEN
+   !     time = time-1.
+   !     iday = iday+1
+   !    ENDIF
+
+
+  !-----------------------------------------------------------------------
+  !   On initialise la pression et la fonction d'Exner :
+  !   --------------------------------------------------
+
+  dq(:,:,:)=0.
+  CALL pression ( ip1jmp1, ap, bp, ps, p       )
+  if (pressure_exner) then
+    CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
+  else
+    CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
+  endif
+
+  !-----------------------------------------------------------------------
+  !   Debut de l'integration temporelle:
+  !   ----------------------------------
+
+   1   CONTINUE ! Matsuno Forward step begins here
+
+  !   date: (NB: date remains unchanged for Backward step)
+  !   -----
+
+  jD_cur = jD_ref + day_ini - day_ref +                             &
+        (itau+1)/day_step
+  jH_cur = jH_ref + start_time +                                    &
+        mod(itau+1,day_step)/float(day_step)
+  jD_cur = jD_cur + int(jH_cur)
+  jH_cur = jH_cur - int(jH_cur)
+
+  call check_isotopes_seq(q,ip1jmp1,'leapfrog 321')
+
+#ifdef CPP_IOIPSL
+  if (ok_guide) then
+    call guide_main(itau,ucov,vcov,teta,q,masse,ps)
+  endif
+#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
+
+#ifdef CPP_IOIPSL
+       CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis, &
+             dtvr, itau)
+#endif
+
+
+     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
+#ifdef CPP_IOIPSL
+  !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"
+  !
+#endif
+  ! #endif of #ifdef CPP_IOIPSL
+#ifdef CPP_PHYS
+     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  )
+#endif
+   ! 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
-         itaufin   = -nday
+        CALL exner_milieu(ip1jmp1,ps,p,pks,pk,pkf)
       endif
-      itaufinp1 = itaufin +1
-      itau = 0
-      physic=.true.
-      if (iflag_phys==0.or.iflag_phys==2) physic=.false.
-
-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 (pressure_exner) then
-        CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
-      else
-        CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
-      endif
-
-c-----------------------------------------------------------------------
-c   Debut de l'integration temporelle:
-c   ----------------------------------
-
-   1  CONTINUE ! Matsuno Forward step begins here
-
-c   date: (NB: date remains unchanged for Backward step)
-c   -----
-
-      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')
-
-#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  ....
-! 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
-
-c-----------------------------------------------------------------------
-
-c   date: (NB: only leapfrog step requires recomputing date)
-c   -----
-
-      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)
+
+     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
-
-
-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,dissip_period ).EQ.0.AND..NOT.forward ) 
-     s        apdiss = .TRUE.
-         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward 
-     s          .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 )
-     s        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')
-
-c-----------------------------------------------------------------------
-c   calcul des tendances dynamiques:
-c   --------------------------------
-
-      ! 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 )
-
-
-c-----------------------------------------------------------------------
-c   calcul des tendances advection des traceurs (dont l'humidite)
-c   -------------------------------------------------------------
-
-      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
-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 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')
-
-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 (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 ! 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
 
-           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 
-
-c rajout debug
-c       lafin = .true.
-
-
-c   Inbterface avec les routines de phylmd (phymars ... )
-c   -----------------------------------------------------
-
-c+jld
-
-c  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))
+           ! ! Ehouarn: re-compute geopotential for outputs
+           CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
+
+           IF (ok_dynzon) THEN
+#ifdef CPP_IOIPSL
+             CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav, &
+                   ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
+#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 (ip_ebil_dyn.ge.1 )
-c-jld
-#ifdef CPP_IOIPSL
-cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ 
-cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
-c        IF (first) THEN
-c         first=.false.
-c#include "ini_paramLMDZ_dyn.h"
-c        ENDIF
-c
-c#include "write_paramLMDZ_dyn.h"
-c
-#endif
-! #endif of #ifdef CPP_IOIPSL
-#ifdef CPP_PHYS
-         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  )
-#endif
-c      ajout des tendances physiques:
-c      ------------------------------
-          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 ! 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
+#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
-
-         IF (ok_strato) THEN
-           CALL top_bound( vcov,ucov,teta,masse,dtphys)
-         ENDIF
-       
-c
-c  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)
-
-
-c-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')
-
-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
-
-        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)
-
-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  ......
-
-      call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
-
-      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
+         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
-                  iav=0
-               ENDIF
-               
-!              ! Ehouarn: re-compute geopotential for outputs
-               CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
-
-               IF (ok_dynzon) THEN
-#ifdef CPP_IOIPSL
-                 CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
-     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
-#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))
-
-            call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
-
-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
-              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
-              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
+   ! 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
+#ifdef CPP_IOIPSL
+             CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav, &
+                   ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
+#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
+           ! 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
+#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 (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)
-
-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)
-
-            call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
-
-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
- 
-              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 
-#ifdef CPP_IOIPSL
-                 CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
-     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
-#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
-                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
-              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."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
+            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
