Sources/dyn3d/leapfrog.f

module leapfrog_m

  IMPLICIT NONE

contains

  SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q)

    ! From dyn3d/leapfrog.F, version 1.6, 2005/04/13 08:58:34 revision 616
    ! Authors: P. Le Van, L. Fairhead, F. Hourdin

    ! Intégration temporelle du modèle : Matsuno-leapfrog scheme.

    use addfi_m, only: addfi
    use bilan_dyn_m, only: bilan_dyn
    use caladvtrac_m, only: caladvtrac
    use caldyn_m, only: caldyn
    USE calfis_m, ONLY: calfis
    USE comconst, ONLY: daysec, dtvr
    USE comgeom, ONLY: aire_2d, apoln, apols
    USE disvert_m, ONLY: ap, bp
    USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
         iflag_phys, iecri
    USE conf_guide_m, ONLY: ok_guide
    USE dimens_m, ONLY: iim, jjm, llm, nqmx
    use dissip_m, only: dissip
    USE dynetat0_m, ONLY: day_ini
    use dynredem1_m, only: dynredem1
    USE exner_hyb_m, ONLY: exner_hyb
    use filtreg_m, only: filtreg
    use fluxstokenc_m, only: fluxstokenc
    use geopot_m, only: geopot
    USE guide_m, ONLY: guide
    use inidissip_m, only: idissip
    use integrd_m, only: integrd
    use nr_util, only: assert
    USE pressure_var, ONLY: p3d
    USE temps, ONLY: itau_dyn
    use writedynav_m, only: writedynav
    use writehist_m, only: writehist

    ! Variables dynamiques:
    REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
    REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant

    REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
    ! potential temperature

    REAL, intent(inout):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol, en Pa
    REAL, intent(inout):: masse(:, :, :) ! (iim + 1, jjm + 1, llm) masse d'air
    REAL, intent(in):: phis(:, :) ! (iim + 1, jjm + 1) surface geopotential

    REAL, intent(inout):: q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nqmx)
    ! mass fractions of advected fields

    ! Local:

    ! Variables dynamiques:

    REAL pks(iim + 1, jjm + 1) ! exner au sol
    REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
    REAL pkf(iim + 1, jjm + 1, llm) ! exner filtr\'e au milieu des couches
    REAL phi(iim + 1, jjm + 1, llm) ! geopotential
    REAL w(iim + 1, jjm + 1, llm) ! vitesse verticale

    ! Variables dynamiques intermediaire pour le transport
    ! Flux de masse :
    REAL pbaru(iim + 1, jjm + 1, llm), pbarv(iim + 1, jjm, llm)

    ! Variables dynamiques au pas - 1
    REAL vcovm1(iim + 1, jjm, llm), ucovm1(iim + 1, jjm + 1, llm)
    REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
    REAL massem1(iim + 1, jjm + 1, llm)

    ! Tendances dynamiques
    REAL dv((iim + 1) * jjm, llm), dudyn(iim + 1, jjm + 1, llm)
    REAL dteta(iim + 1, jjm + 1, llm)
    real dp((iim + 1) * (jjm + 1))

    ! Tendances de la dissipation :
    REAL dvdis(iim + 1, jjm, llm), dudis(iim + 1, jjm + 1, llm)
    REAL dtetadis(iim + 1, jjm + 1, llm)

    ! Tendances physiques
    REAL dvfi(iim + 1, jjm, llm), dufi(iim + 1, jjm + 1, llm)
    REAL dtetafi(iim + 1, jjm + 1, llm), dqfi(iim + 1, jjm + 1, llm, nqmx)

    ! Variables pour le fichier histoire
    INTEGER itau ! index of the time step of the dynamics, starts at 0
    INTEGER itaufin
    REAL time ! time of day, as a fraction of day length
    real finvmaold(iim + 1, jjm + 1, llm)
    INTEGER l

    ! Variables test conservation \'energie
    REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)

    REAL vcont((iim + 1) * jjm, llm), ucont((iim + 1) * (jjm + 1), llm)
    logical leapf
    real dt ! time step, in s

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

    print *, "Call sequence information: leapfrog"
    call assert(shape(ucov) == (/iim + 1, jjm + 1, llm/), "leapfrog")

    itaufin = nday * day_step
    ! "day_step" is a multiple of "iperiod", therefore so is "itaufin".

    ! On initialise la pression et la fonction d'Exner :
    forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
    CALL exner_hyb(ps, p3d, pks, pk, pkf)

    time_integration: do itau = 0, itaufin - 1
       leapf = mod(itau, iperiod) /= 0
       if (leapf) then
          dt = 2 * dtvr
       else
          ! Matsuno
          dt = dtvr
          if (ok_guide) call guide(itau, ucov, vcov, teta, q(:, :, :, 1), ps)
          vcovm1 = vcov
          ucovm1 = ucov
          tetam1 = teta
          massem1 = masse
          psm1 = ps
          finvmaold = masse
          CALL filtreg(finvmaold, direct = .false., intensive = .false.)
       end if

       ! Calcul des tendances dynamiques:
       CALL geopot(teta, pk, pks, phis, phi)
       CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
            dudyn, dv, dteta, dp, w, pbaru, pbarv, &
            conser = MOD(itau, iconser) == 0)

       CALL caladvtrac(q, pbaru, pbarv, p3d, masse, teta, pk)

       ! Stokage du flux de masse pour traceurs offline:
       IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
            dtvr, itau)

       ! Int\'egrations dynamique et traceurs:
       CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, dteta, &
            dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, finvmaold, dt, &
            leapf)

       forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
       CALL exner_hyb(ps, p3d, pks, pk, pkf)

       if (.not. leapf) then
          ! Matsuno backward
          ! Calcul des tendances dynamiques:
          CALL geopot(teta, pk, pks, phis, phi)
          CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
               phi, dudyn, dv, dteta, dp, w, pbaru, pbarv, conser = .false.)

          ! integrations dynamique et traceurs:
          CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, &
               dteta, dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, &
               finvmaold, dtvr, leapf=.false.)

          forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
          CALL exner_hyb(ps, p3d, pks, pk, pkf)
       end if

       IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
          ! Calcul des tendances physiques :
          time = REAL(mod(itau, day_step)) / day_step
          IF (time > 1.) time = time - 1.
          CALL calfis(itau * dtvr / daysec + day_ini, time, ucov, vcov, teta, &
               q, pk, phis, phi, w, dufi, dvfi, dtetafi, dqfi, &
               lafin = itau + 1 == itaufin)

          CALL addfi(ucov, vcov, teta, q, dufi, dvfi, dtetafi, dqfi)
       ENDIF

       IF (MOD(itau + 1, idissip) == 0) THEN
          ! Dissipation horizontale et verticale des petites \'echelles

          ! calcul de l'\'energie cin\'etique 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 ajoute la tendance due \`a la transformation \'energie
          ! cin\'etique en \'energie thermique par la dissipation
          call covcont(llm, ucov, vcov, ucont, vcont)
          call enercin(vcov, ucov, vcont, ucont, ecin)
          dtetadis = dtetadis + (ecin0 - ecin) / pk
          teta = teta + dtetadis

          ! Calcul de la valeur moyenne aux p\^oles :
          forall (l = 1: llm)
             teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
                  / apoln
             teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
                  * teta(:iim, jjm + 1, l)) / apols
          END forall
       END IF

       IF (MOD(itau + 1, iperiod) == 0) THEN
          ! \'Ecriture du fichier histoire moyenne:
          CALL writedynav(vcov, ucov, teta, pk, phi, q, masse, ps, phis, &
               time = itau + 1)
          call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
               q(:, :, :, 1))
       ENDIF

       IF (MOD(itau + 1, iecri * day_step) == 0) THEN
          CALL geopot(teta, pk, pks, phis, phi)
          CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps)
       END IF
    end do time_integration

    CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
         itau = itau_dyn + itaufin)

    ! Calcul des tendances dynamiques:
    CALL geopot(teta, pk, pks, phis, phi)
    CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
         dudyn, dv, dteta, dp, w, pbaru, pbarv, &
         conser = MOD(itaufin, iconser) == 0)

  END SUBROUTINE leapfrog

end module leapfrog_m
1	guez	3	module leapfrog_m
2
3			IMPLICIT NONE
4
5			contains
6
7	guez	128	SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q)
8	guez	3
9	guez	67	! From dyn3d/leapfrog.F, version 1.6, 2005/04/13 08:58:34 revision 616
10	guez	27	! Authors: P. Le Van, L. Fairhead, F. Hourdin
11	guez	3
12	guez	129	! Intégration temporelle du modèle : Matsuno-leapfrog scheme.
13
14	guez	37	use addfi_m, only: addfi
15	guez	40	use bilan_dyn_m, only: bilan_dyn
16			use caladvtrac_m, only: caladvtrac
17	guez	43	use caldyn_m, only: caldyn
18	guez	26	USE calfis_m, ONLY: calfis
19	guez	103	USE comconst, ONLY: daysec, dtvr
20	guez	29	USE comgeom, ONLY: aire_2d, apoln, apols
21	guez	66	USE disvert_m, ONLY: ap, bp
22	guez	57	USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
23	guez	115	iflag_phys, iecri
24			USE conf_guide_m, ONLY: ok_guide
25	guez	29	USE dimens_m, ONLY: iim, jjm, llm, nqmx
26	guez	47	use dissip_m, only: dissip
27	guez	26	USE dynetat0_m, ONLY: day_ini
28	guez	27	use dynredem1_m, only: dynredem1
29	guez	26	USE exner_hyb_m, ONLY: exner_hyb
30	guez	27	use filtreg_m, only: filtreg
31	guez	69	use fluxstokenc_m, only: fluxstokenc
32	guez	43	use geopot_m, only: geopot
33	guez	26	USE guide_m, ONLY: guide
34			use inidissip_m, only: idissip
35	guez	32	use integrd_m, only: integrd
36	guez	55	use nr_util, only: assert
37	guez	26	USE pressure_var, ONLY: p3d
38	guez	28	USE temps, ONLY: itau_dyn
39	guez	56	use writedynav_m, only: writedynav
40	guez	69	use writehist_m, only: writehist
41	guez	3
42	guez	10	! Variables dynamiques:
43	guez	55	REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
44			REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
45	guez	43
46			REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
47			! potential temperature
48
49	guez	45	REAL, intent(inout):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol, en Pa
50	guez	70	REAL, intent(inout):: masse(:, :, :) ! (iim + 1, jjm + 1, llm) masse d'air
51	guez	69	REAL, intent(in):: phis(:, :) ! (iim + 1, jjm + 1) surface geopotential
52	guez	40
53			REAL, intent(inout):: q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nqmx)
54			! mass fractions of advected fields
55
56	guez	91	! Local:
57	guez	10
58			! Variables dynamiques:
59
60	guez	70	REAL pks(iim + 1, jjm + 1) ! exner au sol
61	guez	29	REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
62	guez	90	REAL pkf(iim + 1, jjm + 1, llm) ! exner filtr\'e au milieu des couches
63	guez	47	REAL phi(iim + 1, jjm + 1, llm) ! geopotential
64	guez	91	REAL w(iim + 1, jjm + 1, llm) ! vitesse verticale
65	guez	3
66	guez	55	! Variables dynamiques intermediaire pour le transport
67			! Flux de masse :
68	guez	91	REAL pbaru(iim + 1, jjm + 1, llm), pbarv(iim + 1, jjm, llm)
69	guez	3
70	guez	56	! Variables dynamiques au pas - 1
71	guez	55	REAL vcovm1(iim + 1, jjm, llm), ucovm1(iim + 1, jjm + 1, llm)
72	guez	29	REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
73	guez	67	REAL massem1(iim + 1, jjm + 1, llm)
74	guez	3
75	guez	56	! Tendances dynamiques
76	guez	88	REAL dv((iim + 1) * jjm, llm), dudyn(iim + 1, jjm + 1, llm)
77	guez	71	REAL dteta(iim + 1, jjm + 1, llm)
78	guez	55	real dp((iim + 1) * (jjm + 1))
79	guez	3
80	guez	56	! Tendances de la dissipation :
81	guez	55	REAL dvdis(iim + 1, jjm, llm), dudis(iim + 1, jjm + 1, llm)
82	guez	29	REAL dtetadis(iim + 1, jjm + 1, llm)
83	guez	3
84	guez	56	! Tendances physiques
85	guez	91	REAL dvfi(iim + 1, jjm, llm), dufi(iim + 1, jjm + 1, llm)
86			REAL dtetafi(iim + 1, jjm + 1, llm), dqfi(iim + 1, jjm + 1, llm, nqmx)
87	guez	3
88	guez	56	! Variables pour le fichier histoire
89	guez	22	INTEGER itau ! index of the time step of the dynamics, starts at 0
90	guez	27	INTEGER itaufin
91	guez	20	REAL time ! time of day, as a fraction of day length
92	guez	67	real finvmaold(iim + 1, jjm + 1, llm)
93	guez	33	INTEGER l
94	guez	3
95	guez	90	! Variables test conservation \'energie
96	guez	29	REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
97	guez	43
98	guez	55	REAL vcont((iim + 1) * jjm, llm), ucont((iim + 1) * (jjm + 1), llm)
99	guez	33	logical leapf
100	guez	91	real dt ! time step, in s
101	guez	3
102			!---------------------------------------------------
103
104			print *, "Call sequence information: leapfrog"
105	guez	55	call assert(shape(ucov) == (/iim + 1, jjm + 1, llm/), "leapfrog")
106	guez	3
107			itaufin = nday * day_step
108	guez	62	! "day_step" is a multiple of "iperiod", therefore so is "itaufin".
109	guez	30
110	guez	3	! On initialise la pression et la fonction d'Exner :
111	guez	37	forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
112	guez	10	CALL exner_hyb(ps, p3d, pks, pk, pkf)
113	guez	3
114	guez	40	time_integration: do itau = 0, itaufin - 1
115	guez	33	leapf = mod(itau, iperiod) /= 0
116			if (leapf) then
117			dt = 2 * dtvr
118			else
119			! Matsuno
120			dt = dtvr
121	guez	108	if (ok_guide) call guide(itau, ucov, vcov, teta, q(:, :, :, 1), ps)
122	guez	33	vcovm1 = vcov
123			ucovm1 = ucov
124			tetam1 = teta
125			massem1 = masse
126			psm1 = ps
127			finvmaold = masse
128	guez	107	CALL filtreg(finvmaold, direct = .false., intensive = .false.)
129	guez	33	end if
130	guez	30
131			! Calcul des tendances dynamiques:
132	guez	70	CALL geopot(teta, pk, pks, phis, phi)
133	guez	30	CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
134	guez	128	dudyn, dv, dteta, dp, w, pbaru, pbarv, &
135	guez	78	conser = MOD(itau, iconser) == 0)
136	guez	30
137	guez	71	CALL caladvtrac(q, pbaru, pbarv, p3d, masse, teta, pk)
138	guez	33
139	guez	30	! Stokage du flux de masse pour traceurs offline:
140			IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
141			dtvr, itau)
142
143	guez	90	! Int\'egrations dynamique et traceurs:
144	guez	47	CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, dteta, &
145			dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, finvmaold, dt, &
146			leapf)
147	guez	30
148	guez	97	forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
149			CALL exner_hyb(ps, p3d, pks, pk, pkf)
150
151	guez	33	if (.not. leapf) then
152			! Matsuno backward
153	guez	27	! Calcul des tendances dynamiques:
154	guez	70	CALL geopot(teta, pk, pks, phis, phi)
155	guez	33	CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
156	guez	128	phi, dudyn, dv, dteta, dp, w, pbaru, pbarv, conser = .false.)
157	guez	3
158			! integrations dynamique et traceurs:
159	guez	47	CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, &
160			dteta, dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, &
161			finvmaold, dtvr, leapf=.false.)
162	guez	97
163			forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
164			CALL exner_hyb(ps, p3d, pks, pk, pkf)
165	guez	33	end if
166	guez	3
167	guez	33	IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
168	guez	128	! Calcul des tendances physiques :
169			time = REAL(mod(itau, day_step)) / day_step
170	guez	33	IF (time > 1.) time = time - 1.
171	guez	128	CALL calfis(itau * dtvr / daysec + day_ini, time, ucov, vcov, teta, &
172			q, pk, phis, phi, w, dufi, dvfi, dtetafi, dqfi, &
173			lafin = itau + 1 == itaufin)
174	guez	3
175	guez	91	CALL addfi(ucov, vcov, teta, q, dufi, dvfi, dtetafi, dqfi)
176	guez	33	ENDIF
177	guez	3
178	guez	33	IF (MOD(itau + 1, idissip) == 0) THEN
179	guez	90	! Dissipation horizontale et verticale des petites \'echelles
180	guez	3
181	guez	90	! calcul de l'\'energie cin\'etique avant dissipation
182	guez	33	call covcont(llm, ucov, vcov, ucont, vcont)
183			call enercin(vcov, ucov, vcont, ucont, ecin0)
184	guez	3
185	guez	33	! dissipation
186			CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
187	guez	55	ucov = ucov + dudis
188			vcov = vcov + dvdis
189	guez	3
190	guez	90	! On ajoute la tendance due \`a la transformation \'energie
191			! cin\'etique en \'energie thermique par la dissipation
192	guez	33	call covcont(llm, ucov, vcov, ucont, vcont)
193			call enercin(vcov, ucov, vcont, ucont, ecin)
194	guez	56	dtetadis = dtetadis + (ecin0 - ecin) / pk
195	guez	55	teta = teta + dtetadis
196	guez	3
197	guez	90	! Calcul de la valeur moyenne aux p\^oles :
198	guez	33	forall (l = 1: llm)
199			teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
200			/ apoln
201			teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
202			* teta(:iim, jjm + 1, l)) / apols
203			END forall
204			END IF
205	guez	3
206	guez	33	IF (MOD(itau + 1, iperiod) == 0) THEN
207	guez	90	! \'Ecriture du fichier histoire moyenne:
208	guez	61	CALL writedynav(vcov, ucov, teta, pk, phi, q, masse, ps, phis, &
209			time = itau + 1)
210	guez	40	call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
211	guez	57	q(:, :, :, 1))
212	guez	33	ENDIF
213	guez	68
214			IF (MOD(itau + 1, iecri * day_step) == 0) THEN
215	guez	70	CALL geopot(teta, pk, pks, phis, phi)
216	guez	69	CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps)
217	guez	68	END IF
218	guez	40	end do time_integration
219	guez	3
220	guez	30	CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
221	guez	67	itau = itau_dyn + itaufin)
222	guez	30
223			! Calcul des tendances dynamiques:
224	guez	70	CALL geopot(teta, pk, pks, phis, phi)
225	guez	30	CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
226	guez	128	dudyn, dv, dteta, dp, w, pbaru, pbarv, &
227	guez	67	conser = MOD(itaufin, iconser) == 0)
228	guez	56
229	guez	3	END SUBROUTINE leapfrog
230
231			end module leapfrog_m