/[lmdze]/trunk/dyn3d/leapfrog.f
ViewVC logotype

Annotation of /trunk/dyn3d/leapfrog.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 64 - (hide annotations)
Wed Aug 29 14:47:17 2012 UTC (11 years, 8 months ago) by guez
Original Path: trunk/libf/dyn3d/leapfrog.f90
File size: 9024 byte(s)
Removed variable lstardis in module comdissnew and procedures gradiv
and nxgrarot. lstardir had to be true. gradiv and nxgrarot were called
if lstardis was false. Removed argument iter of procedure
filtreg. iter had to be 1. gradiv and nxgrarot called filtreg with
iter == 2.

Moved procedure flxsetup into module yoecumf. Module yoecumf is only
used in program units of directory Conflx, moved it there.

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 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     USE comconst, ONLY: daysec, dtphys, dtvr
19 guez 29 USE comgeom, ONLY: aire_2d, apoln, apols
20 guez 26 USE comvert, ONLY: ap, bp
21 guez 57 USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
22     iflag_phys, ok_guide
23 guez 29 USE dimens_m, ONLY: iim, jjm, llm, nqmx
24 guez 47 use dissip_m, only: dissip
25 guez 26 USE dynetat0_m, ONLY: day_ini
26 guez 27 use dynredem1_m, only: dynredem1
27 guez 26 USE exner_hyb_m, ONLY: exner_hyb
28 guez 27 use filtreg_m, only: filtreg
29 guez 43 use geopot_m, only: geopot
30 guez 26 USE guide_m, ONLY: guide
31     use inidissip_m, only: idissip
32 guez 32 use integrd_m, only: integrd
33 guez 55 use nr_util, only: assert
34 guez 26 USE pressure_var, ONLY: p3d
35 guez 28 USE temps, ONLY: itau_dyn
36 guez 56 use writedynav_m, only: writedynav
37 guez 3
38 guez 10 ! Variables dynamiques:
39 guez 55 REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
40     REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
41 guez 43
42     REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
43     ! potential temperature
44    
45 guez 45 REAL, intent(inout):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol, en Pa
46 guez 55 REAL masse((iim + 1) * (jjm + 1), llm) ! masse d'air
47     REAL phis((iim + 1) * (jjm + 1)) ! geopotentiel au sol
48 guez 40
49     REAL, intent(inout):: q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nqmx)
50     ! mass fractions of advected fields
51    
52 guez 24 REAL, intent(in):: time_0
53 guez 10
54     ! Variables local to the procedure:
55    
56     ! Variables dynamiques:
57    
58 guez 55 REAL pks((iim + 1) * (jjm + 1)) ! exner au sol
59 guez 29 REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
60 guez 55 REAL pkf((iim + 1) * (jjm + 1), llm) ! exner filt.au milieu des couches
61 guez 47 REAL phi(iim + 1, jjm + 1, llm) ! geopotential
62 guez 55 REAL w((iim + 1) * (jjm + 1), llm) ! vitesse verticale
63 guez 3
64 guez 55 ! Variables dynamiques intermediaire pour le transport
65     ! Flux de masse :
66     REAL pbaru((iim + 1) * (jjm + 1), llm), pbarv((iim + 1) * jjm, llm)
67 guez 3
68 guez 56 ! Variables dynamiques au pas - 1
69 guez 55 REAL vcovm1(iim + 1, jjm, llm), ucovm1(iim + 1, jjm + 1, llm)
70 guez 29 REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
71 guez 55 REAL massem1((iim + 1) * (jjm + 1), llm)
72 guez 3
73 guez 56 ! Tendances dynamiques
74 guez 55 REAL dv((iim + 1) * jjm, llm), dudyn((iim + 1) * (jjm + 1), llm)
75     REAL dteta(iim + 1, jjm + 1, llm), dq((iim + 1) * (jjm + 1), llm, nqmx)
76     real dp((iim + 1) * (jjm + 1))
77 guez 3
78 guez 56 ! Tendances de la dissipation :
79 guez 55 REAL dvdis(iim + 1, jjm, llm), dudis(iim + 1, jjm + 1, llm)
80 guez 29 REAL dtetadis(iim + 1, jjm + 1, llm)
81 guez 3
82 guez 56 ! Tendances physiques
83 guez 55 REAL dvfi((iim + 1) * jjm, llm), dufi((iim + 1) * (jjm + 1), llm)
84     REAL dtetafi(iim + 1, jjm + 1, llm), dqfi((iim + 1) * (jjm + 1), llm, nqmx)
85     real dpfi((iim + 1) * (jjm + 1))
86 guez 3
87 guez 56 ! Variables pour le fichier histoire
88 guez 3
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 55 real finvmaold((iim + 1) * (jjm + 1), llm)
93 guez 33 INTEGER l
94 guez 3 REAL rdayvrai, rdaym_ini
95    
96 guez 24 ! Variables test conservation energie
97 guez 29 REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
98 guez 43
99 guez 55 REAL vcont((iim + 1) * jjm, llm), ucont((iim + 1) * (jjm + 1), llm)
100 guez 33 logical leapf
101     real dt
102 guez 3
103     !---------------------------------------------------
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 24 dq = 0.
112 guez 33
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     if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
125     call guide(itau, ucov, vcov, teta, q, masse, ps)
126     vcovm1 = vcov
127     ucovm1 = ucov
128     tetam1 = teta
129     massem1 = masse
130     psm1 = ps
131     finvmaold = masse
132 guez 64 CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE.)
133 guez 33 end if
134 guez 30
135     ! Calcul des tendances dynamiques:
136 guez 55 CALL geopot((iim + 1) * (jjm + 1), teta, pk, pks, phis, phi)
137 guez 30 CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
138 guez 47 dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
139     conser=MOD(itau, iconser)==0)
140 guez 30
141     ! Calcul des tendances advection des traceurs (dont l'humidité)
142     CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
143 guez 33
144 guez 30 ! Stokage du flux de masse pour traceurs offline:
145     IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
146     dtvr, itau)
147    
148 guez 46 ! Integrations dynamique et traceurs:
149 guez 47 CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, dteta, &
150     dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, finvmaold, dt, &
151     leapf)
152 guez 30
153 guez 33 if (.not. leapf) then
154     ! Matsuno backward
155 guez 37 forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
156 guez 33 CALL exner_hyb(ps, p3d, pks, pk, pkf)
157 guez 30
158 guez 27 ! Calcul des tendances dynamiques:
159 guez 55 CALL geopot((iim + 1) * (jjm + 1), teta, pk, pks, phis, phi)
160 guez 33 CALL caldyn(itau + 1, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
161 guez 47 phi, dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
162     conser=.false.)
163 guez 3
164     ! integrations dynamique et traceurs:
165 guez 47 CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, &
166     dteta, dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, &
167     finvmaold, dtvr, leapf=.false.)
168 guez 33 end if
169 guez 3
170 guez 33 IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
171     ! calcul des tendances physiques:
172 guez 3
173 guez 37 forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
174 guez 33 CALL exner_hyb(ps, p3d, pks, pk, pkf)
175 guez 3
176 guez 33 rdaym_ini = itau * dtvr / daysec
177     rdayvrai = rdaym_ini + day_ini
178     time = REAL(mod(itau, day_step)) / day_step + time_0
179     IF (time > 1.) time = time - 1.
180 guez 3
181 guez 37 CALL calfis(rdayvrai, time, ucov, vcov, teta, q, masse, ps, pk, &
182 guez 47 phis, phi, dudyn, dv, dq, w, dufi, dvfi, dtetafi, dqfi, dpfi, &
183 guez 44 lafin=itau+1==itaufin)
184 guez 3
185 guez 33 ! ajout des tendances physiques:
186     CALL addfi(nqmx, dtphys, ucov, vcov, teta, q, ps, dufi, dvfi, &
187     dtetafi, dqfi, dpfi)
188     ENDIF
189 guez 3
190 guez 37 forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
191 guez 33 CALL exner_hyb(ps, p3d, pks, pk, pkf)
192 guez 3
193 guez 33 IF (MOD(itau + 1, idissip) == 0) THEN
194 guez 55 ! Dissipation horizontale et verticale des petites échelles
195 guez 3
196 guez 55 ! calcul de l'énergie cinétique avant dissipation
197 guez 33 call covcont(llm, ucov, vcov, ucont, vcont)
198     call enercin(vcov, ucov, vcont, ucont, ecin0)
199 guez 3
200 guez 33 ! dissipation
201     CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
202 guez 55 ucov = ucov + dudis
203     vcov = vcov + dvdis
204 guez 3
205 guez 55 ! On ajoute la tendance due à la transformation énergie
206     ! cinétique en énergie thermique par la dissipation
207 guez 33 call covcont(llm, ucov, vcov, ucont, vcont)
208     call enercin(vcov, ucov, vcont, ucont, ecin)
209 guez 56 dtetadis = dtetadis + (ecin0 - ecin) / pk
210 guez 55 teta = teta + dtetadis
211 guez 3
212 guez 33 ! Calcul de la valeur moyenne aux pôles :
213     forall (l = 1: llm)
214     teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
215     / apoln
216     teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
217     * teta(:iim, jjm + 1, l)) / apols
218     END forall
219 guez 3
220 guez 33 ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
221     ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &
222     / apols
223     END IF
224 guez 3
225 guez 33 IF (MOD(itau + 1, iperiod) == 0) THEN
226 guez 40 ! Écriture du fichier histoire moyenne:
227 guez 61 CALL writedynav(vcov, ucov, teta, pk, phi, q, masse, ps, phis, &
228     time = itau + 1)
229 guez 40 call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
230 guez 57 q(:, :, :, 1))
231 guez 33 ENDIF
232 guez 40 end do time_integration
233 guez 3
234 guez 30 CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
235     itau=itau_dyn+itaufin)
236    
237     ! Calcul des tendances dynamiques:
238 guez 55 CALL geopot((iim + 1) * (jjm + 1), teta, pk, pks, phis, phi)
239 guez 30 CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
240 guez 47 dudyn, dv, dteta, dp, w, pbaru, pbarv, time_0, &
241     conser=MOD(itaufin, iconser)==0)
242 guez 56
243 guez 3 END SUBROUTINE leapfrog
244    
245     end module leapfrog_m

  ViewVC Help
Powered by ViewVC 1.1.21