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
    ! Auteurs : P. Le Van, L. Fairhead, F. Hourdin

    USE dimens_m, ONLY : iim, llm, nqmx
    USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1
    USE comconst, ONLY : daysec, dtphys, dtvr
    USE comvert, ONLY : ap, bp
    USE conf_gcm_m, ONLY : day_step, iconser, idissip, iperiod, iphysiq, &
         nday, offline, periodav
    USE logic, ONLY : iflag_phys, ok_guide
    USE comgeom, ONLY : aire, apoln, apols
    USE temps, ONLY : dt, itaufin
    USE dynetat0_m, ONLY : day_ini
    USE iniprint, ONLY : prt_level
    USE com_io_dyn, ONLY : histaveid
    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 ps(ip1jmp1) ! pression au sol, en Pa

    REAL masse(ip1jmp1, llm) ! masse d'air
    REAL phis(ip1jmp1) ! geopotentiel au sol
    REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields
    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 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
    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.) &
            call guide(itau, ucov, vcov, teta, q, masse, ps)
       vcovm1 = vcov
       ucovm1 = ucov
       tetam1 = teta
       massem1 = masse
       psm1 = ps
       forward = .TRUE.
       leapf = .FALSE.
       dt = dtvr
       finvmaold = masse
       CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)

       do
          ! gestion des appels de la physique et des dissipations:
          apphys = MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0
          conser = MOD(itau, iconser) == 0
          apdiss = MOD(itau + 1, idissip) == 0

          ! 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)

          IF (forward .OR. leapf) THEN
             ! calcul des tendances advection des traceurs (dont l'humidite)
             CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
             IF (offline) THEN
                ! 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)

          IF (apphys) THEN
             ! calcul des tendances physiques:
             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

             ! 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)

          IF (apdiss) THEN
             ! dissipation horizontale et verticale des petites echelles:

             ! 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 = SUM(tppn) / apoln
                tps = SUM(tpps) / 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 = SUM(tppn) / apoln
             tps = SUM(tpps) / 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

          IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN
             ! ecriture du fichier histoire moyenne:
             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)
          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	guez	24	! Auteurs : P. Le Van, L. Fairhead, F. Hourdin
11	guez	3
12	guez	25	USE dimens_m, ONLY : iim, llm, nqmx
13			USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1
14			USE comconst, ONLY : daysec, dtphys, dtvr
15			USE comvert, ONLY : ap, bp
16			USE conf_gcm_m, ONLY : day_step, iconser, idissip, iperiod, iphysiq, &
17			nday, offline, periodav
18			USE logic, ONLY : iflag_phys, ok_guide
19			USE comgeom, ONLY : aire, apoln, apols
20			USE temps, ONLY : dt, itaufin
21			USE dynetat0_m, ONLY : day_ini
22			USE iniprint, ONLY : prt_level
23			USE com_io_dyn, ONLY : histaveid
24			USE calfis_m, ONLY : calfis
25			USE exner_hyb_m, ONLY : exner_hyb
26			USE guide_m, ONLY : guide
27			USE pression_m, ONLY : pression
28			USE pressure_var, ONLY : p3d
29	guez	3
30	guez	10	! Variables dynamiques:
31	guez	3	REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
32			REAL teta(ip1jmp1, llm) ! temperature potentielle
33	guez	10	REAL ps(ip1jmp1) ! pression au sol, en Pa
34	guez	25
35	guez	10	REAL masse(ip1jmp1, llm) ! masse d'air
36			REAL phis(ip1jmp1) ! geopotentiel au sol
37	guez	25	REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields
38	guez	24	REAL, intent(in):: time_0
39	guez	10
40			! Variables local to the procedure:
41
42			! Variables dynamiques:
43
44	guez	3	REAL pks(ip1jmp1) ! exner au sol
45			REAL pk(ip1jmp1, llm) ! exner au milieu des couches
46			REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches
47			REAL phi(ip1jmp1, llm) ! geopotential
48			REAL w(ip1jmp1, llm) ! vitesse verticale
49
50			! variables dynamiques intermediaire pour le transport
51			REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse
52
53			! variables dynamiques au pas - 1
54			REAL vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)
55			REAL tetam1(ip1jmp1, llm), psm1(ip1jmp1)
56			REAL massem1(ip1jmp1, llm)
57
58			! tendances dynamiques
59			REAL dv(ip1jm, llm), du(ip1jmp1, llm)
60			REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)
61
62			! tendances de la dissipation
63			REAL dvdis(ip1jm, llm), dudis(ip1jmp1, llm)
64			REAL dtetadis(ip1jmp1, llm)
65
66			! tendances physiques
67			REAL dvfi(ip1jm, llm), dufi(ip1jmp1, llm)
68			REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)
69
70			! variables pour le fichier histoire
71
72			REAL tppn(iim), tpps(iim), tpn, tps
73
74	guez	22	INTEGER itau ! index of the time step of the dynamics, starts at 0
75	guez	3	INTEGER iday ! jour julien
76	guez	20	REAL time ! time of day, as a fraction of day length
77	guez	10	real finvmaold(ip1jmp1, llm)
78	guez	3	LOGICAL :: lafin=.false.
79			INTEGER ij, l
80
81			REAL rdayvrai, rdaym_ini
82
83	guez	24	! Variables test conservation energie
84	guez	3	REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)
85			! Tendance de la temp. potentiel d (theta) / d t due a la
86			! tansformation d'energie cinetique en energie thermique
87			! cree par la dissipation
88			REAL dtetaecdt(ip1jmp1, llm)
89			REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)
90			CHARACTER*15 ztit
91	guez	10	INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.
92	guez	3
93	guez	10	logical:: dissip_conservative = .true.
94	guez	12	logical forward, leapf, apphys, conser, apdiss
95	guez	3
96			!---------------------------------------------------
97
98			print *, "Call sequence information: leapfrog"
99
100			itaufin = nday * day_step
101			itau = 0
102			iday = day_ini
103			time = time_0
104	guez	24	dq = 0.
105	guez	3	! On initialise la pression et la fonction d'Exner :
106	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
107			CALL exner_hyb(ps, p3d, pks, pk, pkf)
108	guez	3
109			! Debut de l'integration temporelle:
110	guez	10	outer_loop:do
111	guez	25	if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
112			call guide(itau, ucov, vcov, teta, q, masse, ps)
113			vcovm1 = vcov
114			ucovm1 = ucov
115			tetam1 = teta
116			massem1 = masse
117			psm1 = ps
118	guez	3	forward = .TRUE.
119			leapf = .FALSE.
120			dt = dtvr
121	guez	25	finvmaold = masse
122	guez	3	CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)
123
124			do
125			! gestion des appels de la physique et des dissipations:
126	guez	25	apphys = MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0
127			conser = MOD(itau, iconser) == 0
128			apdiss = MOD(itau + 1, idissip) == 0
129	guez	3
130			! calcul des tendances dynamiques:
131			CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
132			CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
133			conser, du, dv, dteta, dp, w, pbaru, pbarv, &
134			time + iday - day_ini)
135
136			IF (forward .OR. leapf) THEN
137	guez	25	! calcul des tendances advection des traceurs (dont l'humidite)
138	guez	10	CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
139	guez	3	IF (offline) THEN
140	guez	25	! Stokage du flux de masse pour traceurs off-line
141	guez	3	CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &
142			itau)
143			ENDIF
144			ENDIF
145
146			! integrations dynamique et traceurs:
147			CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
148	guez	12	dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &
149			finvmaold, leapf)
150	guez	3
151			IF (apphys) THEN
152	guez	25	! calcul des tendances physiques:
153	guez	3	IF (itau + 1 == itaufin) lafin = .TRUE.
154
155	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
156			CALL exner_hyb(ps, p3d, pks, pk, pkf)
157	guez	3
158			rdaym_ini = itau * dtvr / daysec
159			rdayvrai = rdaym_ini + day_ini
160
161			! Diagnostique de conservation de l'énergie : initialisation
162			IF (ip_ebil_dyn >= 1) THEN
163			ztit='bil dyn'
164	guez	10	CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
165			teta, q(:, :, 1), q(:, :, 2))
166	guez	3	ENDIF
167
168	guez	23	CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &
169	guez	10	masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
170	guez	13	dufi, dvfi, dtetafi, dqfi, dpfi)
171	guez	3
172			! ajout des tendances physiques:
173			CALL addfi(nqmx, dtphys, &
174			ucov, vcov, teta, q, ps, &
175			dufi, dvfi, dtetafi, dqfi, dpfi)
176
177			! Diagnostique de conservation de l'énergie : difference
178			IF (ip_ebil_dyn >= 1) THEN
179			ztit = 'bil phys'
180	guez	10	CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
181	guez	3	teta, q(:, :, 1), q(:, :, 2))
182			ENDIF
183			ENDIF
184
185	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
186			CALL exner_hyb(ps, p3d, pks, pk, pkf)
187	guez	3
188	guez	25	IF (apdiss) THEN
189			! dissipation horizontale et verticale des petites echelles:
190	guez	3
191			! calcul de l'energie cinetique avant dissipation
192			call covcont(llm, ucov, vcov, ucont, vcont)
193			call enercin(vcov, ucov, vcont, ucont, ecin0)
194
195			! dissipation
196	guez	10	CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
197	guez	3	ucov=ucov + dudis
198			vcov=vcov + dvdis
199
200			if (dissip_conservative) then
201			! On rajoute la tendance due a la transform. Ec -> E
202			! therm. cree lors de la dissipation
203			call covcont(llm, ucov, vcov, ucont, vcont)
204			call enercin(vcov, ucov, vcont, ucont, ecin)
205			dtetaecdt= (ecin0 - ecin) / pk
206			dtetadis=dtetadis + dtetaecdt
207			endif
208			teta=teta + dtetadis
209
210			! Calcul de la valeur moyenne, unique de h aux poles .....
211			DO l = 1, llm
212			DO ij = 1, iim
213			tppn(ij) = aire(ij) * teta(ij, l)
214			tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)
215			ENDDO
216	guez	25	tpn = SUM(tppn) / apoln
217			tps = SUM(tpps) / apols
218	guez	3
219			DO ij = 1, iip1
220			teta(ij, l) = tpn
221			teta(ij + ip1jm, l) = tps
222			ENDDO
223			ENDDO
224
225			DO ij = 1, iim
226			tppn(ij) = aire(ij) * ps(ij)
227			tpps(ij) = aire(ij + ip1jm) * ps(ij + ip1jm)
228			ENDDO
229	guez	25	tpn = SUM(tppn) / apoln
230			tps = SUM(tpps) / apols
231	guez	3
232			DO ij = 1, iip1
233			ps(ij) = tpn
234			ps(ij + ip1jm) = tps
235			ENDDO
236			END IF
237
238			! fin de l'intégration dynamique et physique pour le pas "itau"
239			! préparation du pas d'intégration suivant
240
241			! schema matsuno + leapfrog
242			IF (forward .OR. leapf) THEN
243			itau = itau + 1
244			iday = day_ini + itau / day_step
245			time = REAL(itau - (iday - day_ini) * day_step) / day_step &
246			+ time_0
247			IF (time > 1.) THEN
248			time = time - 1.
249			iday = iday + 1
250			ENDIF
251			ENDIF
252
253	guez	24	IF (itau == itaufin + 1) exit outer_loop
254	guez	3
255			IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN
256	guez	25	! ecriture du fichier histoire moyenne:
257	guez	3	CALL writedynav(histaveid, nqmx, itau, vcov, &
258			ucov, teta, pk, phi, q, masse, ps, phis)
259			call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
260			ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
261			ENDIF
262
263			IF (itau == itaufin) THEN
264	guez	23	CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps)
265	guez	3	ENDIF
266
267			! gestion de l'integration temporelle:
268			IF (MOD(itau, iperiod) == 0) exit
269			IF (MOD(itau - 1, iperiod) == 0) THEN
270			IF (forward) THEN
271			! fin du pas forward et debut du pas backward
272			forward = .FALSE.
273			leapf = .FALSE.
274			ELSE
275			! fin du pas backward et debut du premier pas leapfrog
276			leapf = .TRUE.
277			dt = 2. * dtvr
278			END IF
279			ELSE
280	guez	20	! pas leapfrog
281	guez	3	leapf = .TRUE.
282			dt = 2. * dtvr
283			END IF
284			end do
285	guez	10	end do outer_loop
286	guez	3
287			END SUBROUTINE leapfrog
288
289			end module leapfrog_m