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 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	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, time_0)
10
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	! ... Possibilité de choisir le schéma pour l'advection de
25	! 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	use dimens_m, only: iim, jjm, llm, nqmx
31	use paramet_m, only: ip1jmp1, ip1jm, ijmllm, ijp1llm, jjp1, iip1, iip2
32	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	use logic, only: ok_guide, iflag_phys
37	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	use pressure_var, only: p3d
48
49	integer nq
50
51	! Variables dynamiques:
52	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	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	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	real finvmaold(ip1jmp1, llm)
101
102	LOGICAL :: lafin=.false.
103	INTEGER ij, l
104
105	REAL rdayvrai, rdaym_ini
106	LOGICAL:: callinigrads = .true.
107
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	INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.
117
118	logical:: dissip_conservative = .true.
119	LOGICAL:: prem = .true.
120	logical forward, leapf, apphys, conser, apdiss
121
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	CALL pression(ip1jmp1, ap, bp, ps, p3d)
140	CALL exner_hyb(ps, p3d, pks, pk, pkf)
141
142	! Debut de l'integration temporelle:
143	outer_loop:do
144	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	CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
187	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	dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &
197	finvmaold, leapf)
198
199	! calcul des tendances physiques:
200
201	IF (apphys) THEN
202	IF (itau + 1 == itaufin) lafin = .TRUE.
203
204	CALL pression(ip1jmp1, ap, bp, ps, p3d)
205	CALL exner_hyb(ps, p3d, pks, pk, pkf)
206
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	CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
216	teta, q(:, :, 1), q(:, :, 2))
217	ENDIF
218
219	CALL calfis(nq, lafin, rdayvrai, time, ucov, vcov, teta, q, &
220	masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
221	dufi, dvfi, dtetafi, dqfi, dpfi)
222
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	CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
232	teta, q(:, :, 1), q(:, :, 2))
233	ENDIF
234	ENDIF
235
236	CALL pression(ip1jmp1, ap, bp, ps, p3d)
237	CALL exner_hyb(ps, p3d, pks, pk, pkf)
238
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	CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
248	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	IF (itau == itaufinp1) exit outer_loop
307
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	CALL dynredem1("restart.nc", 0., vcov, ucov, teta, q, masse, ps)
321	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	end do outer_loop
344
345	END SUBROUTINE leapfrog
346
347	end module leapfrog_m