libf/dyn3d/leapfrog.f90

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	guez	3	module leapfrog_m
2
3			IMPLICIT NONE
4
5			contains
6
7	guez	23	SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)
8	guez	3
9			! From dyn3d/leapfrog.F, version 1.6 2005/04/13 08:58:34
10
11	guez	24	! Version du 10/01/98, avec coordonnées verticales hybrides, avec
12			! nouveaux opérateurs dissipation "*" (gradiv2, divgrad2, nxgraro2)
13	guez	3
14	guez	24	! Auteurs : P. Le Van, L. Fairhead, F. Hourdin
15			! Objet: nouvelle grille
16	guez	3
17	guez	24	! 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	guez	3
21	guez	24	! Possibilité de choisir le schéma pour l'advection de
22			! q, en modifiant iadv dans "traceur.def".
23	guez	3
24	guez	24	! Pour Van-Leer + vapeur d'eau saturée, iadv(1)=4.
25	guez	3	! Pour Van-Leer iadv=10
26
27	guez	23	use dimens_m, only: iim, llm, nqmx
28	guez	24	use paramet_m, only: ip1jmp1, ip1jm, ijp1llm, jjp1, iip1
29	guez	3	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	guez	12	use logic, only: ok_guide, iflag_phys
34	guez	3	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	guez	10	use pressure_var, only: p3d
45	guez	3
46	guez	10	! Variables dynamiques:
47	guez	3	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	guez	10	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	guez	24	REAL, intent(in):: time_0
55	guez	10
56			! Variables local to the procedure:
57
58			! Variables dynamiques:
59
60	guez	3	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	guez	22	INTEGER itau ! index of the time step of the dynamics, starts at 0
91	guez	3	INTEGER iday ! jour julien
92	guez	20	REAL time ! time of day, as a fraction of day length
93	guez	3
94			REAL SSUM
95	guez	10	real finvmaold(ip1jmp1, llm)
96	guez	3
97			LOGICAL :: lafin=.false.
98			INTEGER ij, l
99
100			REAL rdayvrai, rdaym_ini
101
102	guez	24	! Variables test conservation energie
103	guez	3	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	guez	10	INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.
111	guez	3
112	guez	10	logical:: dissip_conservative = .true.
113	guez	12	logical forward, leapf, apphys, conser, apdiss
114	guez	3
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	guez	24	dq = 0.
129	guez	3	! On initialise la pression et la fonction d'Exner :
130	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
131			CALL exner_hyb(ps, p3d, pks, pk, pkf)
132	guez	3
133			! Debut de l'integration temporelle:
134	guez	10	outer_loop:do
135	guez	24	if (ok_guide.and.(itaufin - itau - 1) * dtvr > 21600.) then
136	guez	3	call guide(itau, ucov, vcov, teta, q, masse, ps)
137			else
138			IF (prt_level > 9) print *, &
139	guez	24	'Attention : on ne guide pas les 6 dernières heures.'
140	guez	3	endif
141
142	guez	24	CALL SCOPY(ip1jm * llm, vcov, 1, vcovm1, 1)
143	guez	3	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	guez	10	CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
178	guez	3	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	guez	12	dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &
188			finvmaold, leapf)
189	guez	3
190			! calcul des tendances physiques:
191
192			IF (apphys) THEN
193			IF (itau + 1 == itaufin) lafin = .TRUE.
194
195	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
196			CALL exner_hyb(ps, p3d, pks, pk, pkf)
197	guez	3
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	guez	10	CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
207			teta, q(:, :, 1), q(:, :, 2))
208	guez	3	ENDIF
209
210	guez	23	CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &
211	guez	10	masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
212	guez	13	dufi, dvfi, dtetafi, dqfi, dpfi)
213	guez	3
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	guez	10	CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
223	guez	3	teta, q(:, :, 1), q(:, :, 2))
224			ENDIF
225			ENDIF
226
227	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
228			CALL exner_hyb(ps, p3d, pks, pk, pkf)
229	guez	3
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	guez	10	CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
239	guez	3	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	guez	24	IF (itau == itaufin + 1) exit outer_loop
298	guez	3
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	guez	23	CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps)
312	guez	3	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	guez	20	! pas leapfrog
330	guez	3	leapf = .TRUE.
331			dt = 2. * dtvr
332			END IF
333			end do
334	guez	10	end do outer_loop
335	guez	3
336			END SUBROUTINE leapfrog
337
338			end module leapfrog_m