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

Contents of /trunk/dyn3d/leapfrog.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 115 - (show annotations)
Fri Sep 19 17:36:20 2014 UTC (9 years, 7 months ago) by guez
File size: 8704 byte(s) 
Extracted code from tau2alpha for first call into new procedure
init_tau2alpha. dxdys, dxdyu, dxdyv are now local variables if guide
computed by init_tau2alpha. This allows us to remove terrible argument
type of tau2alpha: we just give to tau2alpha the right dxdy and
rlat.

In module conf_guide_m, changed default values of tau_min_*, because
0.02 is too small for the default daystep = 240, iperiod = 5. Changed
default values of guide_[uv] to false. Moved variable ok_guide from
conf_gcm_m to conf_guide_m, ok_guide is no longer an input parameter,
it is computed from guide_*. Had then to move test on ok_guide and
day_step from conf_gcm_m to conf_guide_m. Added checks on input
nudging parameters in procedure conf_guide. Upgraded variable factt to
module conf_guide_m in order to check nudging parameters. Bug fix:
variable guide_q was not in namelist conf_guide_nml.

Removed unused variables aire_min, aire_max of MODULE guide_m.

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