libf/dyn3d/leapfrog.f90

module leapfrog_m

  ! This module is clean: no C preprocessor directive, no include line.

  IMPLICIT NONE

contains

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

    ! From dyn3d/leapfrog.F, version 1.6 2005/04/13 08:58:34

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

    ! ... Possibilité de choisir le schéma 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 

    use dimens_m, only: iim, jjm, llm, nqmx
    use paramet_m, only: ip1jmp1, ip1jm, ijmllm, ijp1llm, jjp1, iip1, iip2
    use comconst, only: dtvr, daysec, dtphys
    use comvert, only: ap, bp
    use conf_gcm_m, only: day_step, iconser, idissip, iphysiq, iperiod, nday, &
         offline, periodav
    use logic, only: ok_guide, iflag_phys
    use comgeom
    use serre
    use temps, only: itaufin, day_ini, dt
    use iniprint, only: prt_level
    use com_io_dyn
    use ener
    use calfis_m, only: calfis
    use exner_hyb_m, only: exner_hyb
    use guide_m, only: guide
    use pression_m, only: pression
    use pressure_var, only: p3d

    integer, intent(in):: nq

    ! Variables dynamiques:
    REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
    REAL teta(ip1jmp1, llm) ! temperature potentielle 
    REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields
    REAL ps(ip1jmp1) ! pression au sol, en Pa
    REAL masse(ip1jmp1, llm) ! masse d'air
    REAL phis(ip1jmp1) ! geopotentiel au sol

    REAL time_0

    ! Variables local to the procedure:

    ! Variables dynamiques:

    REAL pks(ip1jmp1) ! exner au sol
    REAL pk(ip1jmp1, llm) ! exner au milieu des couches
    REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches
    REAL phi(ip1jmp1, llm) ! geopotential
    REAL w(ip1jmp1, llm) ! vitesse verticale

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

    ! variables pour le fichier histoire

    REAL tppn(iim), tpps(iim), tpn, tps

    INTEGER itau, itaufinp1
    INTEGER iday ! jour julien
    REAL time ! Heure de la journee en fraction d'1 jour

    REAL SSUM
    real finvmaold(ip1jmp1, llm)

    LOGICAL :: lafin=.false.
    INTEGER ij, l

    REAL rdayvrai, rdaym_ini
    LOGICAL:: callinigrads = .true.

    !+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)
    CHARACTER*15 ztit
    INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.

    logical:: dissip_conservative = .true.
    LOGICAL:: prem = .true.
    logical forward, leapf, apphys, conser, apdiss

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

    print *, "Call sequence information: leapfrog"

    itaufin = nday * day_step
    itaufinp1 = itaufin + 1

    itau = 0
    iday = day_ini
    time = time_0
    IF (time > 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, p3d)
    CALL exner_hyb(ps, p3d, pks, pk, pkf)

    ! Debut de l'integration temporelle:
    outer_loop:do
       if (ok_guide.and.(itaufin - itau - 1) * dtvr > 21600) then
          call guide(itau, ucov, vcov, teta, q, masse, ps)
       else
          IF (prt_level > 9) print *, &
               'Attention : on ne guide pas les 6 dernieres heures.'
       endif

       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

       CALL SCOPY(ijp1llm, masse, 1, finvmaold, 1)
       CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)

       do
          ! gestion des appels de la physique et des dissipations:

          apphys = .FALSE.
          conser = .FALSE.
          apdiss = .FALSE.

          IF (MOD(itau, iconser) == 0) conser = .TRUE.
          IF (MOD(itau + 1, idissip) == 0) apdiss = .TRUE.
          IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) apphys=.TRUE.

          ! calcul des tendances dynamiques:

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

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

          ! calcul des tendances advection des traceurs (dont l'humidite)

          IF (forward .OR. leapf) THEN
             CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
             IF (offline) THEN
                !maf stokage du flux de masse pour traceurs OFF-LINE
                CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &
                     itau)
             ENDIF
          ENDIF

          ! integrations dynamique et traceurs:
          CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
               dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &
               finvmaold, leapf)

          ! calcul des tendances physiques:

          IF (apphys) THEN
             IF (itau + 1 == itaufin) lafin = .TRUE.

             CALL pression(ip1jmp1, ap, bp, ps, p3d)
             CALL exner_hyb(ps, p3d, pks, pk, pkf)

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

             ! Interface avec les routines de phylmd (phymars ...)

             ! Diagnostique de conservation de l'énergie : initialisation
             IF (ip_ebil_dyn >= 1) THEN 
                ztit='bil dyn'
                CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
                     teta, q(:, :, 1), q(:, :, 2))
             ENDIF

             CALL calfis(nq, lafin, rdayvrai, time, ucov, vcov, teta, q, &
                  masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
                  dufi, dvfi, dtetafi, dqfi, dpfi)

             ! ajout des tendances physiques:
             CALL addfi(nqmx, dtphys, &
                  ucov, vcov, teta, q, ps, &
                  dufi, dvfi, dtetafi, dqfi, dpfi)

             ! Diagnostique de conservation de l'énergie : difference
             IF (ip_ebil_dyn >= 1) THEN 
                ztit = 'bil phys'
                CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
                     teta, q(:, :, 1), q(:, :, 2))
             ENDIF
          ENDIF

          CALL pression(ip1jmp1, ap, bp, ps, p3d)
          CALL exner_hyb(ps, p3d, pks, pk, pkf)

          ! 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, p3d, dvdis, dudis, dtetadis)
             ucov=ucov + dudis
             vcov=vcov + dvdis

             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
                dtetadis=dtetadis + dtetaecdt
             endif
             teta=teta + dtetadis

             ! 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

             DO ij = 1, iim
                tppn(ij) = aire(ij) * ps(ij)
                tpps(ij) = aire(ij + ip1jm) * ps(ij + ip1jm)
             ENDDO
             tpn = SSUM(iim, tppn, 1) / apoln
             tps = SSUM(iim, tpps, 1) / apols

             DO ij = 1, iip1
                ps(ij) = tpn
                ps(ij + ip1jm) = tps
             ENDDO

          END IF

          ! fin de l'intégration dynamique et physique pour le pas "itau"
          ! préparation du pas d'intégration suivant

          ! 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 > 1.) THEN
                time = time - 1.
                iday = iday + 1
             ENDIF
          ENDIF

          IF (itau == itaufinp1) exit outer_loop

          ! ecriture du fichier histoire moyenne:

          ! Comment out the following calls when you do not want the output
          ! files "dyn_hist_ave.nc" and "dynzon.nc"
          IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN
             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

          IF (itau == itaufin) THEN
             CALL dynredem1("restart.nc", 0., vcov, ucov, teta, q, masse, ps)
             CLOSE(99)
          ENDIF

          ! gestion de l'integration temporelle:

          IF (MOD(itau, iperiod) == 0) exit
          IF (MOD(itau - 1, iperiod) == 0) THEN
             IF (forward) THEN
                ! fin du pas forward et debut du pas backward
                forward = .FALSE.
                leapf = .FALSE.
             ELSE
                ! fin du pas backward et debut du premier pas leapfrog
                leapf = .TRUE.
                dt = 2. * dtvr
             END IF
          ELSE
             ! ...... pas leapfrog .....
             leapf = .TRUE.
             dt = 2. * dtvr
          END IF
       end do
    end do outer_loop

  END SUBROUTINE leapfrog

end module leapfrog_m
1	guez	3	module leapfrog_m
2
3			! This module is clean: no C preprocessor directive, no include line.
4
5			IMPLICIT NONE
6
7			contains
8
9	guez	13	SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, nq, q, time_0)
10	guez	3
11			! From dyn3d/leapfrog.F, version 1.6 2005/04/13 08:58:34
12
13			! Version du 10/01/98, avec coordonnees verticales hybrides, avec
14			! nouveaux operat. dissipation * (gradiv2, divgrad2, nxgraro2)
15
16			! Auteur: P. Le Van /L. Fairhead/F.Hourdin
17			! Objet:
18			! GCM LMD nouvelle grille
19
20			! ... Dans inigeom, nouveaux calculs pour les elongations cu, cv
21			! et possibilite d'appeler une fonction f(y) a derivee tangente
22			! hyperbolique a la place de la fonction a derivee sinusoidale.
23
24	guez	10	! ... Possibilité de choisir le schéma pour l'advection de
25	guez	3	! q, en modifiant iadv dans "traceur.def" (10/02) .
26
27			! Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron, 10/99)
28			! Pour Van-Leer iadv=10
29
30	guez	10	use dimens_m, only: iim, jjm, llm, nqmx
31			use paramet_m, only: ip1jmp1, ip1jm, ijmllm, ijp1llm, jjp1, iip1, iip2
32	guez	3	use comconst, only: dtvr, daysec, dtphys
33			use comvert, only: ap, bp
34			use conf_gcm_m, only: day_step, iconser, idissip, iphysiq, iperiod, nday, &
35			offline, periodav
36	guez	12	use logic, only: ok_guide, iflag_phys
37	guez	3	use comgeom
38			use serre
39			use temps, only: itaufin, day_ini, dt
40			use iniprint, only: prt_level
41			use com_io_dyn
42			use ener
43			use calfis_m, only: calfis
44			use exner_hyb_m, only: exner_hyb
45			use guide_m, only: guide
46			use pression_m, only: pression
47	guez	10	use pressure_var, only: p3d
48	guez	3
49	guez	18	integer, intent(in):: nq
50	guez	3
51	guez	10	! Variables dynamiques:
52	guez	3	REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
53			REAL teta(ip1jmp1, llm) ! temperature potentielle
54			REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields
55	guez	10	REAL ps(ip1jmp1) ! pression au sol, en Pa
56			REAL masse(ip1jmp1, llm) ! masse d'air
57			REAL phis(ip1jmp1) ! geopotentiel au sol
58
59			REAL time_0
60
61			! Variables local to the procedure:
62
63			! Variables dynamiques:
64
65	guez	3	REAL pks(ip1jmp1) ! exner au sol
66			REAL pk(ip1jmp1, llm) ! exner au milieu des couches
67			REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches
68			REAL phi(ip1jmp1, llm) ! geopotential
69			REAL w(ip1jmp1, llm) ! vitesse verticale
70
71			! variables dynamiques intermediaire pour le transport
72			REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse
73
74			! variables dynamiques au pas - 1
75			REAL vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)
76			REAL tetam1(ip1jmp1, llm), psm1(ip1jmp1)
77			REAL massem1(ip1jmp1, llm)
78
79			! tendances dynamiques
80			REAL dv(ip1jm, llm), du(ip1jmp1, llm)
81			REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)
82
83			! tendances de la dissipation
84			REAL dvdis(ip1jm, llm), dudis(ip1jmp1, llm)
85			REAL dtetadis(ip1jmp1, llm)
86
87			! tendances physiques
88			REAL dvfi(ip1jm, llm), dufi(ip1jmp1, llm)
89			REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)
90
91			! variables pour le fichier histoire
92
93			REAL tppn(iim), tpps(iim), tpn, tps
94
95			INTEGER itau, itaufinp1
96			INTEGER iday ! jour julien
97			REAL time ! Heure de la journee en fraction d'1 jour
98
99			REAL SSUM
100	guez	10	real finvmaold(ip1jmp1, llm)
101	guez	3
102			LOGICAL :: lafin=.false.
103			INTEGER ij, l
104
105			REAL rdayvrai, rdaym_ini
106	guez	10	LOGICAL:: callinigrads = .true.
107	guez	3
108			!+jld variables test conservation energie
109			REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)
110			! Tendance de la temp. potentiel d (theta) / d t due a la
111			! tansformation d'energie cinetique en energie thermique
112			! cree par la dissipation
113			REAL dtetaecdt(ip1jmp1, llm)
114			REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)
115			CHARACTER*15 ztit
116	guez	10	INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.
117	guez	3
118	guez	10	logical:: dissip_conservative = .true.
119			LOGICAL:: prem = .true.
120	guez	12	logical forward, leapf, apphys, conser, apdiss
121	guez	3
122			!---------------------------------------------------
123
124			print *, "Call sequence information: leapfrog"
125
126			itaufin = nday * day_step
127			itaufinp1 = itaufin + 1
128
129			itau = 0
130			iday = day_ini
131			time = time_0
132			IF (time > 1.) THEN
133			time = time - 1.
134			iday = iday + 1
135			ENDIF
136
137			! On initialise la pression et la fonction d'Exner :
138			dq=0.
139	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
140			CALL exner_hyb(ps, p3d, pks, pk, pkf)
141	guez	3
142			! Debut de l'integration temporelle:
143	guez	10	outer_loop:do
144	guez	3	if (ok_guide.and.(itaufin - itau - 1) * dtvr > 21600) then
145			call guide(itau, ucov, vcov, teta, q, masse, ps)
146			else
147			IF (prt_level > 9) print *, &
148			'Attention : on ne guide pas les 6 dernieres heures.'
149			endif
150
151			CALL SCOPY(ijmllm, vcov, 1, vcovm1, 1)
152			CALL SCOPY(ijp1llm, ucov, 1, ucovm1, 1)
153			CALL SCOPY(ijp1llm, teta, 1, tetam1, 1)
154			CALL SCOPY(ijp1llm, masse, 1, massem1, 1)
155			CALL SCOPY(ip1jmp1, ps, 1, psm1, 1)
156
157			forward = .TRUE.
158			leapf = .FALSE.
159			dt = dtvr
160
161			CALL SCOPY(ijp1llm, masse, 1, finvmaold, 1)
162			CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)
163
164			do
165			! gestion des appels de la physique et des dissipations:
166
167			apphys = .FALSE.
168			conser = .FALSE.
169			apdiss = .FALSE.
170
171			IF (MOD(itau, iconser) == 0) conser = .TRUE.
172			IF (MOD(itau + 1, idissip) == 0) apdiss = .TRUE.
173			IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) apphys=.TRUE.
174
175			! calcul des tendances dynamiques:
176
177			CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
178
179			CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
180			conser, du, dv, dteta, dp, w, pbaru, pbarv, &
181			time + iday - day_ini)
182
183			! calcul des tendances advection des traceurs (dont l'humidite)
184
185			IF (forward .OR. leapf) THEN
186	guez	10	CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
187	guez	3	IF (offline) THEN
188			!maf stokage du flux de masse pour traceurs OFF-LINE
189			CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &
190			itau)
191			ENDIF
192			ENDIF
193
194			! integrations dynamique et traceurs:
195			CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
196	guez	12	dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &
197			finvmaold, leapf)
198	guez	3
199			! calcul des tendances physiques:
200
201			IF (apphys) THEN
202			IF (itau + 1 == itaufin) lafin = .TRUE.
203
204	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
205			CALL exner_hyb(ps, p3d, pks, pk, pkf)
206	guez	3
207			rdaym_ini = itau * dtvr / daysec
208			rdayvrai = rdaym_ini + day_ini
209
210			! Interface avec les routines de phylmd (phymars ...)
211
212			! Diagnostique de conservation de l'énergie : initialisation
213			IF (ip_ebil_dyn >= 1) THEN
214			ztit='bil dyn'
215	guez	10	CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
216			teta, q(:, :, 1), q(:, :, 2))
217	guez	3	ENDIF
218
219			CALL calfis(nq, lafin, rdayvrai, time, ucov, vcov, teta, q, &
220	guez	10	masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
221	guez	13	dufi, dvfi, dtetafi, dqfi, dpfi)
222	guez	3
223			! ajout des tendances physiques:
224			CALL addfi(nqmx, dtphys, &
225			ucov, vcov, teta, q, ps, &
226			dufi, dvfi, dtetafi, dqfi, dpfi)
227
228			! Diagnostique de conservation de l'énergie : difference
229			IF (ip_ebil_dyn >= 1) THEN
230			ztit = 'bil phys'
231	guez	10	CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
232	guez	3	teta, q(:, :, 1), q(:, :, 2))
233			ENDIF
234			ENDIF
235
236	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
237			CALL exner_hyb(ps, p3d, pks, pk, pkf)
238	guez	3
239			! dissipation horizontale et verticale des petites echelles:
240
241			IF (apdiss) THEN
242			! calcul de l'energie cinetique avant dissipation
243			call covcont(llm, ucov, vcov, ucont, vcont)
244			call enercin(vcov, ucov, vcont, ucont, ecin0)
245
246			! dissipation
247	guez	10	CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
248	guez	3	ucov=ucov + dudis
249			vcov=vcov + dvdis
250
251			if (dissip_conservative) then
252			! On rajoute la tendance due a la transform. Ec -> E
253			! therm. cree lors de la dissipation
254			call covcont(llm, ucov, vcov, ucont, vcont)
255			call enercin(vcov, ucov, vcont, ucont, ecin)
256			dtetaecdt= (ecin0 - ecin) / pk
257			dtetadis=dtetadis + dtetaecdt
258			endif
259			teta=teta + dtetadis
260
261			! Calcul de la valeur moyenne, unique de h aux poles .....
262
263			DO l = 1, llm
264			DO ij = 1, iim
265			tppn(ij) = aire(ij) * teta(ij, l)
266			tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)
267			ENDDO
268			tpn = SSUM(iim, tppn, 1) / apoln
269			tps = SSUM(iim, tpps, 1) / apols
270
271			DO ij = 1, iip1
272			teta(ij, l) = tpn
273			teta(ij + ip1jm, l) = tps
274			ENDDO
275			ENDDO
276
277			DO ij = 1, iim
278			tppn(ij) = aire(ij) * ps(ij)
279			tpps(ij) = aire(ij + ip1jm) * ps(ij + ip1jm)
280			ENDDO
281			tpn = SSUM(iim, tppn, 1) / apoln
282			tps = SSUM(iim, tpps, 1) / apols
283
284			DO ij = 1, iip1
285			ps(ij) = tpn
286			ps(ij + ip1jm) = tps
287			ENDDO
288
289			END IF
290
291			! fin de l'intégration dynamique et physique pour le pas "itau"
292			! préparation du pas d'intégration suivant
293
294			! schema matsuno + leapfrog
295			IF (forward .OR. leapf) THEN
296			itau = itau + 1
297			iday = day_ini + itau / day_step
298			time = REAL(itau - (iday - day_ini) * day_step) / day_step &
299			+ time_0
300			IF (time > 1.) THEN
301			time = time - 1.
302			iday = iday + 1
303			ENDIF
304			ENDIF
305
306	guez	10	IF (itau == itaufinp1) exit outer_loop
307	guez	3
308			! ecriture du fichier histoire moyenne:
309
310			! Comment out the following calls when you do not want the output
311			! files "dyn_hist_ave.nc" and "dynzon.nc"
312			IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN
313			CALL writedynav(histaveid, nqmx, itau, vcov, &
314			ucov, teta, pk, phi, q, masse, ps, phis)
315			call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
316			ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
317			ENDIF
318
319			IF (itau == itaufin) THEN
320	guez	5	CALL dynredem1("restart.nc", 0., vcov, ucov, teta, q, masse, ps)
321	guez	3	CLOSE(99)
322			ENDIF
323
324			! gestion de l'integration temporelle:
325
326			IF (MOD(itau, iperiod) == 0) exit
327			IF (MOD(itau - 1, iperiod) == 0) THEN
328			IF (forward) THEN
329			! fin du pas forward et debut du pas backward
330			forward = .FALSE.
331			leapf = .FALSE.
332			ELSE
333			! fin du pas backward et debut du premier pas leapfrog
334			leapf = .TRUE.
335			dt = 2. * dtvr
336			END IF
337			ELSE
338			! ...... pas leapfrog .....
339			leapf = .TRUE.
340			dt = 2. * dtvr
341			END IF
342			end do
343	guez	10	end do outer_loop
344	guez	3
345			END SUBROUTINE leapfrog
346
347			end module leapfrog_m