trunk/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: dtvr
    USE comgeom, ONLY: aire_2d, apoln, apols
    use covcont_m, only: covcont
    USE disvert_m, ONLY: ap, bp
    USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, &
         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 enercin_m, only: enercin
    USE exner_hyb_m, ONLY: exner_hyb
    use filtreg_scal_m, only: filtreg_scal
    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 temps, ONLY: itau_dyn
    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 interm\'ediaires 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), du(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
    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

    REAL p3d(iim + 1, jjm + 1, llm + 1) ! pressure at layer interfaces, in Pa
    ! ("p3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
    ! for interface "l")

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

    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 = pk
    CALL filtreg_scal(pkf, direct = .true., intensive = .true.)

    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
       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, &
            du, dv, dteta, dp, w, pbaru, pbarv, &
            conser = MOD(itau, iconser) == 0)

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

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

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

       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, du, dv, dteta, dp, w, pbaru, pbarv, conser = .false.)

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

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

       IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys) THEN
          CALL calfis(ucov, vcov, teta, q, p3d, pk, phis, phi, w, dufi, dvfi, &
               dtetafi, dqfi, dayvrai = itau / day_step + day_ini, &
               time = REAL(mod(itau, day_step)) / day_step, &
               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
          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(vcov, ucov, teta, pk, phi, q, masse, ps, itau)
       END IF
    end do time_integration

    CALL dynredem1(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, &
         du, 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	178	USE comconst, ONLY: dtvr
20	guez	29	USE comgeom, ONLY: aire_2d, apoln, apols
21	guez	161	use covcont_m, only: covcont
22	guez	66	USE disvert_m, ONLY: ap, bp
23	guez	224	USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, &
24	guez	115	iflag_phys, iecri
25			USE conf_guide_m, ONLY: ok_guide
26	guez	29	USE dimens_m, ONLY: iim, jjm, llm, nqmx
27	guez	47	use dissip_m, only: dissip
28	guez	26	USE dynetat0_m, ONLY: day_ini
29	guez	27	use dynredem1_m, only: dynredem1
30	guez	207	use enercin_m, only: enercin
31	guez	26	USE exner_hyb_m, ONLY: exner_hyb
32	guez	137	use filtreg_scal_m, only: filtreg_scal
33	guez	43	use geopot_m, only: geopot
34	guez	26	USE guide_m, ONLY: guide
35			use inidissip_m, only: idissip
36	guez	32	use integrd_m, only: integrd
37	guez	55	use nr_util, only: assert
38	guez	28	USE temps, ONLY: itau_dyn
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	91	! Local:
56	guez	10
57			! Variables dynamiques:
58
59	guez	70	REAL pks(iim + 1, jjm + 1) ! exner au sol
60	guez	29	REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
61	guez	90	REAL pkf(iim + 1, jjm + 1, llm) ! exner filtr\'e au milieu des couches
62	guez	47	REAL phi(iim + 1, jjm + 1, llm) ! geopotential
63	guez	91	REAL w(iim + 1, jjm + 1, llm) ! vitesse verticale
64	guez	3
65	guez	260	! Variables dynamiques interm\'ediaires pour le transport
66	guez	55	! Flux de masse :
67	guez	91	REAL pbaru(iim + 1, jjm + 1, llm), pbarv(iim + 1, jjm, llm)
68	guez	3
69	guez	56	! Variables dynamiques au pas - 1
70	guez	55	REAL vcovm1(iim + 1, jjm, llm), ucovm1(iim + 1, jjm + 1, llm)
71	guez	29	REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
72	guez	67	REAL massem1(iim + 1, jjm + 1, llm)
73	guez	3
74	guez	56	! Tendances dynamiques
75	guez	260	REAL dv((iim + 1) * jjm, llm), du(iim + 1, jjm + 1, llm)
76	guez	71	REAL dteta(iim + 1, jjm + 1, llm)
77	guez	252	real dp(iim + 1, jjm + 1)
78	guez	3
79	guez	56	! Tendances de la dissipation :
80	guez	55	REAL dvdis(iim + 1, jjm, llm), dudis(iim + 1, jjm + 1, llm)
81	guez	29	REAL dtetadis(iim + 1, jjm + 1, llm)
82	guez	3
83	guez	56	! Tendances physiques
84	guez	91	REAL dvfi(iim + 1, jjm, llm), dufi(iim + 1, jjm + 1, llm)
85			REAL dtetafi(iim + 1, jjm + 1, llm), dqfi(iim + 1, jjm + 1, llm, nqmx)
86	guez	3
87	guez	56	! Variables pour le fichier histoire
88	guez	22	INTEGER itau ! index of the time step of the dynamics, starts at 0
89	guez	27	INTEGER itaufin
90	guez	33	INTEGER l
91	guez	3
92	guez	90	! Variables test conservation \'energie
93	guez	29	REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
94	guez	43
95	guez	55	REAL vcont((iim + 1) * jjm, llm), ucont((iim + 1) * (jjm + 1), llm)
96	guez	33	logical leapf
97	guez	91	real dt ! time step, in s
98	guez	3
99	guez	212	REAL p3d(iim + 1, jjm + 1, llm + 1) ! pressure at layer interfaces, in Pa
100	guez	162	! ("p3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
101			! for interface "l")
102
103	guez	3	!---------------------------------------------------
104
105			print *, "Call sequence information: leapfrog"
106	guez	55	call assert(shape(ucov) == (/iim + 1, jjm + 1, llm/), "leapfrog")
107	guez	3
108			itaufin = nday * day_step
109	guez	62	! "day_step" is a multiple of "iperiod", therefore so is "itaufin".
110	guez	30
111	guez	3	! On initialise la pression et la fonction d'Exner :
112	guez	37	forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
113	guez	137	CALL exner_hyb(ps, p3d, pks, pk)
114			pkf = pk
115			CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
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			end if
131	guez	30
132			! Calcul des tendances dynamiques:
133	guez	70	CALL geopot(teta, pk, pks, phis, phi)
134	guez	30	CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
135	guez	260	du, dv, dteta, dp, w, pbaru, pbarv, &
136	guez	78	conser = MOD(itau, iconser) == 0)
137	guez	30
138	guez	71	CALL caladvtrac(q, pbaru, pbarv, p3d, masse, teta, pk)
139	guez	33
140	guez	90	! Int\'egrations dynamique et traceurs:
141	guez	260	CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &
142	guez	161	dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, dt, leapf)
143	guez	30
144	guez	97	forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
145	guez	137	CALL exner_hyb(ps, p3d, pks, pk)
146			pkf = pk
147			CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
148	guez	97
149	guez	33	if (.not. leapf) then
150			! Matsuno backward
151	guez	27	! Calcul des tendances dynamiques:
152	guez	70	CALL geopot(teta, pk, pks, phis, phi)
153	guez	33	CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
154	guez	260	phi, du, dv, dteta, dp, w, pbaru, pbarv, conser = .false.)
155	guez	3
156			! integrations dynamique et traceurs:
157	guez	260	CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
158	guez	161	dteta, dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, dtvr, &
159			leapf=.false.)
160	guez	97
161			forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
162	guez	137	CALL exner_hyb(ps, p3d, pks, pk)
163			pkf = pk
164			CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
165	guez	33	end if
166	guez	3
167	guez	208	IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys) THEN
168	guez	162	CALL calfis(ucov, vcov, teta, q, p3d, pk, phis, phi, w, dufi, dvfi, &
169	guez	154	dtetafi, dqfi, dayvrai = itau / day_step + day_ini, &
170			time = REAL(mod(itau, day_step)) / day_step, &
171	guez	128	lafin = itau + 1 == itaufin)
172	guez	3
173	guez	91	CALL addfi(ucov, vcov, teta, q, dufi, dvfi, dtetafi, dqfi)
174	guez	33	ENDIF
175	guez	3
176	guez	33	IF (MOD(itau + 1, idissip) == 0) THEN
177	guez	90	! Dissipation horizontale et verticale des petites \'echelles
178	guez	3
179	guez	90	! calcul de l'\'energie cin\'etique avant dissipation
180	guez	33	call covcont(llm, ucov, vcov, ucont, vcont)
181			call enercin(vcov, ucov, vcont, ucont, ecin0)
182	guez	3
183	guez	33	! dissipation
184			CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
185	guez	55	ucov = ucov + dudis
186			vcov = vcov + dvdis
187	guez	3
188	guez	90	! On ajoute la tendance due \`a la transformation \'energie
189			! cin\'etique en \'energie thermique par la dissipation
190	guez	33	call covcont(llm, ucov, vcov, ucont, vcont)
191			call enercin(vcov, ucov, vcont, ucont, ecin)
192	guez	56	dtetadis = dtetadis + (ecin0 - ecin) / pk
193	guez	55	teta = teta + dtetadis
194	guez	3
195	guez	90	! Calcul de la valeur moyenne aux p\^oles :
196	guez	33	forall (l = 1: llm)
197			teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
198			/ apoln
199	guez	212	teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm + 1) &
200	guez	33	* teta(:iim, jjm + 1, l)) / apols
201			END forall
202			END IF
203	guez	3
204	guez	33	IF (MOD(itau + 1, iperiod) == 0) THEN
205	guez	40	call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
206	guez	57	q(:, :, :, 1))
207	guez	33	ENDIF
208	guez	68
209			IF (MOD(itau + 1, iecri * day_step) == 0) THEN
210	guez	70	CALL geopot(teta, pk, pks, phis, phi)
211	guez	261	CALL writehist(vcov, ucov, teta, pk, phi, q, masse, ps, itau)
212	guez	68	END IF
213	guez	40	end do time_integration
214	guez	3
215	guez	157	CALL dynredem1(vcov, ucov, teta, q, masse, ps, itau = itau_dyn + itaufin)
216	guez	30
217			! Calcul des tendances dynamiques:
218	guez	70	CALL geopot(teta, pk, pks, phis, phi)
219	guez	30	CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
220	guez	260	du, dv, dteta, dp, w, pbaru, pbarv, &
221	guez	67	conser = MOD(itaufin, iconser) == 0)
222	guez	56
223	guez	3	END SUBROUTINE leapfrog
224
225			end module leapfrog_m