Sources/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, 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, llm, nqmx
    use paramet_m, only: ip1jmp1, ip1jm, ijmllm, 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 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 itaufinp1
    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

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