trunk/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

    ! Version du 10/01/98, avec coordonnées verticales hybrides, avec
    ! nouveaux opérateurs dissipation "*" (gradiv2, divgrad2, nxgraro2)

    ! Auteurs : P. Le Van, L. Fairhead, F. Hourdin
    ! Objet: nouvelle grille

    ! Dans "inigeom", nouveaux calculs pour les élongations cu, cv
    ! et possibilité d'appeler une fonction f(y) à dérivée tangente
    ! hyperbolique à la place de la fonction à dérivée sinusoïdale.

    ! Possibilité de choisir le schéma pour l'advection de
    ! q, en modifiant iadv dans "traceur.def".

    ! Pour Van-Leer + vapeur d'eau saturée, iadv(1)=4.
    ! Pour Van-Leer iadv=10 

    use dimens_m, only: iim, llm, nqmx
    use paramet_m, only: ip1jmp1, ip1jm, ijp1llm, jjp1, iip1
    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

    ! 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, intent(in):: 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 ! index of the time step of the dynamics, starts at 0
    INTEGER iday ! jour julien
    REAL time ! time of day, as a fraction of day length

    REAL SSUM
    real finvmaold(ip1jmp1, llm)

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

    REAL rdayvrai, rdaym_ini

    ! 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 forward, leapf, apphys, conser, apdiss

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

    print *, "Call sequence information: leapfrog"

    itaufin = nday * day_step
    itau = 0
    iday = day_ini
    time = time_0
    IF (time > 1.) THEN
       time = time - 1.
       iday = iday + 1
    ENDIF

    dq = 0.
    ! On initialise la pression et la fonction d'Exner :
    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 dernières heures.'
       endif

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