libf/dyn3d/leapfrog.f90

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, 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

    integer nq
    REAL, intent(in):: 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.
    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(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	guez	3	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	guez	10	! ... Possibilité de choisir le schéma pour l'advection de
26	guez	3	! 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	guez	10	use dimens_m, only: iim, jjm, llm, nqmx
32			use paramet_m, only: ip1jmp1, ip1jm, ijmllm, ijp1llm, jjp1, iip1, iip2
33	guez	3	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	guez	12	use logic, only: ok_guide, iflag_phys
38	guez	3	use comgeom
39			use serre
40			use temps, only: itaufin, day_ini, dt
41			use iniprint, only: prt_level
42			use com_io_dyn
43			use ener
44			use calfis_m, only: calfis
45			use exner_hyb_m, only: exner_hyb
46			use guide_m, only: guide
47			use pression_m, only: pression
48	guez	10	use pressure_var, only: p3d
49	guez	3
50			integer nq
51	guez	12	REAL, intent(in):: clesphy0(:)
52	guez	3
53	guez	10	! Variables dynamiques:
54	guez	3	REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
55			REAL teta(ip1jmp1, llm) ! temperature potentielle
56			REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields
57	guez	10	REAL ps(ip1jmp1) ! pression au sol, en Pa
58			REAL masse(ip1jmp1, llm) ! masse d'air
59			REAL phis(ip1jmp1) ! geopotentiel au sol
60
61			REAL time_0
62
63			! Variables local to the procedure:
64
65			! Variables dynamiques:
66
67	guez	3	REAL pks(ip1jmp1) ! exner au sol
68			REAL pk(ip1jmp1, llm) ! exner au milieu des couches
69			REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches
70			REAL phi(ip1jmp1, llm) ! geopotential
71			REAL w(ip1jmp1, llm) ! vitesse verticale
72
73			! variables dynamiques intermediaire pour le transport
74			REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse
75
76			! variables dynamiques au pas - 1
77			REAL vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)
78			REAL tetam1(ip1jmp1, llm), psm1(ip1jmp1)
79			REAL massem1(ip1jmp1, llm)
80
81			! tendances dynamiques
82			REAL dv(ip1jm, llm), du(ip1jmp1, llm)
83			REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)
84
85			! tendances de la dissipation
86			REAL dvdis(ip1jm, llm), dudis(ip1jmp1, llm)
87			REAL dtetadis(ip1jmp1, llm)
88
89			! tendances physiques
90			REAL dvfi(ip1jm, llm), dufi(ip1jmp1, llm)
91			REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)
92
93			! variables pour le fichier histoire
94
95			REAL tppn(iim), tpps(iim), tpn, tps
96
97			INTEGER itau, itaufinp1
98			INTEGER iday ! jour julien
99			REAL time ! Heure de la journee en fraction d'1 jour
100
101			REAL SSUM
102	guez	10	real finvmaold(ip1jmp1, llm)
103	guez	3
104			LOGICAL :: lafin=.false.
105			INTEGER ij, l
106
107			REAL rdayvrai, rdaym_ini
108	guez	10	LOGICAL:: callinigrads = .true.
109	guez	3
110			!+jld variables test conservation energie
111			REAL ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)
112			! Tendance de la temp. potentiel d (theta) / d t due a la
113			! tansformation d'energie cinetique en energie thermique
114			! cree par la dissipation
115			REAL dtetaecdt(ip1jmp1, llm)
116			REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm)
117			CHARACTER*15 ztit
118	guez	10	INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag.
119	guez	3
120	guez	10	logical:: dissip_conservative = .true.
121			LOGICAL:: prem = .true.
122	guez	12	logical forward, leapf, apphys, conser, apdiss
123	guez	3
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	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
142			CALL exner_hyb(ps, p3d, pks, pk, pkf)
143	guez	3
144			! Debut de l'integration temporelle:
145	guez	10	outer_loop:do
146	guez	3	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			conser = .FALSE.
171			apdiss = .FALSE.
172
173			IF (MOD(itau, iconser) == 0) conser = .TRUE.
174			IF (MOD(itau + 1, idissip) == 0) apdiss = .TRUE.
175			IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) apphys=.TRUE.
176
177			! calcul des tendances dynamiques:
178
179			CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
180
181			CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
182			conser, du, dv, dteta, dp, w, pbaru, pbarv, &
183			time + iday - day_ini)
184
185			! calcul des tendances advection des traceurs (dont l'humidite)
186
187			IF (forward .OR. leapf) THEN
188	guez	10	CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
189	guez	3	IF (offline) THEN
190			!maf stokage du flux de masse pour traceurs OFF-LINE
191			CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &
192			itau)
193			ENDIF
194			ENDIF
195
196			! integrations dynamique et traceurs:
197			CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
198	guez	12	dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &
199			finvmaold, leapf)
200	guez	3
201			! calcul des tendances physiques:
202
203			IF (apphys) THEN
204			IF (itau + 1 == itaufin) lafin = .TRUE.
205
206	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
207			CALL exner_hyb(ps, p3d, pks, pk, pkf)
208	guez	3
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	guez	10	CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
218			teta, q(:, :, 1), q(:, :, 2))
219	guez	3	ENDIF
220
221			CALL calfis(nq, lafin, rdayvrai, time, ucov, vcov, teta, q, &
222	guez	10	masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
223	guez	3	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	guez	10	CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, &
234	guez	3	teta, q(:, :, 1), q(:, :, 2))
235			ENDIF
236			ENDIF
237
238	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
239			CALL exner_hyb(ps, p3d, pks, pk, pkf)
240	guez	3
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	guez	10	CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
250	guez	3	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	guez	10	IF (itau == itaufinp1) exit outer_loop
309	guez	3
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	guez	5	CALL dynredem1("restart.nc", 0., vcov, ucov, teta, q, masse, ps)
323	guez	3	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	guez	10	end do outer_loop
346	guez	3
347			END SUBROUTINE leapfrog
348
349			end module leapfrog_m