trunk/dyn3d/leapfrog.f

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, clesphy0, &
       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, apdiss, apphys, conser, forward, iflag_phys, &
         leapf, statcl
    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 nq
    REAL clesphy0(:)

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

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

    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.
          statcl = .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)

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