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

Contents of /trunk/dyn3d/leapfrog.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 157 - (show annotations)
Mon Jul 20 16:01:49 2015 UTC (8 years, 10 months ago) by guez
Original Path: trunk/Sources/dyn3d/leapfrog.f
File size: 8618 byte(s) 
Just encapsulated SUBROUTINE vlsplt in a module and cleaned it.

In procedure vlx, local variables dxqu and adxqu only need indices
iip2:ip1jm. Otherwise, just cleaned vlx.

Procedures dynredem0 and dynredem1 no longer have argument fichnom,
they just operate on a file named "restart.nc". The programming
guideline here is that gcm should not be more complex than it needs by
itself, other programs (ce0l etc.) just have to adapt to gcm. So ce0l
now creates files "restart.nc" and "restartphy.nc".

In order to facilitate decentralizing the writing of "restartphy.nc",
created a procedure phyredem0 out of phyredem. phyredem0 creates the
NetCDF header of "restartphy.nc" while phyredem writes the NetCDF
variables. As the global attribute itau_phy needs to be filled in
phyredem0, at the beginnig of the run, we must compute its value
instead of just using itap. So we have a dummy argument lmt_pas of
phyredem0. Also, the ncid of "startphy.nc" is upgraded from local
variable of phyetat0 to dummy argument. phyetat0 no longer closes
"startphy.nc".

Following the same decentralizing objective, the ncid of "restart.nc"
is upgraded from local variable of dynredem0 to module variable of
dynredem0_m. "restart.nc" is not closed at the end of dynredem0 nor
opened at the beginning of dynredem1.

In procedure etat0, instead of creating many vectors of size klon
which will be filled with zeroes, just create one array null_array.

In procedure phytrac, instead of writing trs(: 1) to a text file,
write it to "restartphy.nc" (following LMDZ). This is better because
now trs(: 1) is next to its coordinates. We can write to
"restartphy.nc" from phytrac directly, and not add trs(: 1) to the
long list of variables in physiq, thanks to the decentralizing of
"restartphy.nc".

In procedure phyetat0, we no longer write to standard output the
minimum and maximum values of read arrays. It is ok to check input and
abort on invalid values but just printing statistics on input seems too
much useless computation and out of place clutter.
1 module leapfrog_m
2
3 IMPLICIT NONE
4
5 contains
6
7 SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q)
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
12 ! Intégration temporelle du modèle : Matsuno-leapfrog scheme.
13
14 use addfi_m, only: addfi
15 use bilan_dyn_m, only: bilan_dyn
16 use caladvtrac_m, only: caladvtrac
17 use caldyn_m, only: caldyn
18 USE calfis_m, ONLY: calfis
19 USE comconst, ONLY: daysec, dtvr
20 USE comgeom, ONLY: aire_2d, apoln, apols
21 USE disvert_m, ONLY: ap, bp
22 USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline, &
23 iflag_phys, iecri
24 USE conf_guide_m, ONLY: ok_guide
25 USE dimens_m, ONLY: iim, jjm, llm, nqmx
26 use dissip_m, only: dissip
27 USE dynetat0_m, ONLY: day_ini
28 use dynredem1_m, only: dynredem1
29 USE exner_hyb_m, ONLY: exner_hyb
30 use filtreg_scal_m, only: filtreg_scal
31 use fluxstokenc_m, only: fluxstokenc
32 use geopot_m, only: geopot
33 USE guide_m, ONLY: guide
34 use inidissip_m, only: idissip
35 use integrd_m, only: integrd
36 use nr_util, only: assert
37 USE pressure_var, ONLY: p3d
38 USE temps, ONLY: itau_dyn
39 use writedynav_m, only: writedynav
40 use writehist_m, only: writehist
41
42 ! Variables dynamiques:
43 REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
44 REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
45
46 REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
47 ! potential temperature
48
49 REAL, intent(inout):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol, en Pa
50 REAL, intent(inout):: masse(:, :, :) ! (iim + 1, jjm + 1, llm) masse d'air
51 REAL, intent(in):: phis(:, :) ! (iim + 1, jjm + 1) surface geopotential
52
53 REAL, intent(inout):: q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nqmx)
54 ! mass fractions of advected fields
55
56 ! Local:
57
58 ! Variables dynamiques:
59
60 REAL pks(iim + 1, jjm + 1) ! exner au sol
61 REAL pk(iim + 1, jjm + 1, llm) ! exner au milieu des couches
62 REAL pkf(iim + 1, jjm + 1, llm) ! exner filtr\'e au milieu des couches
63 REAL phi(iim + 1, jjm + 1, llm) ! geopotential
64 REAL w(iim + 1, jjm + 1, llm) ! vitesse verticale
65
66 ! Variables dynamiques intermediaire pour le transport
67 ! Flux de masse :
68 REAL pbaru(iim + 1, jjm + 1, llm), pbarv(iim + 1, jjm, llm)
69
70 ! Variables dynamiques au pas - 1
71 REAL vcovm1(iim + 1, jjm, llm), ucovm1(iim + 1, jjm + 1, llm)
72 REAL tetam1(iim + 1, jjm + 1, llm), psm1(iim + 1, jjm + 1)
73 REAL massem1(iim + 1, jjm + 1, llm)
74
75 ! Tendances dynamiques
76 REAL dv((iim + 1) * jjm, llm), dudyn(iim + 1, jjm + 1, llm)
77 REAL dteta(iim + 1, jjm + 1, llm)
78 real dp((iim + 1) * (jjm + 1))
79
80 ! Tendances de la dissipation :
81 REAL dvdis(iim + 1, jjm, llm), dudis(iim + 1, jjm + 1, llm)
82 REAL dtetadis(iim + 1, jjm + 1, llm)
83
84 ! Tendances physiques
85 REAL dvfi(iim + 1, jjm, llm), dufi(iim + 1, jjm + 1, llm)
86 REAL dtetafi(iim + 1, jjm + 1, llm), dqfi(iim + 1, jjm + 1, llm, nqmx)
87
88 ! Variables pour le fichier histoire
89 INTEGER itau ! index of the time step of the dynamics, starts at 0
90 INTEGER itaufin
91 real finvmaold(iim + 1, jjm + 1, llm)
92 INTEGER l
93
94 ! Variables test conservation \'energie
95 REAL ecin(iim + 1, jjm + 1, llm), ecin0(iim + 1, jjm + 1, llm)
96
97 REAL vcont((iim + 1) * jjm, llm), ucont((iim + 1) * (jjm + 1), llm)
98 logical leapf
99 real dt ! time step, in s
100
101 !---------------------------------------------------
102
103 print *, "Call sequence information: leapfrog"
104 call assert(shape(ucov) == (/iim + 1, jjm + 1, llm/), "leapfrog")
105
106 itaufin = nday * day_step
107 ! "day_step" is a multiple of "iperiod", therefore so is "itaufin".
108
109 ! On initialise la pression et la fonction d'Exner :
110 forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
111 CALL exner_hyb(ps, p3d, pks, pk)
112 pkf = pk
113 CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
114
115 time_integration: do itau = 0, itaufin - 1
116 leapf = mod(itau, iperiod) /= 0
117 if (leapf) then
118 dt = 2 * dtvr
119 else
120 ! Matsuno
121 dt = dtvr
122 if (ok_guide) call guide(itau, ucov, vcov, teta, q(:, :, :, 1), ps)
123 vcovm1 = vcov
124 ucovm1 = ucov
125 tetam1 = teta
126 massem1 = masse
127 psm1 = ps
128 finvmaold = masse
129 CALL filtreg_scal(finvmaold, direct = .false., intensive = .false.)
130 end if
131
132 ! Calcul des tendances dynamiques:
133 CALL geopot(teta, pk, pks, phis, phi)
134 CALL caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
135 dudyn, dv, dteta, dp, w, pbaru, pbarv, &
136 conser = MOD(itau, iconser) == 0)
137
138 CALL caladvtrac(q, pbaru, pbarv, p3d, masse, teta, pk)
139
140 ! Stokage du flux de masse pour traceurs offline:
141 IF (offline) CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
142 dtvr, itau)
143
144 ! Int\'egrations dynamique et traceurs:
145 CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, dteta, &
146 dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, finvmaold, dt, &
147 leapf)
148
149 forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
150 CALL exner_hyb(ps, p3d, pks, pk)
151 pkf = pk
152 CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
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, conser = .false.)
160
161 ! integrations dynamique et traceurs:
162 CALL integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, &
163 dteta, dp, vcov, ucov, teta, q(:, :, :, :2), ps, masse, &
164 finvmaold, dtvr, leapf=.false.)
165
166 forall (l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
167 CALL exner_hyb(ps, p3d, pks, pk)
168 pkf = pk
169 CALL filtreg_scal(pkf, direct = .true., intensive = .true.)
170 end if
171
172 IF (MOD(itau + 1, iphysiq) == 0 .AND. iflag_phys /= 0) THEN
173 CALL calfis(ucov, vcov, teta, q, pk, phis, phi, w, dufi, dvfi, &
174 dtetafi, dqfi, dayvrai = itau / day_step + day_ini, &
175 time = REAL(mod(itau, day_step)) / day_step, &
176 lafin = itau + 1 == itaufin)
177
178 CALL addfi(ucov, vcov, teta, q, dufi, dvfi, dtetafi, dqfi)
179 ENDIF
180
181 IF (MOD(itau + 1, idissip) == 0) THEN
182 ! Dissipation horizontale et verticale des petites \'echelles
183
184 ! calcul de l'\'energie cin\'etique avant dissipation
185 call covcont(llm, ucov, vcov, ucont, vcont)
186 call enercin(vcov, ucov, vcont, ucont, ecin0)
187
188 ! dissipation
189 CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis)
190 ucov = ucov + dudis
191 vcov = vcov + dvdis
192
193 ! On ajoute la tendance due \`a la transformation \'energie
194 ! cin\'etique en \'energie thermique par la dissipation
195 call covcont(llm, ucov, vcov, ucont, vcont)
196 call enercin(vcov, ucov, vcont, ucont, ecin)
197 dtetadis = dtetadis + (ecin0 - ecin) / pk
198 teta = teta + dtetadis
199
200 ! Calcul de la valeur moyenne aux p\^oles :
201 forall (l = 1: llm)
202 teta(:, 1, l) = SUM(aire_2d(:iim, 1) * teta(:iim, 1, l)) &
203 / apoln
204 teta(:, jjm + 1, l) = SUM(aire_2d(:iim, jjm+1) &
205 * teta(:iim, jjm + 1, l)) / apols
206 END forall
207 END IF
208
209 IF (MOD(itau + 1, iperiod) == 0) THEN
210 ! \'Ecriture du fichier histoire moyenne:
211 CALL writedynav(vcov, ucov, teta, pk, phi, q, masse, ps, phis, &
212 time = itau + 1)
213 call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
214 q(:, :, :, 1))
215 ENDIF
216
217 IF (MOD(itau + 1, iecri * day_step) == 0) THEN
218 CALL geopot(teta, pk, pks, phis, phi)
219 CALL writehist(itau, vcov, ucov, teta, phi, masse, ps)
220 END IF
221 end do time_integration
222
223 CALL dynredem1(vcov, ucov, teta, q, masse, ps, itau = itau_dyn + itaufin)
224
225 ! Calcul des tendances dynamiques:
226 CALL geopot(teta, pk, pks, phis, phi)
227 CALL caldyn(itaufin, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, &
228 dudyn, dv, dteta, dp, w, pbaru, pbarv, &
229 conser = MOD(itaufin, iconser) == 0)
230
231 END SUBROUTINE leapfrog
232
233 end module leapfrog_m
  ViewVC Help
Powered by ViewVC 1.1.21