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

Contents of /trunk/dyn3d/leapfrog.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 32 - (show annotations)
Tue Apr 6 17:52:58 2010 UTC (14 years ago) by guez
Original Path: trunk/libf/dyn3d/leapfrog.f90
File size: 10090 byte(s)
Split "stringop.f90" into single-procedure files. Gathered files in directory
"IOIPSL/Stringop".

Split "flincom.f90" into "flincom.f90" and "flinget.f90". Removed
unused procedures from module "flincom". Removed unused argument
"filename" of procedure "flinopen_nozoom".

Removed unused files.

Split "grid_change.f90" into "grid_change.f90" and
"gr_phy_write_3d.f90".

Removed unused procedures from modules "calendar", "ioipslmpp",
"grid_atob", "gath_cpl" and "getincom". Removed unused procedures in
files "ppm3d.f" and "thermcell.f".

Split "mathelp.f90" into "mathelp.f90" and "mathop.f90".

Removed unused variable "dpres" of module "comvert".

Use argument "itau" instead of local variables "iadvtr" and "first" to
control algorithm in procedure "fluxstokenc".

Removed unused arguments of procedure "integrd".

Removed useless computations at the end of procedure "leapfrog".

Merged common block "matrfil" into module "parafilt".

1 module leapfrog_m
2
3 IMPLICIT NONE
4
5 contains
6
7 SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)
8
9 ! From dyn3d/leapfrog.F, version 1.6, 2005/04/13 08:58:34
10 ! Authors: P. Le Van, L. Fairhead, F. Hourdin
11 ! schema matsuno + leapfrog
12
13 USE calfis_m, ONLY: calfis
14 USE com_io_dyn, ONLY: histaveid
15 USE comconst, ONLY: daysec, dtphys, dtvr
16 USE comgeom, ONLY: aire_2d, apoln, apols
17 USE comvert, ONLY: ap, bp
18 USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
19 periodav
20 USE dimens_m, ONLY: iim, jjm, llm, nqmx
21 USE dynetat0_m, ONLY: day_ini
22 use dynredem1_m, only: dynredem1
23 USE exner_hyb_m, ONLY: exner_hyb
24 use filtreg_m, only: filtreg
25 USE guide_m, ONLY: guide
26 use inidissip_m, only: idissip
27 use integrd_m, only: integrd
28 USE logic, ONLY: iflag_phys, ok_guide
29 USE paramet_m, ONLY: ip1jmp1
30 USE pression_m, ONLY: pression
31 USE pressure_var, ONLY: p3d
32 USE temps, ONLY: itau_dyn
33
34 ! Variables dynamiques:
35 REAL, intent(inout):: vcov((iim + 1) * jjm, llm) ! vent covariant
36 REAL, intent(inout):: ucov(ip1jmp1, llm) ! vent covariant
37 REAL, intent(inout):: teta(iim + 1, jjm + 1, llm) ! potential temperature
38 REAL ps(iim + 1, jjm + 1) ! pression au sol, en Pa
39
40 REAL masse(ip1jmp1, llm) ! masse d'air
41 REAL phis(ip1jmp1) ! geopotentiel au sol
42 REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields
43 REAL, intent(in):: time_0
44
45 ! Variables local to the procedure:
46
47 ! Variables dynamiques:
48
49 REAL pks(ip1jmp1) ! exner au sol
50 REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
51 REAL pkf(ip1jmp1, llm) ! exner filt.au milieu des couches
52 REAL phi(ip1jmp1, llm) ! geopotential
53 REAL w(ip1jmp1, llm) ! vitesse verticale
54
55 ! variables dynamiques intermediaire pour le transport
56 REAL pbaru(ip1jmp1, llm), pbarv((iim + 1) * jjm, llm) !flux de masse
57
58 ! variables dynamiques au pas - 1
59 REAL vcovm1((iim + 1) * jjm, llm), ucovm1(ip1jmp1, llm)
60 REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
61 REAL massem1(ip1jmp1, llm)
62
63 ! tendances dynamiques
64 REAL dv((iim + 1) * jjm, llm), du(ip1jmp1, llm)
65 REAL dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqmx), dp(ip1jmp1)
66
67 ! tendances de la dissipation
68 REAL dvdis((iim + 1) * jjm, llm), dudis(ip1jmp1, llm)
69 REAL dtetadis(iim + 1, jjm + 1, llm)
70
71 ! tendances physiques
72 REAL dvfi((iim + 1) * jjm, llm), dufi(ip1jmp1, llm)
73 REAL dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqmx), dpfi(ip1jmp1)
74
75 ! variables pour le fichier histoire
76
77 INTEGER itau ! index of the time step of the dynamics, starts at 0
78 INTEGER itaufin
79 INTEGER iday ! jour julien
80 REAL time ! time of day, as a fraction of day length
81 real finvmaold(ip1jmp1, llm)
82 LOGICAL:: lafin=.false.
83 INTEGER i, j, l
84
85 REAL rdayvrai, rdaym_ini
86
87 ! Variables test conservation energie
88 REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
89 ! Tendance de la temp. potentiel d (theta) / d t due a la
90 ! tansformation d'energie cinetique en energie thermique
91 ! cree par la dissipation
92 REAL dtetaecdt(iim + 1, jjm + 1, llm)
93 REAL vcont((iim + 1) * jjm, llm), ucont(ip1jmp1, llm)
94
95 !---------------------------------------------------
96
97 print *, "Call sequence information: leapfrog"
98
99 itaufin = nday * day_step
100 ! "day_step" is a multiple of "iperiod", therefore "itaufin" is one too
101
102 itau = 0
103 iday = day_ini
104 time = time_0
105 dq = 0.
106 ! On initialise la pression et la fonction d'Exner :
107 CALL pression(ip1jmp1, ap, bp, ps, p3d)
108 CALL exner_hyb(ps, p3d, pks, pk, pkf)
109
110 ! Début de l'integration temporelle :
111 period_loop:do i = 1, itaufin / iperiod
112 ! {"itau" is a multiple of "iperiod"}
113
114 ! 1. Matsuno forward:
115
116 if (ok_guide .and. (itaufin - itau - 1) * dtvr > 21600.) &
117 call guide(itau, ucov, vcov, teta, q, masse, ps)
118 vcovm1 = vcov
119 ucovm1 = ucov
120 tetam1 = teta
121 massem1 = masse
122 psm1 = ps
123 finvmaold = masse
124 CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)
125
126 ! Calcul des tendances dynamiques:
127 CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
128 CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
129 MOD(itau, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
130 time + iday - day_ini)
131
132 ! Calcul des tendances advection des traceurs (dont l'humidité)
133 CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
134 ! Stokage du flux de masse pour traceurs offline:
135 IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
136 dtvr, itau)
137
138 ! integrations dynamique et traceurs:
139 CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &
140 dp, vcov, ucov, teta, q, ps, masse, finvmaold, .false., &
141 dtvr)
142
143 CALL pression(ip1jmp1, ap, bp, ps, p3d)
144 CALL exner_hyb(ps, p3d, pks, pk, pkf)
145
146 ! 2. Matsuno backward:
147
148 itau = itau + 1
149 iday = day_ini + itau / day_step
150 time = REAL(itau - (iday - day_ini) * day_step) / day_step + time_0
151 IF (time > 1.) THEN
152 time = time - 1.
153 iday = iday + 1
154 ENDIF
155
156 ! Calcul des tendances dynamiques:
157 CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
158 CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
159 .false., du, dv, dteta, dp, w, pbaru, pbarv, time + iday - day_ini)
160
161 ! integrations dynamique et traceurs:
162 CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &
163 dp, vcov, ucov, teta, q, ps, masse, finvmaold, .false., &
164 dtvr)
165
166 CALL pression(ip1jmp1, ap, bp, ps, p3d)
167 CALL exner_hyb(ps, p3d, pks, pk, pkf)
168
169 ! 3. Leapfrog:
170
171 leapfrog_loop: do j = 1, iperiod - 1
172 ! Calcul des tendances dynamiques:
173 CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
174 CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
175 .false., du, dv, dteta, dp, w, pbaru, pbarv, &
176 time + iday - day_ini)
177
178 ! Calcul des tendances advection des traceurs (dont l'humidité)
179 CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk)
180 ! Stokage du flux de masse pour traceurs off-line:
181 IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
182 dtvr, itau)
183
184 ! integrations dynamique et traceurs:
185 CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
186 dteta, dp, vcov, ucov, teta, q, ps, masse, &
187 finvmaold, .true., 2 * dtvr)
188
189 IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
190 ! calcul des tendances physiques:
191 IF (itau + 1 == itaufin) lafin = .TRUE.
192
193 CALL pression(ip1jmp1, ap, bp, ps, p3d)
194 CALL exner_hyb(ps, p3d, pks, pk, pkf)
195
196 rdaym_ini = itau * dtvr / daysec
197 rdayvrai = rdaym_ini + day_ini
198
199 CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, &
200 masse, ps, pk, phis, phi, du, dv, dteta, dq, w, &
201 dufi, dvfi, dtetafi, dqfi, dpfi)
202
203 ! ajout des tendances physiques:
204 CALL addfi(nqmx, dtphys, ucov, vcov, teta, q, ps, dufi, dvfi, &
205 dtetafi, dqfi, dpfi)
206 ENDIF
207
208 CALL pression(ip1jmp1, ap, bp, ps, p3d)
209 CALL exner_hyb(ps, p3d, pks, pk, pkf)
210
211 IF (MOD(itau + 1, idissip) == 0) THEN
212 ! dissipation horizontale et verticale des petites echelles:
213
214 ! calcul de l'energie cinetique avant dissipation
215 call covcont(llm, ucov, vcov, ucont, vcont)
216 call enercin(vcov, ucov, vcont, ucont, ecin0)
217
218 ! dissipation
219 CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
220 ucov=ucov + dudis
221 vcov=vcov + dvdis
222
223 ! On rajoute la tendance due à la transformation Ec -> E
224 ! thermique créée lors de la dissipation
225 call covcont(llm, ucov, vcov, ucont, vcont)
226 call enercin(vcov, ucov, vcont, ucont, ecin)
227 dtetaecdt= (ecin0 - ecin) / pk
228 dtetadis=dtetadis + dtetaecdt
229 teta=teta + dtetadis
230
231 ! Calcul de la valeur moyenne aux pôles :
232 forall (l = 1: llm)
233 teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
234 / apoln
235 teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
236 * teta(:iim, jjm + 1, l)) / apols
237 END forall
238
239 ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
240 ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm+1) * ps(:iim, jjm + 1)) &
241 / apols
242 END IF
243
244 itau = itau + 1
245 iday = day_ini + itau / day_step
246 time = REAL(itau - (iday - day_ini) * day_step) / day_step + time_0
247 IF (time > 1.) THEN
248 time = time - 1.
249 iday = iday + 1
250 ENDIF
251
252 IF (MOD(itau, iperiod) == 0) THEN
253 ! ecriture du fichier histoire moyenne:
254 CALL writedynav(histaveid, nqmx, itau, vcov, &
255 ucov, teta, pk, phi, q, masse, ps, phis)
256 call bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
257 ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
258 ENDIF
259 end do leapfrog_loop
260 end do period_loop
261
262 ! {itau == itaufin}
263 CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps, &
264 itau=itau_dyn+itaufin)
265
266 ! Calcul des tendances dynamiques:
267 CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
268 CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
269 MOD(itaufin, iconser) == 0, du, dv, dteta, dp, w, pbaru, pbarv, &
270 time + iday - day_ini)
271
272 END SUBROUTINE leapfrog
273
274 end module leapfrog_m

  ViewVC Help
Powered by ViewVC 1.1.21