--- trunk/libf/dyn3d/leapfrog.f90	2012/10/02 15:50:56	67
+++ trunk/libf/dyn3d/leapfrog.f90	2013/06/24 15:39:52	70
@@ -19,13 +19,14 @@
     USE comgeom, ONLY: aire_2d, apoln, apols
     USE disvert_m, ONLY: ap, bp
     USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
-         iflag_phys, ok_guide
+         iflag_phys, ok_guide, iecri
     USE dimens_m, ONLY: iim, jjm, llm, nqmx
     use dissip_m, only: dissip
     USE dynetat0_m, ONLY: day_ini
     use dynredem1_m, only: dynredem1
     USE exner_hyb_m, ONLY: exner_hyb
     use filtreg_m, only: filtreg
+    use fluxstokenc_m, only: fluxstokenc
     use geopot_m, only: geopot
     USE guide_m, ONLY: guide
     use inidissip_m, only: idissip
@@ -34,6 +35,7 @@
     USE pressure_var, ONLY: p3d
     USE temps, ONLY: itau_dyn
     use writedynav_m, only: writedynav
+    use writehist_m, only: writehist
 
     ! Variables dynamiques:
     REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
@@ -43,8 +45,8 @@
     ! potential temperature
 
     REAL, intent(inout):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol, en Pa
-    REAL masse(:, :, :) ! (iim + 1, jjm + 1, llm) masse d'air
-    REAL phis(:, :) ! (iim + 1, jjm + 1) geopotentiel au sol
+    REAL, intent(inout):: masse(:, :, :) ! (iim + 1, jjm + 1, llm) masse d'air
+    REAL, intent(in):: phis(:, :) ! (iim + 1, jjm + 1) surface geopotential
 
     REAL, intent(inout):: q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nqmx)
     ! mass fractions of advected fields
@@ -55,7 +57,7 @@
 
     ! Variables dynamiques:
 
-    REAL pks((iim + 1) * (jjm + 1)) ! exner au sol
+    REAL pks(iim + 1, jjm + 1) ! exner au sol
     REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
     REAL pkf(iim + 1, jjm + 1, llm) ! exner filtré au milieu des couches
     REAL phi(iim + 1, jjm + 1, llm) ! geopotential
@@ -133,7 +135,7 @@
        end if
 
        ! Calcul des tendances dynamiques:
-       CALL geopot((iim + 1) * (jjm + 1), teta, pk, pks, phis, phi)
+       CALL geopot(teta, pk, pks, phis, phi)
        CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
             dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
             conser=MOD(itau, iconser)==0)
@@ -145,7 +147,7 @@
        IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
             dtvr, itau)
 
-       ! Integrations dynamique et traceurs:
+       ! Intégrations dynamique et traceurs:
        CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, dteta, &
             dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, finvmaold, dt, &
             leapf)
@@ -156,7 +158,7 @@
           CALL exner_hyb(ps, p3d, pks, pk, pkf)
 
           ! Calcul des tendances dynamiques:
-          CALL geopot((iim + 1) * (jjm + 1), teta, pk, pks, phis, phi)
+          CALL geopot(teta, pk, pks, phis, phi)
           CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
                phi, dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
                conser=.false.)
@@ -168,7 +170,7 @@
        end if
 
        IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
-          ! calcul des tendances physiques:
+          ! Calcul des tendances physiques:
 
           forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
           CALL exner_hyb(ps, p3d, pks, pk, pkf)
@@ -182,9 +184,9 @@
                phis, phi, dudyn, dv, dq, w, dufi, dvfi, dtetafi, dqfi, dpfi, &
                lafin = itau + 1 == itaufin)
 
-          ! ajout des tendances physiques:
-          CALL addfi(nqmx, dtphys, ucov, vcov, teta, q, ps, dufi, dvfi, &
-               dtetafi, dqfi, dpfi)
+          ! Ajout des tendances physiques:
+          CALL addfi(nqmx, ucov, vcov, teta, q, ps, dufi, dvfi, dtetafi, &
+               dqfi, dpfi)
        ENDIF
 
        forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
@@ -229,13 +231,18 @@
           call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
                q(:, :, :, 1))
        ENDIF
+
+       IF (MOD(itau + 1, iecri * day_step) == 0) THEN
+          CALL geopot(teta, pk, pks, phis, phi)
+          CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps)
+       END IF
     end do time_integration
 
     CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
          itau = itau_dyn + itaufin)
 
     ! Calcul des tendances dynamiques:
-    CALL geopot((iim + 1) * (jjm + 1), teta, pk, pks, phis, phi)
+    CALL geopot(teta, pk, pks, phis, phi)
     CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
          dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
          conser = MOD(itaufin, iconser) == 0)