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

Annotation of /trunk/Sources/dyn3d/leapfrog.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 56 - (hide annotations)
Tue Jan 10 19:02:02 2012 UTC (12 years, 4 months ago) by guez
Original Path: trunk/libf/dyn3d/leapfrog.f90
File size: 9064 byte(s)
Imported "writehist.f" from LMDZ.

Moved module variable "histaveid" from "com_io_dyn" to "initdynav_m".

In "inithist", access directly module variables from "com_io_dyn"
instead of going through the arguments. Copying from LMDZ, write "u"
and scalar variables to separate files. Create a new variable for the
new file in "com_io_dyn". Copying from LMDZ, change the vertical axes
of the three files.

Removed some useless initializations in "dissip".

In "bilan_dyn", removed useless variable "time". Avoiding the
approximate test on "dt_cum" being a multiple of "dt_app", just
compute "ncum" from known usage of "bilan_dyn" and compute "dt_cum"
from "ncum". Change "periodav" from real to integer in
"conf_gcm_m". Since "day_step" is required to be a multiple of
"iperiod", so is "ncum".

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 56 USE conf_gcm_m, ONLY: day_step, iconser, iperiod, iphysiq, nday, offline
22 guez 29 USE dimens_m, ONLY: iim, jjm, llm, nqmx
23 guez 47 use dissip_m, only: dissip
24 guez 26 USE dynetat0_m, ONLY: day_ini
25 guez 27 use dynredem1_m, only: dynredem1
26 guez 26 USE exner_hyb_m, ONLY: exner_hyb
27 guez 27 use filtreg_m, only: filtreg
28 guez 43 use geopot_m, only: geopot
29 guez 26 USE guide_m, ONLY: guide
30     use inidissip_m, only: idissip
31 guez 32 use integrd_m, only: integrd
32 guez 26 USE logic, ONLY: iflag_phys, ok_guide
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 30 ! "day_step" is a multiple of "iperiod", therefore "itaufin" is one too
110    
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     CALL filtreg(finvmaold, jjm + 1, llm, - 2, 2, .TRUE., 1)
133     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 56 CALL writedynav(nqmx, itau + 1, vcov, ucov, teta, pk, phi, q, &
228     masse, ps, phis)
229 guez 40 call bilan_dyn(ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, &
230 guez 56 q(:, :, :, 1), dt_app = dtvr * iperiod)
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