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, 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	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	! ... Possibilité de choisir le schéma pour l'advection de
26	! 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	use dimens_m, only: iim, jjm, llm, nqmx
32	use paramet_m, only: ip1jmp1, ip1jm, ijmllm, ijp1llm, jjp1, iip1, iip2
33	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	use pressure_var, only: p3d
50
51	integer nq
52	REAL clesphy0(:)
53
54	! Variables dynamiques:
55	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	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	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	real finvmaold(ip1jmp1, llm)
104
105	LOGICAL :: lafin=.false.
106	INTEGER ij, l
107
108	REAL rdayvrai, rdaym_ini
109	LOGICAL:: callinigrads = .true.
110
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	INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.
120
121	logical:: dissip_conservative = .true.
122	LOGICAL:: prem = .true.
123
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	CALL pression(ip1jmp1, ap, bp, ps, p3d)
142	CALL exner_hyb(ps, p3d, pks, pk, pkf)
143
144	! Debut de l'integration temporelle:
145	outer_loop:do
146	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	CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
190	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	CALL pression(ip1jmp1, ap, bp, ps, p3d)
207	CALL exner_hyb(ps, p3d, pks, pk, pkf)
208
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	CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
218	teta, q(:, :, 1), q(:, :, 2))
219	ENDIF
220
221	CALL calfis(nq, lafin, rdayvrai, time, ucov, vcov, teta, q, &
222	masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
223	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	CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
234	teta, q(:, :, 1), q(:, :, 2))
235	ENDIF
236	ENDIF
237
238	CALL pression(ip1jmp1, ap, bp, ps, p3d)
239	CALL exner_hyb(ps, p3d, pks, pk, pkf)
240
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	CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
250	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	IF (itau == itaufinp1) exit outer_loop
309
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	CALL dynredem1("restart.nc", 0., vcov, ucov, teta, q, masse, ps)
323	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	end do outer_loop
346
347	END SUBROUTINE leapfrog
348
349	end module leapfrog_m