Index: LMDZ4/trunk/libf/dyn3d/leapfrog.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/leapfrog.F	(revision 1060)
+++ LMDZ4/trunk/libf/dyn3d/leapfrog.F	(revision 1146)
@@ -2,10 +2,13 @@
 c
 c
-      SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0,
+      SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
      &                    time_0)
 
 
 cIM : pour sortir les param. du modele dans un fis. netcdf 110106
-      USE IOIPSL 
+#ifdef CPP_IOIPSL
+      use IOIPSL
+#endif
+      USE infotrac
 
       IMPLICIT NONE
@@ -56,13 +59,8 @@
 #include "com_io_dyn.h"
 #include "iniprint.h"
-#include "advtrac.h"
-c#include "tracstoke.h"
-
 #include "academic.h"
 
 ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
 ! #include "clesphys.h"
-
-      integer nq
 
       INTEGER         longcles
@@ -76,5 +74,5 @@
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
       REAL teta(ip1jmp1,llm)                 ! temperature potentielle 
-      REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
+      REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
       REAL ps(ip1jmp1)                       ! pression  au sol
       REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
@@ -97,5 +95,5 @@
 c   tendances dynamiques
       REAL dv(ip1jm,llm),du(ip1jmp1,llm)
-      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
+      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
 
 c   tendances de la dissipation
@@ -105,5 +103,5 @@
 c   tendances physiques
       REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
-      REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
+      REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
 
 c   variables pour le fichier histoire
@@ -165,5 +163,5 @@
 
       character*80 dynhist_file, dynhistave_file
-      character*20 modname
+      character(len=20) :: modname
       character*80 abort_message
 
@@ -182,5 +180,5 @@
       PARAMETER (testita = 9)
 
-      logical , parameter :: flag_verif = .false.
+      logical , parameter :: flag_verif = .true.
       
 
@@ -190,5 +188,6 @@
       itaufin   = nday*day_step
       itaufinp1 = itaufin +1
-
+      modname="leapfrog"
+      
 
       itau = 0
@@ -220,6 +219,6 @@
         call guide(itau,ucov,vcov,teta,q,masse,ps)
       else
-        IF(prt_level>9)WRITE(*,*)'attention on ne guide pas les ',
-     .    '6 dernieres heures'
+        IF(prt_level>9)WRITE(lunout,*)'leapfrog: attention on ne ',
+     .    'guide pas les 6 dernieres heures'
       endif
 #endif
@@ -230,4 +229,6 @@
 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 )
@@ -245,5 +246,6 @@
       CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
 
-      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
+! Ehouarn: what is this for? zqmin & zqmax are not used anyway ...
+!      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
 
    2  CONTINUE
@@ -305,7 +307,7 @@
 
 
-         ENDIF
-c
-      ENDIF
+         ENDIF ! of IF (offline)
+c
+      ENDIF ! of IF( forward. OR . leapf )
 
 
@@ -353,24 +355,28 @@
 c   -----------------------------------------------------
 
-#ifdef CPP_PHYS
 c+jld
 
 c  Diagnostique de conservation de l'énergie : initialisation
-      IF (ip_ebil_dyn.ge.1 ) THEN 
+         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 
+! 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
-       first=.false.
+         IF (first) THEN
+          first=.false.
 #include "ini_paramLMDZ_dyn.h"
-      ENDIF
+         ENDIF
 c
 #include "write_paramLMDZ_dyn.h"
 c
-
-        CALL calfis( nq, lafin ,rdayvrai,time  ,
+#endif
+! #endif of #ifdef CPP_IOIPSL
+         CALL calfis( lafin ,rdayvrai,time  ,
      $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
      $               du,dv,dteta,dq,
@@ -378,24 +384,26 @@
      $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
 
-       IF (ok_strato) THEN
-         CALL top_bound( vcov,ucov,teta, dufi,dvfi,dtetafi)
-       ENDIF
+         IF (ok_strato) THEN
+           CALL top_bound( vcov,ucov,teta, dufi,dvfi,dtetafi)
+         ENDIF
        
 c      ajout des tendances physiques:
 c      ------------------------------
-          CALL addfi( nqmx, dtphys, leapf, forward   ,
+          CALL addfi( 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 
+         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 
-#endif
+          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 physics" case
+      IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
 c   Calcul academique de la physique = Rappel Newtonien + friction 
 c   --------------------------------------------------------------
@@ -475,5 +483,5 @@
 
 
-      END IF
+      END IF ! of IF(apdiss)
 
 c ajout debug
@@ -509,11 +517,16 @@
             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(80,*) 'ucov',ucov
+                write(81,*) 'vcov',vcov
+                write(82,*) 'teta',teta
+                write(83,*) 'ps',ps
+                write(84,*) 'q',q
                 WRITE(85,*) 'q1 = ',q(:,:,1)
                 WRITE(86,*) 'q3 = ',q(:,:,3)
+                write(90) ucov
+                write(91) vcov
+                write(92) teta
+                write(93) ps
+                write(94) q
               endif
 
@@ -532,10 +545,13 @@
                   iav=0
                ENDIF
+               
+               IF (ok_dynzon) THEN
 #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
+                  CALL writedynav(histaveid, 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
+               END IF
 
             ENDIF
@@ -548,31 +564,35 @@
 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
+              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
+c             CALL writehist(histid,histvid,itau,vcov, 
+c     &                      ucov,teta,phi,q,masse,ps,phis)
+#endif
+! For some Grads outputs of fields
+             if (output_grads_dyn) then
 #include "write_grads_dyn.h"
-#endif
-
-
-            ENDIF
+             endif
+
+            ENDIF ! of IF(MOD(itau,iecri).EQ.0)
 
             IF(itau.EQ.itaufin) THEN
 
 
-#ifdef CPP_IOIPSL
-       CALL dynredem1("restart.nc",0.0,
-     ,                     vcov,ucov,teta,q,nqmx,masse,ps)
-#endif
+              if (planet_type.eq."earth") then
+#ifdef CPP_EARTH
+! Write an Earth-format restart file
+                CALL dynredem1("restart.nc",0.0,
+     &                         vcov,ucov,teta,q,masse,ps)
+#endif
+              endif ! of if (planet_type.eq."earth")
 
               CLOSE(99)
-            ENDIF
+            ENDIF ! of IF (itau.EQ.itaufin)
 
 c-----------------------------------------------------------------------
@@ -596,6 +616,6 @@
                         leapf =  .TRUE.
                         dt  =  2.*dtvr
-                        GO TO 2
-                   END IF
+                        GO TO 2 
+                   END IF ! of IF (forward)
             ELSE
 
@@ -605,7 +625,8 @@
                  dt  = 2.*dtvr
                  GO TO 2
-            END IF
-
-      ELSE
+            END IF ! of IF (MOD(itau,iperiod).EQ.0)
+                   !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
+
+      ELSE ! of IF (.not.purmats)
 
 c       ........................................................
@@ -630,7 +651,7 @@
                GO TO 2
 
-            ELSE
-
-            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
+            ELSE ! of IF(forward)
+
+              IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
                IF(itau.EQ.itaufin) THEN
                   iav=1
@@ -638,44 +659,51 @@
                   iav=0
                ENDIF
+
+               IF (ok_dynzon) THEN 
 #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
+                  CALL writedynav(histaveid, 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
+               END IF
+
+              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
-       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
+                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
+c               CALL writehist( histid, histvid, itau,vcov , 
+c    &                           ucov,teta,phi,q,masse,ps,phis)
+#endif
+! For some Grads outputs
+                if (output_grads_dyn) then
 #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
+                endif
+
+              ENDIF ! of IF(MOD(itau,iecri         ).EQ.0) 
+
+              IF(itau.EQ.itaufin) THEN
+                if (planet_type.eq."earth") then
+#ifdef CPP_EARTH
+                  CALL dynredem1("restart.nc",0.0,
+     &                           vcov,ucov,teta,q,masse,ps)
+#endif
+                endif ! of if (planet_type.eq."earth")
+              ENDIF ! of IF(itau.EQ.itaufin)
+
+              forward = .TRUE.
+              GO TO  1
+
+            ENDIF ! of IF (forward)
+
+      END IF ! of IF(.not.purmats)
 
       STOP
