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

    REAL, intent(in):: time_0

    ! 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
    REAL rdayvrai, rdaym_ini

    ! 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, time_0, &
            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, time_0, &
               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:

          rdaym_ini = itau * dtvr / daysec
          rdayvrai = rdaym_ini + day_ini
          time = REAL(mod(itau, day_step)) / day_step + time_0
          IF (time > 1.) time = time - 1.

          CALL calfis(rdayvrai, time, ucov, vcov, teta, q, pk, phis, phi, w, &
               dufi, dvfi, dtetafi, dqfi, lafin = itau + 1 == itaufin)

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