/[lmdze]/trunk/Sources/phylmd/phytrac.f
ViewVC logotype

Contents of /trunk/Sources/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 191 - (show annotations)
Mon May 9 19:56:28 2016 UTC (8 years ago) by guez
File size: 14879 byte(s)
Extracted the call to read_comdissnew out of conf_gcm.

Made ok_instan a variable of module clesphys, itau_phy a variable of
module phyetat0_m, nid_ins a variable of module ini_histins_m, itap a
variable of new module time_phylmdz, so that histwrite_phy can be
called from any procedure without the need to cascade those variables
into that procedure. Made itau_w a variable of module time_phylmdz so
that it is computed only once per time step of physics.

Extracted variables of module clesphys which were in namelist
conf_phys_nml into their own namelist, clesphys_nml, and created
procedure read_clesphys reading clesphys_nml, to avoid side effect.

No need for double precision in procedure getso4fromfile. Assume there
is a single variable for the whole year in the NetCDF file instead of
one variable per month.

Created generic procedure histwrite_phy and removed procedure
write_histins, following LMDZ. histwrite_phy has only two arguments,
can be called from anywhere, and should manage the logic of writing or
not writing into various history files with various operations. So the
test on ok_instan goes inside histwrite_phy.

Test for raz_date in phyetat0 instead of physiq to avoid side effect.

Created procedure increment_itap to avoid side effect.

Removed unnecessary differences between procedures readsulfate and
readsulfate_pi.

1 module phytrac_m
2
3 IMPLICIT none
4
5 private
6 public phytrac
7
8 contains
9
10 SUBROUTINE phytrac(lmt_pas, julien, gmtime, firstcal, lafin, pdtphys, &
11 t_seri, paprs, pplay, pmfu, pmfd, pde_u, pen_d, coefh, fm_therm, &
12 entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, &
13 mp, upwd, dnwd, tr_seri, zmasse, ncid_startphy)
14
15 ! From phylmd/phytrac.F, version 1.15 2006/02/21 08:08:30 (SVN revision 679)
16
17 ! Authors: Fr\'ed\'eric Hourdin, Abderrahmane Idelkadi, Marie-Alice
18 ! Foujols, Olivia
19
20 ! Objet : moniteur g\'en\'eral des tendances des traceurs
21
22 ! L'appel de "phytrac" se fait avec "nqmx - 2" donc nous avons
23 ! bien les vrais traceurs, sans la vapeur d'eau ni l'eau liquide.
24
25 ! Modifications pour les traceurs :
26 ! - uniformisation des parametrisations dans phytrac
27 ! - stockage des moyennes des champs n\'ecessaires en mode traceur off-line
28
29 use abort_gcm_m, only: abort_gcm
30 use clesphys, only: ecrit_tra
31 use clesphys2, only: conv_emanuel
32 use cltrac_m, only: cltrac
33 use cltracrn_m, only: cltracrn
34 use ctherm, only: iflag_thermals
35 use cvltr_m, only: cvltr
36 use dimens_m, only: llm, nqmx
37 use dimphy, only: klon
38 use indicesol, only: nbsrf
39 use initrrnpb_m, only: initrrnpb
40 use minmaxqfi_m, only: minmaxqfi
41 use netcdf, only: NF90_FILL_float
42 use netcdf95, only: nf95_inq_varid, nf95_get_var, nf95_put_var
43 use nflxtr_m, only: nflxtr
44 use nr_util, only: assert
45 use o3_chem_m, only: o3_chem
46 use phyetat0_m, only: rlat
47 use phyredem0_m, only: ncid_restartphy
48 use press_coefoz_m, only: press_coefoz
49 use radiornpb_m, only: radiornpb
50 use regr_pr_comb_coefoz_m, only: regr_pr_comb_coefoz
51 use SUPHEC_M, only: rg
52 use time_phylmdz, only: itap
53
54 integer, intent(in):: lmt_pas ! number of time steps of "physics" per day
55 integer, intent(in):: julien !jour julien, 1 <= julien <= 360
56 real, intent(in):: gmtime ! heure de la journ\'ee en fraction de jour
57 logical, intent(in):: firstcal ! first call to "calfis"
58 logical, intent(in):: lafin ! fin de la physique
59 real, intent(in):: pdtphys ! pas d'integration pour la physique (s)
60 real, intent(in):: t_seri(klon, llm) ! temperature, in K
61
62 real, intent(in):: paprs(klon, llm+1)
63 ! (pression pour chaque inter-couche, en Pa)
64
65 real, intent(in):: pplay(klon, llm)
66 ! (pression pour le mileu de chaque couche, en Pa)
67
68 ! convection:
69
70 REAL, intent(in):: pmfu(klon, llm) ! flux de masse dans le panache montant
71
72 REAL, intent(in):: pmfd(klon, llm)
73 ! flux de masse dans le panache descendant
74
75 REAL pde_u(klon, llm) ! flux detraine dans le panache montant
76 REAL pen_d(klon, llm) ! flux entraine dans le panache descendant
77 REAL coefh(klon, llm) ! coeff melange couche limite
78
79 ! thermiques:
80 real fm_therm(klon, llm+1), entr_therm(klon, llm)
81
82 ! Couche limite:
83 REAL yu1(klon) ! vents au premier niveau
84 REAL yv1(klon) ! vents au premier niveau
85
86 ! Arguments n\'ecessaires pour les sources et puits de traceur :
87 real, intent(in):: ftsol(klon, nbsrf) ! Temperature du sol (surf)(Kelvin)
88 real pctsrf(klon, nbsrf) ! Pourcentage de sol f(nature du sol)
89
90 ! Lessivage pour le on-line
91 REAL frac_impa(klon, llm) ! fraction d'aerosols impactes
92 REAL frac_nucl(klon, llm) ! fraction d'aerosols nuclees
93
94 ! Kerry Emanuel
95 real, intent(in):: da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
96 REAL, intent(in):: upwd(klon, llm) ! saturated updraft mass flux
97 REAL, intent(in):: dnwd(klon, llm) ! saturated downdraft mass flux
98
99 real, intent(inout):: tr_seri(:, :, :) ! (klon, llm, nqmx - 2)
100 ! (mass fractions of tracers, excluding water, at mid-layers)
101
102 real, intent(in):: zmasse(:, :) ! (klon, llm)
103 ! (column-density of mass of air in a cell, in kg m-2)
104
105 integer, intent(in):: ncid_startphy
106
107 ! Local:
108
109 integer nsplit
110
111 ! TRACEURS
112
113 ! Sources et puits des traceurs:
114
115 ! Pour l'instant seuls les cas du rn et du pb ont ete envisages.
116
117 REAL source(klon) ! a voir lorsque le flux est prescrit
118 !
119 ! Pour la source de radon et son reservoir de sol
120
121 REAL, save:: trs(klon, nqmx - 2) ! Concentration de traceur dans le sol
122
123 REAL masktr(klon, nqmx - 2) ! Masque reservoir de sol traceur
124 ! Masque de l'echange avec la surface
125 ! (1 = reservoir) ou (possible => 1)
126 SAVE masktr
127 REAL fshtr(klon, nqmx - 2) ! Flux surfacique dans le reservoir de sol
128 SAVE fshtr
129 REAL hsoltr(nqmx - 2) ! Epaisseur equivalente du reservoir de sol
130 SAVE hsoltr
131 REAL tautr(nqmx - 2) ! Constante de decroissance radioactive
132 SAVE tautr
133 REAL vdeptr(nqmx - 2) ! Vitesse de depot sec dans la couche Brownienne
134 SAVE vdeptr
135 REAL scavtr(nqmx - 2) ! Coefficient de lessivage
136 SAVE scavtr
137
138 CHARACTER itn
139
140 ! nature du traceur
141
142 logical aerosol(nqmx - 2) ! Nature du traceur
143 ! ! aerosol(it) = true => aerosol
144 ! ! aerosol(it) = false => gaz
145 logical clsol(nqmx - 2) ! couche limite sol calcul\'ee
146 logical radio(nqmx - 2) ! d\'ecroisssance radioactive
147 save aerosol, clsol, radio
148
149 ! convection tiedtke
150 INTEGER i, k, it
151 REAL delp(klon, llm)
152
153 ! Variables liees a l'ecriture de la bande histoire physique
154
155 ! Variables locales pour effectuer les appels en serie
156
157 REAL d_tr(klon, llm), d_trs(klon) ! tendances de traceurs
158 REAL d_tr_cl(klon, llm, nqmx - 2) ! tendance de traceurs couche limite
159
160 REAL d_tr_cv(klon, llm, nqmx - 2)
161 ! tendance de traceurs conv pour chq traceur
162
163 REAL d_tr_th(klon, llm, nqmx - 2) ! la tendance des thermiques
164 REAL d_tr_dec(klon, llm, 2) ! la tendance de la decroissance
165 ! ! radioactive du rn - > pb
166 REAL d_tr_lessi_impa(klon, llm, nqmx - 2) ! la tendance du lessivage
167 ! ! par impaction
168 REAL d_tr_lessi_nucl(klon, llm, nqmx - 2) ! la tendance du lessivage
169 ! ! par nucleation
170 REAL flestottr(klon, llm, nqmx - 2) ! flux de lessivage
171 ! ! dans chaque couche
172
173 real ztra_th(klon, llm)
174 integer isplit, varid
175
176 ! Controls:
177 logical:: couchelimite = .true.
178 logical:: convection = .true.
179 logical:: lessivage = .true.
180 logical, save:: inirnpb
181
182 !--------------------------------------
183
184 call assert(shape(zmasse) == (/klon, llm/), "phytrac zmasse")
185 call assert(shape(tr_seri) == (/klon, llm, nqmx - 2/), "phytrac tr_seri")
186
187 if (firstcal) then
188 print *, 'phytrac: pdtphys = ', pdtphys
189 PRINT *, 'Frequency of tracer output: ecrit_tra = ', ecrit_tra
190 inirnpb = .true.
191
192 ! Initialisation de certaines variables pour le radon et le plomb
193 ! Initialisation du traceur dans le sol (couche limite radonique)
194 trs(:, 2:) = 0.
195
196 call nf95_inq_varid(ncid_startphy, "trs", varid)
197 call nf95_get_var(ncid_startphy, varid, trs(:, 1))
198 if (any(trs(:, 1) == NF90_FILL_float)) call abort_gcm("phytrac", &
199 "some missing values in trs(:, 1)")
200
201 ! Initialisation de la fraction d'aerosols lessivee
202
203 d_tr_lessi_impa = 0.
204 d_tr_lessi_nucl = 0.
205
206 ! Initialisation de la nature des traceurs
207
208 DO it = 1, nqmx - 2
209 aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
210 radio(it) = .FALSE. ! par d\'efaut pas de passage par "radiornpb"
211 clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
212 ENDDO
213
214 if (nqmx >= 5) then
215 call press_coefoz ! read input pressure levels for ozone coefficients
216 end if
217 ENDIF
218
219 if (inirnpb) THEN
220 ! Initialisation du traceur dans le sol (couche limite radonique)
221 radio(1)= .true.
222 radio(2)= .true.
223 clsol(1)= .true.
224 clsol(2)= .true.
225 aerosol(2) = .TRUE. ! le Pb est un aerosol
226 call initrrnpb(pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, scavtr)
227 inirnpb=.false.
228 endif
229
230 if (convection) then
231 ! Calcul de l'effet de la convection
232 DO it=1, nqmx - 2
233 if (conv_emanuel) then
234 call cvltr(pdtphys, da, phi, mp, paprs, tr_seri(:, :, it), upwd, &
235 dnwd, d_tr_cv(:, :, it))
236 else
237 CALL nflxtr(pdtphys, pmfu, pmfd, pde_u, pen_d, paprs, &
238 tr_seri(:, :, it), d_tr_cv(:, :, it))
239 endif
240
241 DO k = 1, llm
242 DO i = 1, klon
243 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cv(i, k, it)
244 ENDDO
245 ENDDO
246 WRITE(unit=itn, fmt='(i1)') it
247 CALL minmaxqfi(tr_seri(:, :, it), 0., 1e33, &
248 'convection, tracer index = ' // itn)
249 ENDDO
250 endif
251
252 ! Calcul de l'effet des thermiques
253
254 do it=1, nqmx - 2
255 do k=1, llm
256 do i=1, klon
257 d_tr_th(i, k, it)=0.
258 tr_seri(i, k, it)=max(tr_seri(i, k, it), 0.)
259 tr_seri(i, k, it)=min(tr_seri(i, k, it), 1e10)
260 enddo
261 enddo
262 enddo
263
264 if (iflag_thermals > 0) then
265 nsplit=10
266 DO it=1, nqmx - 2
267 do isplit=1, nsplit
268 ! Thermiques
269 call dqthermcell(klon, llm, pdtphys/nsplit &
270 , fm_therm, entr_therm, zmasse &
271 , tr_seri(1:klon, 1:llm, it), d_tr, ztra_th)
272
273 do k=1, llm
274 do i=1, klon
275 d_tr(i, k)=pdtphys*d_tr(i, k)/nsplit
276 d_tr_th(i, k, it)=d_tr_th(i, k, it)+d_tr(i, k)
277 tr_seri(i, k, it)=max(tr_seri(i, k, it)+d_tr(i, k), 0.)
278 enddo
279 enddo
280 enddo
281 ENDDO
282 endif
283
284 ! Calcul de l'effet de la couche limite
285
286 if (couchelimite) then
287 DO k = 1, llm
288 DO i = 1, klon
289 delp(i, k) = paprs(i, k)-paprs(i, k+1)
290 ENDDO
291 ENDDO
292
293 ! MAF modif pour tenir compte du cas traceur
294 DO it=1, nqmx - 2
295 if (clsol(it)) then
296 ! couche limite avec quantite dans le sol calculee
297 CALL cltracrn(it, pdtphys, yu1, yv1, coefh, t_seri, ftsol, &
298 pctsrf, tr_seri(:, :, it), trs(:, it), paprs, pplay, delp, &
299 masktr(1, it), fshtr(1, it), hsoltr(it), tautr(it), &
300 vdeptr(it), rlat, d_tr_cl(1, 1, it), d_trs)
301 DO k = 1, llm
302 DO i = 1, klon
303 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
304 ENDDO
305 ENDDO
306
307 trs(:, it) = trs(:, it) + d_trs
308 else
309 ! couche limite avec flux prescrit
310 !MAF provisoire source / traceur a creer
311 DO i=1, klon
312 source(i) = 0. ! pas de source, pour l'instant
313 ENDDO
314
315 CALL cltrac(pdtphys, coefh, t_seri, tr_seri(:, :, it), source, &
316 paprs, pplay, delp, d_tr_cl(1, 1, it))
317 DO k = 1, llm
318 DO i = 1, klon
319 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
320 ENDDO
321 ENDDO
322 endif
323 ENDDO
324 endif
325
326 ! Calcul de l'effet du puits radioactif
327
328 ! MAF il faudrait faire une modification pour passer dans radiornpb
329 ! si radio=true
330 d_tr_dec = radiornpb(tr_seri, pdtphys, tautr)
331 DO it = 1, nqmx - 2
332 if (radio(it)) then
333 tr_seri(:, :, it) = tr_seri(:, :, it) + d_tr_dec(:, :, it)
334 WRITE(unit=itn, fmt='(i1)') it
335 CALL minmaxqfi(tr_seri(:, :, it), 0., 1e33, 'puits rn it='//itn)
336 endif
337 ENDDO
338
339 if (nqmx >= 5) then
340 ! Ozone as a tracer:
341 if (mod(itap - 1, lmt_pas) == 0) then
342 ! Once per day, update the coefficients for ozone chemistry:
343 call regr_pr_comb_coefoz(julien, paprs, pplay)
344 end if
345 call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
346 end if
347
348 ! Calcul de l'effet de la precipitation
349
350 IF (lessivage) THEN
351 d_tr_lessi_nucl = 0.
352 d_tr_lessi_impa = 0.
353 flestottr = 0.
354
355 ! tendance des aerosols nuclees et impactes
356
357 DO it = 1, nqmx - 2
358 IF (aerosol(it)) THEN
359 DO k = 1, llm
360 DO i = 1, klon
361 d_tr_lessi_nucl(i, k, it) = d_tr_lessi_nucl(i, k, it) + &
362 (1 - frac_nucl(i, k))*tr_seri(i, k, it)
363 d_tr_lessi_impa(i, k, it) = d_tr_lessi_impa(i, k, it) + &
364 (1 - frac_impa(i, k))*tr_seri(i, k, it)
365 ENDDO
366 ENDDO
367 ENDIF
368 ENDDO
369
370 ! Mises a jour des traceurs + calcul des flux de lessivage
371 ! Mise a jour due a l'impaction et a la nucleation
372
373 DO it = 1, nqmx - 2
374 IF (aerosol(it)) THEN
375 DO k = 1, llm
376 DO i = 1, klon
377 tr_seri(i, k, it) = tr_seri(i, k, it) * frac_impa(i, k) &
378 * frac_nucl(i, k)
379 ENDDO
380 ENDDO
381 ENDIF
382 ENDDO
383
384 ! Flux lessivage total
385
386 DO it = 1, nqmx - 2
387 DO k = 1, llm
388 DO i = 1, klon
389 flestottr(i, k, it) = flestottr(i, k, it) &
390 - (d_tr_lessi_nucl(i, k, it) + d_tr_lessi_impa(i, k, it)) &
391 * (paprs(i, k)-paprs(i, k+1)) / (RG * pdtphys)
392 ENDDO
393 ENDDO
394 ENDDO
395 ENDIF
396
397 ! Ecriture des sorties
398 call write_histrac(lessivage)
399
400 if (lafin) then
401 call nf95_inq_varid(ncid_restartphy, "trs", varid)
402 call nf95_put_var(ncid_restartphy, varid, trs(:, 1))
403 endif
404
405 contains
406
407 subroutine write_histrac(lessivage)
408
409 ! From phylmd/write_histrac.h, version 1.9 2006/02/21 08:08:30
410
411 use gr_phy_write_m, only: gr_phy_write
412 use histwrite_m, only: histwrite
413 use iniadvtrac_m, only: tname
414 use ini_histins_m, only: nid_ins
415 use phyetat0_m, only: itau_phy
416
417 logical, intent(in):: lessivage
418
419 ! Variables local to the procedure:
420 integer it
421 integer itau_w ! pas de temps ecriture
422
423 !-----------------------------------------------------
424
425 itau_w = itau_phy + itap
426
427 CALL histwrite(nid_ins, "zmasse", itau_w, gr_phy_write(zmasse))
428
429 DO it=1, nqmx - 2
430 CALL histwrite(nid_ins, tname(it+2), itau_w, &
431 gr_phy_write(tr_seri(:, :, it)))
432 if (lessivage) THEN
433 CALL histwrite(nid_ins, "fl"//tname(it+2), itau_w, &
434 gr_phy_write(flestottr(:, :, it)))
435 endif
436 CALL histwrite(nid_ins, "d_tr_th_"//tname(it+2), itau_w, &
437 gr_phy_write(d_tr_th(:, :, it)))
438 CALL histwrite(nid_ins, "d_tr_cv_"//tname(it+2), itau_w, &
439 gr_phy_write(d_tr_cv(:, :, it)))
440 CALL histwrite(nid_ins, "d_tr_cl_"//tname(it+2), itau_w, &
441 gr_phy_write(d_tr_cl(:, :, it)))
442 ENDDO
443
444 end subroutine write_histrac
445
446 END SUBROUTINE phytrac
447
448 end module phytrac_m

  ViewVC Help
Powered by ViewVC 1.1.21