trunk/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: 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
    REAL dt

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

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

          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 à la transformation Ec -> E
             ! thermique créée 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 pôles
             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	28	USE temps, ONLY: 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	28	REAL dt
95	guez	3
96			!---------------------------------------------------
97
98			print *, "Call sequence information: leapfrog"
99
100			itaufin = nday * day_step
101			itau = 0
102			iday = day_ini
103			time = time_0
104	guez	24	dq = 0.
105	guez	3	! On initialise la pression et la fonction d'Exner :
106	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
107			CALL exner_hyb(ps, p3d, pks, pk, pkf)
108	guez	3
109	guez	27	! Début de l'integration temporelle :
110	guez	10	outer_loop:do
111	guez	25	if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
112			call guide(itau, ucov, vcov, teta, q, masse, ps)
113			vcovm1 = vcov
114			ucovm1 = ucov
115			tetam1 = teta
116			massem1 = masse
117			psm1 = ps
118	guez	3	forward = .TRUE.
119			leapf = .FALSE.
120			dt = dtvr
121	guez	25	finvmaold = masse
122	guez	3	CALL filtreg(finvmaold, jjp1, llm, - 2, 2, .TRUE., 1)
123
124			do
125	guez	27	! Calcul des tendances dynamiques:
126	guez	3	CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
127			CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
128	guez	27	MOD(itau, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
129	guez	3	time + iday - day_ini)
130
131			IF (forward .OR. leapf) THEN
132	guez	27	! Calcul des tendances advection des traceurs (dont l'humidité)
133	guez	10	CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
134	guez	3	IF (offline) THEN
135	guez	25	! Stokage du flux de masse pour traceurs off-line
136	guez	3	CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, &
137			itau)
138			ENDIF
139			ENDIF
140
141			! integrations dynamique et traceurs:
142			CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
143	guez	12	dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, &
144	guez	28	finvmaold, leapf, dt)
145	guez	3
146	guez	27	IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
147	guez	25	! calcul des tendances physiques:
148	guez	3	IF (itau + 1 == itaufin) lafin = .TRUE.
149
150	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
151			CALL exner_hyb(ps, p3d, pks, pk, pkf)
152	guez	3
153			rdaym_ini = itau * dtvr / daysec
154			rdayvrai = rdaym_ini + day_ini
155
156	guez	23	CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &
157	guez	10	masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
158	guez	13	dufi, dvfi, dtetafi, dqfi, dpfi)
159	guez	3
160			! ajout des tendances physiques:
161			CALL addfi(nqmx, dtphys, &
162			ucov, vcov, teta, q, ps, &
163			dufi, dvfi, dtetafi, dqfi, dpfi)
164			ENDIF
165
166	guez	10	CALL pression(ip1jmp1, ap, bp, ps, p3d)
167			CALL exner_hyb(ps, p3d, pks, pk, pkf)
168	guez	3
169	guez	27	IF (MOD(itau + 1, idissip) == 0) THEN
170	guez	25	! dissipation horizontale et verticale des petites echelles:
171	guez	3
172			! calcul de l'energie cinetique avant dissipation
173			call covcont(llm, ucov, vcov, ucont, vcont)
174			call enercin(vcov, ucov, vcont, ucont, ecin0)
175
176			! dissipation
177	guez	10	CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
178	guez	3	ucov=ucov + dudis
179			vcov=vcov + dvdis
180
181	guez	28	! On rajoute la tendance due à la transformation Ec -> E
182			! thermique créée lors de la dissipation
183	guez	27	call covcont(llm, ucov, vcov, ucont, vcont)
184			call enercin(vcov, ucov, vcont, ucont, ecin)
185			dtetaecdt= (ecin0 - ecin) / pk
186			dtetadis=dtetadis + dtetaecdt
187	guez	3	teta=teta + dtetadis
188
189	guez	28	! Calcul de la valeur moyenne unique de h aux pôles
190	guez	3	DO l = 1, llm
191			DO ij = 1, iim
192			tppn(ij) = aire(ij) * teta(ij, l)
193			tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)
194			ENDDO
195	guez	25	tpn = SUM(tppn) / apoln
196			tps = SUM(tpps) / apols
197	guez	3
198			DO ij = 1, iip1
199			teta(ij, l) = tpn
200			teta(ij + ip1jm, l) = tps
201			ENDDO
202			ENDDO
203
204			DO ij = 1, iim
205			tppn(ij) = aire(ij) * ps(ij)
206			tpps(ij) = aire(ij + ip1jm) * ps(ij + ip1jm)
207			ENDDO
208	guez	25	tpn = SUM(tppn) / apoln
209			tps = SUM(tpps) / apols
210	guez	3
211			DO ij = 1, iip1
212			ps(ij) = tpn
213			ps(ij + ip1jm) = tps
214			ENDDO
215			END IF
216
217			! fin de l'intégration dynamique et physique pour le pas "itau"
218			! préparation du pas d'intégration suivant
219
220			! schema matsuno + leapfrog
221			IF (forward .OR. leapf) THEN
222			itau = itau + 1
223			iday = day_ini + itau / day_step
224			time = REAL(itau - (iday - day_ini) * day_step) / day_step &
225			+ time_0
226			IF (time > 1.) THEN
227			time = time - 1.
228			iday = iday + 1
229			ENDIF
230			ENDIF
231
232	guez	24	IF (itau == itaufin + 1) exit outer_loop
233	guez	3
234			IF (MOD(itau, iperiod) == 0 .OR. itau == itaufin) THEN
235	guez	25	! ecriture du fichier histoire moyenne:
236	guez	3	CALL writedynav(histaveid, nqmx, itau, vcov, &
237			ucov, teta, pk, phi, q, masse, ps, phis)
238			call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
239			ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
240			ENDIF
241
242			IF (itau == itaufin) THEN
243	guez	27	CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
244			itau=itau_dyn+itaufin)
245	guez	3	ENDIF
246
247			! gestion de l'integration temporelle:
248			IF (MOD(itau, iperiod) == 0) exit
249			IF (MOD(itau - 1, iperiod) == 0) THEN
250			IF (forward) THEN
251			! fin du pas forward et debut du pas backward
252			forward = .FALSE.
253			leapf = .FALSE.
254			ELSE
255			! fin du pas backward et debut du premier pas leapfrog
256			leapf = .TRUE.
257			dt = 2. * dtvr
258			END IF
259			ELSE
260	guez	20	! pas leapfrog
261	guez	3	leapf = .TRUE.
262			dt = 2. * dtvr
263			END IF
264			end do
265	guez	10	end do outer_loop
266	guez	3
267			END SUBROUTINE leapfrog
268
269			end module leapfrog_m