Sources/dyn3d/leapfrog.f

module leapfrog_m

  IMPLICIT NONE

contains

  SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)

    ! From dyn3d/leapfrog.F, version 1.6, 2005/04/13 08:58:34 revision 616
    ! Authors: P. Le Van, L. Fairhead, F. Hourdin
    ! Matsuno-leapfrog scheme.

    use addfi_m, only: addfi
    use bilan_dyn_m, only: bilan_dyn
    use caladvtrac_m, only: caladvtrac
    use caldyn_m, only: caldyn
    USE calfis_m, ONLY: calfis
    USE comconst, ONLY: daysec, dtphys, dtvr
    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, 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
    use integrd_m, only: integrd
    use nr_util, only: assert
    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
    REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant

    REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
    ! potential temperature

    REAL, intent(inout):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol, en Pa
    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

    REAL, intent(in):: time_0

    ! Variables local to the procedure:

    ! Variables dynamiques:

    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
    REAL w((iim + 1) * (jjm + 1), llm) ! vitesse verticale

    ! Variables dynamiques intermediaire pour le transport
    ! Flux de masse :
    REAL pbaru((iim + 1) * (jjm + 1), llm), pbarv((iim + 1) * jjm, llm)

    ! Variables dynamiques au pas - 1
    REAL vcovm1(iim + 1, jjm, llm), ucovm1(iim + 1, jjm + 1, llm)
    REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
    REAL massem1(iim + 1, jjm + 1, llm)

    ! Tendances dynamiques
    REAL dv((iim + 1) * jjm, llm), dudyn((iim + 1) * (jjm + 1), llm)
    REAL dteta(iim + 1, jjm + 1, llm), dq((iim + 1) * (jjm + 1), llm, nqmx)
    real dp((iim + 1) * (jjm + 1))

    ! Tendances de la dissipation :
    REAL dvdis(iim + 1, jjm, llm), dudis(iim + 1, jjm + 1, llm)
    REAL dtetadis(iim + 1, jjm + 1, llm)

    ! Tendances physiques
    REAL dvfi((iim + 1) * jjm, llm), dufi((iim + 1) * (jjm + 1), llm)
    REAL dtetafi(iim + 1, jjm + 1, llm), dqfi((iim + 1) * (jjm + 1), llm, nqmx)
    real dpfi((iim + 1) * (jjm + 1))

    ! Variables pour le fichier histoire

    INTEGER itau ! index of the time step of the dynamics, starts at 0
    INTEGER itaufin
    REAL time ! time of day, as a fraction of day length
    real finvmaold(iim + 1, jjm + 1, llm)
    INTEGER l
    REAL rdayvrai, rdaym_ini

    ! Variables test conservation energie
    REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)

    REAL vcont((iim + 1) * jjm, llm), ucont((iim + 1) * (jjm + 1), llm)
    logical leapf
    real dt

    !---------------------------------------------------

    print *, "Call sequence information: leapfrog"
    call assert(shape(ucov) == (/iim + 1, jjm + 1, llm/), "leapfrog")

    itaufin = nday * day_step
    ! "day_step" is a multiple of "iperiod", therefore so is "itaufin".

    dq = 0.

    ! On initialise la pression et la fonction d'Exner :
    forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
    CALL exner_hyb(ps, p3d, pks, pk, pkf)

    time_integration: do itau = 0, itaufin - 1
       leapf = mod(itau, iperiod) /= 0
       if (leapf) then
          dt = 2 * dtvr
       else
          ! Matsuno
          dt = dtvr
          if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
               call guide(itau, ucov, vcov, teta, q, masse, ps)
          vcovm1 = vcov
          ucovm1 = ucov
          tetam1 = teta
          massem1 = masse
          psm1 = ps
          finvmaold = masse
          CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE.)
       end if

       ! Calcul des tendances dynamiques:
       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)

       ! Calcul des tendances advection des traceurs (dont l'humidité)
       CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)

       ! Stokage du flux de masse pour traceurs offline:
       IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
            dtvr, itau)

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

       if (.not. leapf) then
          ! Matsuno backward
          forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
          CALL exner_hyb(ps, p3d, pks, pk, pkf)

          ! Calcul des tendances dynamiques:
          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.)

          ! integrations dynamique et traceurs:
          CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, &
               dteta, dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, &
               finvmaold, dtvr, leapf=.false.)
       end if

       IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
          ! Calcul des tendances physiques:

          forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
          CALL exner_hyb(ps, p3d, pks, pk, pkf)

          rdaym_ini = itau * dtvr / daysec
          rdayvrai = rdaym_ini + day_ini
          time = REAL(mod(itau, day_step)) / day_step + time_0
          IF (time > 1.) time = time - 1.

          CALL calfis(rdayvrai, time, ucov, vcov, teta, q, masse, ps, pk, &
               phis, phi, dudyn, dv, dq, w, dufi, dvfi, dtetafi, dqfi, dpfi, &
               lafin = itau + 1 == itaufin)

          ! 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
       CALL exner_hyb(ps, p3d, pks, pk, pkf)

       IF (MOD(itau + 1, idissip) == 0) THEN
          ! Dissipation horizontale et verticale des petites échelles

          ! calcul de l'énergie cinétique avant dissipation
          call covcont(llm, ucov, vcov, ucont, vcont)
          call enercin(vcov, ucov, vcont, ucont, ecin0)

          ! dissipation
          CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
          ucov = ucov + dudis
          vcov = vcov + dvdis

          ! On ajoute la tendance due à la transformation énergie
          ! cinétique en énergie thermique par la dissipation
          call covcont(llm, ucov, vcov, ucont, vcont)
          call enercin(vcov, ucov, vcont, ucont, ecin)
          dtetadis = dtetadis + (ecin0 - ecin) / pk
          teta = teta + dtetadis

          ! Calcul de la valeur moyenne aux pôles :
          forall (l = 1: llm)
             teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
                  / apoln
             teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
                  * teta(:iim, jjm + 1, l)) / apols
          END forall

          ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
          ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &
               / apols
       END IF

       IF (MOD(itau + 1, iperiod) == 0) THEN
          ! Écriture du fichier histoire moyenne:
          CALL writedynav(vcov, ucov, teta, pk, phi, q, masse, ps, phis, &
               time = itau + 1)
          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(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)

  END SUBROUTINE leapfrog

end module leapfrog_m
1	guez	3	module leapfrog_m
2
3			IMPLICIT NONE
4
5			contains
6
7	guez	23	SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)
8	guez	3
9	guez	67	! From dyn3d/leapfrog.F, version 1.6, 2005/04/13 08:58:34 revision 616
10	guez	27	! Authors: P. Le Van, L. Fairhead, F. Hourdin
11	guez	36	! Matsuno-leapfrog scheme.
12	guez	3
13	guez	37	use addfi_m, only: addfi
14	guez	40	use bilan_dyn_m, only: bilan_dyn
15			use caladvtrac_m, only: caladvtrac
16	guez	43	use caldyn_m, only: caldyn
17	guez	26	USE calfis_m, ONLY: calfis
18			USE comconst, ONLY: daysec, dtphys, dtvr
19	guez	29	USE comgeom, ONLY: aire_2d, apoln, apols
20	guez	66	USE disvert_m, ONLY: ap, bp
21	guez	57	USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
22	guez	68	iflag_phys, ok_guide, iecri
23	guez	29	USE dimens_m, ONLY: iim, jjm, llm, nqmx
24	guez	47	use dissip_m, only: dissip
25	guez	26	USE dynetat0_m, ONLY: day_ini
26	guez	27	use dynredem1_m, only: dynredem1
27	guez	26	USE exner_hyb_m, ONLY: exner_hyb
28	guez	27	use filtreg_m, only: filtreg
29	guez	69	use fluxstokenc_m, only: fluxstokenc
30	guez	43	use geopot_m, only: geopot
31	guez	26	USE guide_m, ONLY: guide
32			use inidissip_m, only: idissip
33	guez	32	use integrd_m, only: integrd
34	guez	55	use nr_util, only: assert
35	guez	26	USE pressure_var, ONLY: p3d
36	guez	28	USE temps, ONLY: itau_dyn
37	guez	56	use writedynav_m, only: writedynav
38	guez	69	use writehist_m, only: writehist
39	guez	3
40	guez	10	! Variables dynamiques:
41	guez	55	REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
42			REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
43	guez	43
44			REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
45			! potential temperature
46
47	guez	45	REAL, intent(inout):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol, en Pa
48	guez	70	REAL, intent(inout):: masse(:, :, :) ! (iim + 1, jjm + 1, llm) masse d'air
49	guez	69	REAL, intent(in):: phis(:, :) ! (iim + 1, jjm + 1) surface geopotential
50	guez	40
51			REAL, intent(inout):: q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nqmx)
52			! mass fractions of advected fields
53
54	guez	24	REAL, intent(in):: time_0
55	guez	10
56			! Variables local to the procedure:
57
58			! Variables dynamiques:
59
60	guez	70	REAL pks(iim + 1, jjm + 1) ! exner au sol
61	guez	29	REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
62	guez	65	REAL pkf(iim + 1, jjm + 1, llm) ! exner filtré au milieu des couches
63	guez	47	REAL phi(iim + 1, jjm + 1, llm) ! geopotential
64	guez	55	REAL w((iim + 1) * (jjm + 1), llm) ! vitesse verticale
65	guez	3
66	guez	55	! Variables dynamiques intermediaire pour le transport
67			! Flux de masse :
68			REAL pbaru((iim + 1) * (jjm + 1), llm), pbarv((iim + 1) * jjm, llm)
69	guez	3
70	guez	56	! Variables dynamiques au pas - 1
71	guez	55	REAL vcovm1(iim + 1, jjm, llm), ucovm1(iim + 1, jjm + 1, llm)
72	guez	29	REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
73	guez	67	REAL massem1(iim + 1, jjm + 1, llm)
74	guez	3
75	guez	56	! Tendances dynamiques
76	guez	55	REAL dv((iim + 1) * jjm, llm), dudyn((iim + 1) * (jjm + 1), llm)
77			REAL dteta(iim + 1, jjm + 1, llm), dq((iim + 1) * (jjm + 1), llm, nqmx)
78			real dp((iim + 1) * (jjm + 1))
79	guez	3
80	guez	56	! Tendances de la dissipation :
81	guez	55	REAL dvdis(iim + 1, jjm, llm), dudis(iim + 1, jjm + 1, llm)
82	guez	29	REAL dtetadis(iim + 1, jjm + 1, llm)
83	guez	3
84	guez	56	! Tendances physiques
85	guez	55	REAL dvfi((iim + 1) * jjm, llm), dufi((iim + 1) * (jjm + 1), llm)
86			REAL dtetafi(iim + 1, jjm + 1, llm), dqfi((iim + 1) * (jjm + 1), llm, nqmx)
87			real dpfi((iim + 1) * (jjm + 1))
88	guez	3
89	guez	56	! Variables pour le fichier histoire
90	guez	3
91	guez	22	INTEGER itau ! index of the time step of the dynamics, starts at 0
92	guez	27	INTEGER itaufin
93	guez	20	REAL time ! time of day, as a fraction of day length
94	guez	67	real finvmaold(iim + 1, jjm + 1, llm)
95	guez	33	INTEGER l
96	guez	3	REAL rdayvrai, rdaym_ini
97
98	guez	24	! Variables test conservation energie
99	guez	29	REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
100	guez	43
101	guez	55	REAL vcont((iim + 1) * jjm, llm), ucont((iim + 1) * (jjm + 1), llm)
102	guez	33	logical leapf
103			real dt
104	guez	3
105			!---------------------------------------------------
106
107			print *, "Call sequence information: leapfrog"
108	guez	55	call assert(shape(ucov) == (/iim + 1, jjm + 1, llm/), "leapfrog")
109	guez	3
110			itaufin = nday * day_step
111	guez	62	! "day_step" is a multiple of "iperiod", therefore so is "itaufin".
112	guez	30
113	guez	24	dq = 0.
114	guez	33
115	guez	3	! On initialise la pression et la fonction d'Exner :
116	guez	37	forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
117	guez	10	CALL exner_hyb(ps, p3d, pks, pk, pkf)
118	guez	3
119	guez	40	time_integration: do itau = 0, itaufin - 1
120	guez	33	leapf = mod(itau, iperiod) /= 0
121			if (leapf) then
122			dt = 2 * dtvr
123			else
124			! Matsuno
125			dt = dtvr
126			if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
127			call guide(itau, ucov, vcov, teta, q, masse, ps)
128			vcovm1 = vcov
129			ucovm1 = ucov
130			tetam1 = teta
131			massem1 = masse
132			psm1 = ps
133			finvmaold = masse
134	guez	64	CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE.)
135	guez	33	end if
136	guez	30
137			! Calcul des tendances dynamiques:
138	guez	70	CALL geopot(teta, pk, pks, phis, phi)
139	guez	30	CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
140	guez	47	dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
141			conser=MOD(itau, iconser)==0)
142	guez	30
143			! Calcul des tendances advection des traceurs (dont l'humidité)
144			CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
145	guez	33
146	guez	30	! Stokage du flux de masse pour traceurs offline:
147			IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
148			dtvr, itau)
149
150	guez	70	! Intégrations dynamique et traceurs:
151	guez	47	CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, dteta, &
152			dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, finvmaold, dt, &
153			leapf)
154	guez	30
155	guez	33	if (.not. leapf) then
156			! Matsuno backward
157	guez	37	forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
158	guez	33	CALL exner_hyb(ps, p3d, pks, pk, pkf)
159	guez	30
160	guez	27	! Calcul des tendances dynamiques:
161	guez	70	CALL geopot(teta, pk, pks, phis, phi)
162	guez	33	CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
163	guez	47	phi, dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
164			conser=.false.)
165	guez	3
166			! integrations dynamique et traceurs:
167	guez	47	CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, &
168			dteta, dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, &
169			finvmaold, dtvr, leapf=.false.)
170	guez	33	end if
171	guez	3
172	guez	33	IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
173	guez	69	! Calcul des tendances physiques:
174	guez	3
175	guez	37	forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
176	guez	33	CALL exner_hyb(ps, p3d, pks, pk, pkf)
177	guez	3
178	guez	33	rdaym_ini = itau * dtvr / daysec
179			rdayvrai = rdaym_ini + day_ini
180			time = REAL(mod(itau, day_step)) / day_step + time_0
181			IF (time > 1.) time = time - 1.
182	guez	3
183	guez	37	CALL calfis(rdayvrai, time, ucov, vcov, teta, q, masse, ps, pk, &
184	guez	47	phis, phi, dudyn, dv, dq, w, dufi, dvfi, dtetafi, dqfi, dpfi, &
185	guez	67	lafin = itau + 1 == itaufin)
186	guez	3
187	guez	69	! Ajout des tendances physiques:
188	guez	70	CALL addfi(nqmx, ucov, vcov, teta, q, ps, dufi, dvfi, dtetafi, &
189			dqfi, dpfi)
190	guez	33	ENDIF
191	guez	3
192	guez	37	forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
193	guez	33	CALL exner_hyb(ps, p3d, pks, pk, pkf)
194	guez	3
195	guez	33	IF (MOD(itau + 1, idissip) == 0) THEN
196	guez	55	! Dissipation horizontale et verticale des petites échelles
197	guez	3
198	guez	55	! calcul de l'énergie cinétique avant dissipation
199	guez	33	call covcont(llm, ucov, vcov, ucont, vcont)
200			call enercin(vcov, ucov, vcont, ucont, ecin0)
201	guez	3
202	guez	33	! dissipation
203			CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
204	guez	55	ucov = ucov + dudis
205			vcov = vcov + dvdis
206	guez	3
207	guez	55	! On ajoute la tendance due à la transformation énergie
208			! cinétique en énergie thermique par la dissipation
209	guez	33	call covcont(llm, ucov, vcov, ucont, vcont)
210			call enercin(vcov, ucov, vcont, ucont, ecin)
211	guez	56	dtetadis = dtetadis + (ecin0 - ecin) / pk
212	guez	55	teta = teta + dtetadis
213	guez	3
214	guez	33	! Calcul de la valeur moyenne aux pôles :
215			forall (l = 1: llm)
216			teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
217			/ apoln
218			teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
219			* teta(:iim, jjm + 1, l)) / apols
220			END forall
221	guez	3
222	guez	33	ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
223			ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &
224			/ apols
225			END IF
226	guez	3
227	guez	33	IF (MOD(itau + 1, iperiod) == 0) THEN
228	guez	40	! Écriture du fichier histoire moyenne:
229	guez	61	CALL writedynav(vcov, ucov, teta, pk, phi, q, masse, ps, phis, &
230			time = itau + 1)
231	guez	40	call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
232	guez	57	q(:, :, :, 1))
233	guez	33	ENDIF
234	guez	68
235			IF (MOD(itau + 1, iecri * day_step) == 0) THEN
236	guez	70	CALL geopot(teta, pk, pks, phis, phi)
237	guez	69	CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps)
238	guez	68	END IF
239	guez	40	end do time_integration
240	guez	3
241	guez	30	CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
242	guez	67	itau = itau_dyn + itaufin)
243	guez	30
244			! Calcul des tendances dynamiques:
245	guez	70	CALL geopot(teta, pk, pks, phis, phi)
246	guez	30	CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
247	guez	47	dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
248	guez	67	conser = MOD(itaufin, iconser) == 0)
249	guez	56
250	guez	3	END SUBROUTINE leapfrog
251
252			end module leapfrog_m