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

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

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
File size: 15713 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 phytrac_m
2
3 IMPLICIT none
4
5 private
6 public phytrac
7
8 contains
9
10 SUBROUTINE phytrac(itap, 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, pphis, da, &
13 phi, 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: iflag_con
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 ini_histrac_m, only: ini_histrac
40 use initrrnpb_m, only: initrrnpb
41 use minmaxqfi_m, only: minmaxqfi
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
53 integer, intent(in):: itap ! number of calls to "physiq"
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 real, intent(in):: pphis(klon)
95
96 ! Kerry Emanuel
97 real, intent(in):: da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
98 REAL, intent(in):: upwd(klon, llm) ! saturated updraft mass flux
99 REAL, intent(in):: dnwd(klon, llm) ! saturated downdraft mass flux
100
101 real, intent(inout):: tr_seri(:, :, :) ! (klon, llm, nqmx - 2)
102 ! (mass fractions of tracers, excluding water, at mid-layers)
103
104 real, intent(in):: zmasse(:, :) ! (klon, llm)
105 ! (column-density of mass of air in a cell, in kg m-2)
106
107 integer, intent(in):: ncid_startphy
108
109 ! Local:
110
111 integer nsplit
112
113 ! TRACEURS
114
115 ! Sources et puits des traceurs:
116
117 ! Pour l'instant seuls les cas du rn et du pb ont ete envisages.
118
119 REAL source(klon) ! a voir lorsque le flux est prescrit
120 !
121 ! Pour la source de radon et son reservoir de sol
122
123 REAL, save:: trs(klon, nqmx - 2) ! Concentration de traceur dans le sol
124
125 REAL masktr(klon, nqmx - 2) ! Masque reservoir de sol traceur
126 ! Masque de l'echange avec la surface
127 ! (1 = reservoir) ou (possible => 1)
128 SAVE masktr
129 REAL fshtr(klon, nqmx - 2) ! Flux surfacique dans le reservoir de sol
130 SAVE fshtr
131 REAL hsoltr(nqmx - 2) ! Epaisseur equivalente du reservoir de sol
132 SAVE hsoltr
133 REAL tautr(nqmx - 2) ! Constante de decroissance radioactive
134 SAVE tautr
135 REAL vdeptr(nqmx - 2) ! Vitesse de depot sec dans la couche Brownienne
136 SAVE vdeptr
137 REAL scavtr(nqmx - 2) ! Coefficient de lessivage
138 SAVE scavtr
139
140 CHARACTER itn
141 INTEGER, save:: nid_tra
142
143 ! nature du traceur
144
145 logical aerosol(nqmx - 2) ! Nature du traceur
146 ! ! aerosol(it) = true => aerosol
147 ! ! aerosol(it) = false => gaz
148 logical clsol(nqmx - 2) ! couche limite sol calcul\'ee
149 logical radio(nqmx - 2) ! d\'ecroisssance radioactive
150 save aerosol, clsol, radio
151
152 ! convection tiedtke
153 INTEGER i, k, it
154 REAL delp(klon, llm)
155
156 ! Variables liees a l'ecriture de la bande histoire physique
157
158 ! Variables locales pour effectuer les appels en serie
159
160 REAL d_tr(klon, llm), d_trs(klon) ! tendances de traceurs
161 REAL d_tr_cl(klon, llm, nqmx - 2) ! tendance de traceurs couche limite
162
163 REAL d_tr_cv(klon, llm, nqmx - 2)
164 ! tendance de traceurs conv pour chq traceur
165
166 REAL d_tr_th(klon, llm, nqmx - 2) ! la tendance des thermiques
167 REAL d_tr_dec(klon, llm, 2) ! la tendance de la decroissance
168 ! ! radioactive du rn - > pb
169 REAL d_tr_lessi_impa(klon, llm, nqmx - 2) ! la tendance du lessivage
170 ! ! par impaction
171 REAL d_tr_lessi_nucl(klon, llm, nqmx - 2) ! la tendance du lessivage
172 ! ! par nucleation
173 REAL flestottr(klon, llm, nqmx - 2) ! flux de lessivage
174 ! ! dans chaque couche
175
176 real ztra_th(klon, llm)
177 integer isplit, varid
178
179 ! Controls:
180 logical:: couchelimite = .true.
181 logical:: convection = .true.
182 logical:: lessivage = .true.
183 logical, save:: inirnpb
184
185 !--------------------------------------
186
187 call assert(shape(zmasse) == (/klon, llm/), "phytrac zmasse")
188 call assert(shape(tr_seri) == (/klon, llm, nqmx - 2/), "phytrac tr_seri")
189
190 if (firstcal) then
191 print *, 'phytrac: pdtphys = ', pdtphys
192 PRINT *, 'Frequency of tracer output: ecrit_tra = ', ecrit_tra
193 inirnpb = .true.
194
195 ! Initialisation des sorties :
196 call ini_histrac(nid_tra, pdtphys, nqmx - 2, lessivage)
197
198 ! Initialisation de certaines variables pour le radon et le plomb
199 ! Initialisation du traceur dans le sol (couche limite radonique)
200 trs(:, 2:) = 0.
201
202 call nf95_inq_varid(ncid_startphy, "trs", varid)
203 call nf95_get_var(ncid_startphy, varid, trs(:, 1))
204
205 ! Initialisation de la fraction d'aerosols lessivee
206
207 d_tr_lessi_impa = 0.
208 d_tr_lessi_nucl = 0.
209
210 ! Initialisation de la nature des traceurs
211
212 DO it = 1, nqmx - 2
213 aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
214 radio(it) = .FALSE. ! par d\'efaut pas de passage par "radiornpb"
215 clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
216 ENDDO
217
218 if (nqmx >= 5) then
219 call press_coefoz ! read input pressure levels for ozone coefficients
220 end if
221 ENDIF
222
223 if (inirnpb) THEN
224 ! Initialisation du traceur dans le sol (couche limite radonique)
225 radio(1)= .true.
226 radio(2)= .true.
227 clsol(1)= .true.
228 clsol(2)= .true.
229 aerosol(2) = .TRUE. ! le Pb est un aerosol
230 call initrrnpb(pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, scavtr)
231 inirnpb=.false.
232 endif
233
234 if (convection) then
235 ! Calcul de l'effet de la convection
236 DO it=1, nqmx - 2
237 if (iflag_con == 2) then
238 ! Tiedke
239 CALL nflxtr(pdtphys, pmfu, pmfd, pde_u, pen_d, paprs, &
240 tr_seri(:, :, it), d_tr_cv(:, :, it))
241 else if (iflag_con == 3) then
242 ! Emanuel
243 call cvltr(pdtphys, da, phi, mp, paprs, tr_seri(:, :, it), upwd, &
244 dnwd, d_tr_cv(:, :, it))
245 endif
246
247 DO k = 1, llm
248 DO i = 1, klon
249 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cv(i, k, it)
250 ENDDO
251 ENDDO
252 WRITE(unit=itn, fmt='(i1)') it
253 CALL minmaxqfi(tr_seri(:, :, it), 0., 1e33, &
254 'convection, tracer index = ' // itn)
255 ENDDO
256 endif
257
258 ! Calcul de l'effet des thermiques
259
260 do it=1, nqmx - 2
261 do k=1, llm
262 do i=1, klon
263 d_tr_th(i, k, it)=0.
264 tr_seri(i, k, it)=max(tr_seri(i, k, it), 0.)
265 tr_seri(i, k, it)=min(tr_seri(i, k, it), 1e10)
266 enddo
267 enddo
268 enddo
269
270 if (iflag_thermals > 0) then
271 nsplit=10
272 DO it=1, nqmx - 2
273 do isplit=1, nsplit
274 ! Thermiques
275 call dqthermcell(klon, llm, pdtphys/nsplit &
276 , fm_therm, entr_therm, zmasse &
277 , tr_seri(1:klon, 1:llm, it), d_tr, ztra_th)
278
279 do k=1, llm
280 do i=1, klon
281 d_tr(i, k)=pdtphys*d_tr(i, k)/nsplit
282 d_tr_th(i, k, it)=d_tr_th(i, k, it)+d_tr(i, k)
283 tr_seri(i, k, it)=max(tr_seri(i, k, it)+d_tr(i, k), 0.)
284 enddo
285 enddo
286 enddo
287 ENDDO
288 endif
289
290 ! Calcul de l'effet de la couche limite
291
292 if (couchelimite) then
293 DO k = 1, llm
294 DO i = 1, klon
295 delp(i, k) = paprs(i, k)-paprs(i, k+1)
296 ENDDO
297 ENDDO
298
299 ! MAF modif pour tenir compte du cas traceur
300 DO it=1, nqmx - 2
301 if (clsol(it)) then
302 ! couche limite avec quantite dans le sol calculee
303 CALL cltracrn(it, pdtphys, yu1, yv1, coefh, t_seri, ftsol, &
304 pctsrf, tr_seri(:, :, it), trs(:, it), paprs, pplay, delp, &
305 masktr(1, it), fshtr(1, it), hsoltr(it), tautr(it), &
306 vdeptr(it), rlat, d_tr_cl(1, 1, it), d_trs)
307 DO k = 1, llm
308 DO i = 1, klon
309 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
310 ENDDO
311 ENDDO
312
313 trs(:, it) = trs(:, it) + d_trs
314 else
315 ! couche limite avec flux prescrit
316 !MAF provisoire source / traceur a creer
317 DO i=1, klon
318 source(i) = 0. ! pas de source, pour l'instant
319 ENDDO
320
321 CALL cltrac(pdtphys, coefh, t_seri, tr_seri(:, :, it), source, &
322 paprs, pplay, delp, d_tr_cl(1, 1, it))
323 DO k = 1, llm
324 DO i = 1, klon
325 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
326 ENDDO
327 ENDDO
328 endif
329 ENDDO
330 endif
331
332 ! Calcul de l'effet du puits radioactif
333
334 ! MAF il faudrait faire une modification pour passer dans radiornpb
335 ! si radio=true
336 d_tr_dec = radiornpb(tr_seri, pdtphys, tautr)
337 DO it = 1, nqmx - 2
338 if (radio(it)) then
339 tr_seri(:, :, it) = tr_seri(:, :, it) + d_tr_dec(:, :, it)
340 WRITE(unit=itn, fmt='(i1)') it
341 CALL minmaxqfi(tr_seri(:, :, it), 0., 1e33, 'puits rn it='//itn)
342 endif
343 ENDDO
344
345 if (nqmx >= 5) then
346 ! Ozone as a tracer:
347 if (mod(itap - 1, lmt_pas) == 0) then
348 ! Once per day, update the coefficients for ozone chemistry:
349 call regr_pr_comb_coefoz(julien)
350 end if
351 call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
352 end if
353
354 ! Calcul de l'effet de la precipitation
355
356 IF (lessivage) THEN
357 d_tr_lessi_nucl = 0.
358 d_tr_lessi_impa = 0.
359 flestottr = 0.
360
361 ! tendance des aerosols nuclees et impactes
362
363 DO it = 1, nqmx - 2
364 IF (aerosol(it)) THEN
365 DO k = 1, llm
366 DO i = 1, klon
367 d_tr_lessi_nucl(i, k, it) = d_tr_lessi_nucl(i, k, it) + &
368 (1 - frac_nucl(i, k))*tr_seri(i, k, it)
369 d_tr_lessi_impa(i, k, it) = d_tr_lessi_impa(i, k, it) + &
370 (1 - frac_impa(i, k))*tr_seri(i, k, it)
371 ENDDO
372 ENDDO
373 ENDIF
374 ENDDO
375
376 ! Mises a jour des traceurs + calcul des flux de lessivage
377 ! Mise a jour due a l'impaction et a la nucleation
378
379 DO it = 1, nqmx - 2
380 IF (aerosol(it)) THEN
381 DO k = 1, llm
382 DO i = 1, klon
383 tr_seri(i, k, it) = tr_seri(i, k, it) * frac_impa(i, k) &
384 * frac_nucl(i, k)
385 ENDDO
386 ENDDO
387 ENDIF
388 ENDDO
389
390 ! Flux lessivage total
391
392 DO it = 1, nqmx - 2
393 DO k = 1, llm
394 DO i = 1, klon
395 flestottr(i, k, it) = flestottr(i, k, it) &
396 - (d_tr_lessi_nucl(i, k, it) + d_tr_lessi_impa(i, k, it)) &
397 * (paprs(i, k)-paprs(i, k+1)) / (RG * pdtphys)
398 ENDDO
399 ENDDO
400 ENDDO
401 ENDIF
402
403 ! Ecriture des sorties
404 call write_histrac(lessivage, itap, nid_tra)
405
406 if (lafin) then
407 call nf95_inq_varid(ncid_restartphy, "trs", varid)
408 call nf95_put_var(ncid_restartphy, varid, trs(:, 1))
409 endif
410
411 contains
412
413 subroutine write_histrac(lessivage, itap, nid_tra)
414
415 ! From phylmd/write_histrac.h, version 1.9 2006/02/21 08:08:30
416
417 use dimens_m, only: iim, jjm, llm
418 use histsync_m, only: histsync
419 use histwrite_m, only: histwrite
420 use temps, only: itau_phy
421 use iniadvtrac_m, only: tname
422 use comgeomphy, only: airephy
423 use dimphy, only: klon
424 use grid_change, only: gr_phy_write_2d
425 use gr_phy_write_3d_m, only: gr_phy_write_3d
426
427 logical, intent(in):: lessivage
428 integer, intent(in):: itap ! number of calls to "physiq"
429 integer, intent(in):: nid_tra
430
431 ! Variables local to the procedure:
432 integer it
433 integer itau_w ! pas de temps ecriture
434 logical, parameter:: ok_sync = .true.
435
436 !-----------------------------------------------------
437
438 itau_w = itau_phy + itap
439
440 CALL histwrite(nid_tra, "phis", itau_w, gr_phy_write_2d(pphis))
441 CALL histwrite(nid_tra, "aire", itau_w, gr_phy_write_2d(airephy))
442 CALL histwrite(nid_tra, "zmasse", itau_w, gr_phy_write_3d(zmasse))
443
444 DO it=1, nqmx - 2
445 CALL histwrite(nid_tra, tname(it+2), itau_w, &
446 gr_phy_write_3d(tr_seri(:, :, it)))
447 if (lessivage) THEN
448 CALL histwrite(nid_tra, "fl"//tname(it+2), itau_w, &
449 gr_phy_write_3d(flestottr(:, :, it)))
450 endif
451 CALL histwrite(nid_tra, "d_tr_th_"//tname(it+2), itau_w, &
452 gr_phy_write_3d(d_tr_th(:, :, it)))
453 CALL histwrite(nid_tra, "d_tr_cv_"//tname(it+2), itau_w, &
454 gr_phy_write_3d(d_tr_cv(:, :, it)))
455 CALL histwrite(nid_tra, "d_tr_cl_"//tname(it+2), itau_w, &
456 gr_phy_write_3d(d_tr_cl(:, :, it)))
457 ENDDO
458
459 CALL histwrite(nid_tra, "pplay", itau_w, gr_phy_write_3d(pplay))
460 CALL histwrite(nid_tra, "T", itau_w, gr_phy_write_3d(t_seri))
461
462 if (ok_sync) then
463 call histsync(nid_tra)
464 endif
465
466 end subroutine write_histrac
467
468 END SUBROUTINE phytrac
469
470 end module phytrac_m

  ViewVC Help
Powered by ViewVC 1.1.21