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

    USE calfis_m, ONLY: calfis
    USE com_io_dyn, ONLY: histaveid
    USE comconst, ONLY: daysec, dtphys, dtvr
    USE comgeom, ONLY: aire, apoln, apols
    USE comvert, ONLY: ap, bp
    USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
         periodav
    USE dimens_m, ONLY: iim, llm, nqmx
    USE dynetat0_m, ONLY: day_ini
    use dynredem1_m, only: dynredem1
    USE exner_hyb_m, ONLY: exner_hyb
    use filtreg_m, only: filtreg
    USE guide_m, ONLY: guide
    use inidissip_m, only: idissip
    USE logic, ONLY: iflag_phys, ok_guide
    USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1, jjp1
    USE pression_m, ONLY: pression
    USE pressure_var, ONLY: p3d
    USE temps, ONLY: dt, itau_dyn

    ! 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 itaufin
    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)
    logical forward, leapf

    !---------------------------------------------------

    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)

    ! Début 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
          ! Calcul des tendances dynamiques:
          CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
          CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
               MOD(itau, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
               time + iday - day_ini)

          IF (forward .OR. leapf) THEN
             ! Calcul des tendances advection des traceurs (dont l'humidité)
             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 (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) 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

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

          CALL pression(ip1jmp1, ap, bp, ps, p3d)
          CALL exner_hyb(ps, p3d, pks, pk, pkf)

          IF (MOD(itau + 1, idissip) == 0) 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

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