/[lmdze]/trunk/phylmd/phytrac.f90
ViewVC logotype

Contents of /trunk/phylmd/phytrac.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 17 - (show annotations)
Tue Aug 5 13:31:32 2008 UTC (15 years, 9 months ago) by guez
Original Path: trunk/libf/phylmd/phytrac.f90
File size: 17956 byte(s)
Created rule for "compare_sampl_*" files in
"Documentation/Manuel_LMDZE.texfol/Graphiques/GNUmakefile".

Extracted "qcheck", "radiornpb", "minmaxqfi" into separate files.

Read pressure coordinate of ozone coefficients once per run instead of
every day.

Added some "intent" attributes.

Added argument "nq" to "ini_histday". Replaced calls to "gr_fi_ecrit"
by calls to "gr_phy_write_2d". "Sigma_O3_Royer" is written to
"histday.nc" only if "nq >= 4". Moved "ini_histrac" to module
"ini_hist".

Compute "zmasse" in "physiq", pass it to "phytrac".

Removed computations of "pftsol*" and "ppsrf*" in "phytrac".

Do not use variable "rg" from module "YOMCST" in "TLIFT".

1 module phytrac_m
2
3 ! This module is clean: no C preprocessor directive, no include line.
4
5 IMPLICIT none
6
7 private
8 public phytrac
9
10 contains
11
12 SUBROUTINE phytrac(rnpb, itap, lmt_pas, julien, gmtime, firstcal, lafin, &
13 nqmax, pdtphys, u, v, t_seri, paprs, pplay, pmfu, pmfd, pen_u, &
14 pde_u, pen_d, pde_d, coefh, fm_therm, entr_therm, yu1, yv1, ftsol, &
15 pctsrf, frac_impa, frac_nucl, presnivs, pphis, &
16 pphi, albsol, rh, cldfra, rneb, diafra, cldliq, itop_con, &
17 ibas_con, pmflxr, pmflxs, prfl, psfl, da, phi, mp, upwd, dnwd, &
18 tr_seri, zmasse)
19
20 ! From phylmd/phytrac.F, version 1.15 2006/02/21 08:08:30
21
22 ! Authors : Frédéric Hourdin, Abderrahmane Idelkadi, Marie-Alice
23 ! Foujols, Olivia
24 ! Objet : moniteur général des tendances des traceurs
25
26 ! Remarques :
27 ! 1/ L'appel de "phytrac" se fait avec "nq-2" donc nous avons bien
28 ! les vrais traceurs (en nombre "nbtr", sans la vapeur d'eau ni l'eau
29 ! liquide) dans "phytrac".
30 ! 2/ Le choix du radon et du plomb se fait juste avec un "data"
31 ! (peu propre).
32 ! Pourrait-on avoir une variable qui indiquerait le type de traceur ?
33
34 use dimens_m, only: llm
35 use indicesol, only: nbsrf
36 use dimphy, only: klon, nbtr
37 use clesphys, only: ecrit_tra
38 use clesphys2, only: iflag_con
39 use abort_gcm_m, only: abort_gcm
40 use YOMCST, only: rg
41 use ctherm, only: iflag_thermals
42 use regr_pr_comb_coefoz_m, only: regr_pr_comb_coefoz
43 use phyetat0_m, only: rlat
44 use o3_chem_m, only: o3_chem
45 use ini_hist, only: ini_histrac
46 use radiornpb_m, only: radiornpb
47 use minmaxqfi_m, only: minmaxqfi
48 use numer_rec, only: assert
49 use press_coefoz_m, only: press_coefoz
50
51 ! Arguments:
52
53 ! EN ENTREE:
54
55 ! divers:
56
57 logical, intent(in):: rnpb
58
59 integer, intent(in):: nqmax
60 ! (nombre de traceurs auxquels on applique la physique)
61
62 integer, intent(in):: itap ! number of calls to "physiq"
63 integer, intent(in):: lmt_pas ! number of time steps of "physics" per day
64 integer, intent(in):: julien !jour julien, 1 <= julien <= 360
65 integer itop_con(klon)
66 integer ibas_con(klon)
67 real, intent(in):: gmtime ! heure de la journée en fraction de jour
68 real, intent(in):: pdtphys ! pas d'integration pour la physique (s)
69 real, intent(in):: t_seri(klon, llm) ! temperature, in K
70
71 real, intent(inout):: tr_seri(klon, llm, nbtr)
72 ! (mass fractions of tracers, excluding water, at mid-layers)
73
74 real u(klon, llm)
75 real v(klon, llm)
76 real rh(klon, llm) ! humidite relative
77 real cldliq(klon, llm) ! eau liquide nuageuse
78 real cldfra(klon, llm) ! fraction nuageuse (tous les nuages)
79
80 real diafra(klon, llm)
81 ! (fraction nuageuse (convection ou stratus artificiels))
82
83 real rneb(klon, llm) ! fraction nuageuse (grande echelle)
84 real albsol(klon) ! albedo surface
85
86 real, intent(in):: paprs(klon, llm+1)
87 ! (pression pour chaque inter-couche, en Pa)
88
89 real, intent(in):: pplay(klon, llm)
90 ! (pression pour le mileu de chaque couche, en Pa)
91
92 real pphi(klon, llm) ! geopotentiel
93 real pphis(klon)
94 REAL, intent(in):: presnivs(llm)
95 logical, intent(in):: firstcal ! first call to "calfis"
96 logical, intent(in):: lafin ! fin de la physique
97
98 REAL pmflxr(klon, llm+1), pmflxs(klon, llm+1) !--lessivage convection
99 REAL prfl(klon, llm+1), psfl(klon, llm+1) !--lessivage large-scale
100
101 ! convection:
102 REAL pmfu(klon, llm) ! flux de masse dans le panache montant
103 REAL pmfd(klon, llm) ! flux de masse dans le panache descendant
104 REAL pen_u(klon, llm) ! flux entraine dans le panache montant
105
106 ! thermiques:
107
108 real fm_therm(klon, llm+1), entr_therm(klon, llm)
109
110 REAL pde_u(klon, llm) ! flux detraine dans le panache montant
111 REAL pen_d(klon, llm) ! flux entraine dans le panache descendant
112 REAL pde_d(klon, llm) ! flux detraine dans le panache descendant
113 ! KE
114 real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
115 REAL upwd(klon, llm) ! saturated updraft mass flux
116 REAL dnwd(klon, llm) ! saturated downdraft mass flux
117
118 ! Couche limite:
119
120 REAL coefh(klon, llm) ! coeff melange CL
121 REAL yu1(klon) ! vents au premier niveau
122 REAL yv1(klon) ! vents au premier niveau
123
124 ! Lessivage:
125
126 ! pour le ON-LINE
127
128 REAL frac_impa(klon, llm) ! fraction d'aerosols impactes
129 REAL frac_nucl(klon, llm) ! fraction d'aerosols nuclees
130
131 ! Arguments necessaires pour les sources et puits de traceur:
132
133 real ftsol(klon, nbsrf) ! Temperature du sol (surf)(Kelvin)
134 real pctsrf(klon, nbsrf) ! Pourcentage de sol f(nature du sol)
135
136 real, intent(in):: zmasse(:, :) ! (klon, llm)
137 ! (column-density of mass of air in a cell, in kg m-2)
138
139 ! Variables local to the procedure:
140
141 integer nsplit
142
143 ! TRACEURS
144
145 ! Sources et puits des traceurs:
146
147 ! Pour l'instant seuls les cas du rn et du pb ont ete envisages.
148
149 REAL source(klon) ! a voir lorsque le flux est prescrit
150 !
151 ! Pour la source de radon et son reservoir de sol
152
153 REAL, save:: trs(klon, nbtr) ! Concentration de radon dans le sol
154
155 REAL masktr(klon, nbtr) ! Masque reservoir de sol traceur
156 ! Masque de l'echange avec la surface
157 ! (1 = reservoir) ou (possible => 1 )
158 SAVE masktr
159 REAL fshtr(klon, nbtr) ! Flux surfacique dans le reservoir de sol
160 SAVE fshtr
161 REAL hsoltr(nbtr) ! Epaisseur equivalente du reservoir de sol
162 SAVE hsoltr
163 REAL tautr(nbtr) ! Constante de decroissance radioactive
164 SAVE tautr
165 REAL vdeptr(nbtr) ! Vitesse de depot sec dans la couche Brownienne
166 SAVE vdeptr
167 REAL scavtr(nbtr) ! Coefficient de lessivage
168 SAVE scavtr
169
170 CHARACTER itn
171 INTEGER, save:: nid_tra
172
173 ! nature du traceur
174
175 logical aerosol(nbtr) ! Nature du traceur
176 ! ! aerosol(it) = true => aerosol
177 ! ! aerosol(it) = false => gaz
178 logical clsol(nbtr) ! couche limite sol calculée
179 logical radio(nbtr) ! décroisssance radioactive
180 save aerosol, clsol, radio
181
182 ! convection tiedtke
183 INTEGER i, k, it
184 REAL delp(klon, llm)
185
186 ! Variables liees a l'ecriture de la bande histoire physique
187
188 ! Variables locales pour effectuer les appels en serie
189
190 REAL d_tr(klon, llm), d_trs(klon) ! tendances de traceurs
191 REAL d_tr_cl(klon, llm, nbtr) ! tendance de traceurs couche limite
192 REAL d_tr_cv(klon, llm, nbtr) ! tendance de traceurs conv pour chq traceur
193 REAL d_tr_th(klon, llm, nbtr) ! la tendance des thermiques
194 REAL d_tr_dec(klon, llm, 2) ! la tendance de la decroissance
195 ! ! radioactive du rn - > pb
196 REAL d_tr_lessi_impa(klon, llm, nbtr) ! la tendance du lessivage
197 ! ! par impaction
198 REAL d_tr_lessi_nucl(klon, llm, nbtr) ! la tendance du lessivage
199 ! ! par nucleation
200 REAL flestottr(klon, llm, nbtr) ! flux de lessivage
201 ! ! dans chaque couche
202
203 real ztra_th(klon, llm)
204 integer isplit
205
206 ! Controls:
207 logical:: couchelimite = .true.
208 logical:: convection = .true.
209 logical:: lessivage = .true.
210 logical, save:: inirnpb
211
212 !--------------------------------------
213
214 call assert(shape(zmasse) == (/klon, llm/), "phytrac")
215
216 if (firstcal) then
217 print *, 'phytrac: pdtphys = ', pdtphys
218 PRINT *, 'Fréquence de sortie des traceurs : ecrit_tra = ', ecrit_tra
219 if (nbtr < nqmax) call abort_gcm('phytrac', 'See above', 1)
220 inirnpb=rnpb
221
222 ! Initialisation des sorties :
223 call ini_histrac(nid_tra, pdtphys, presnivs, nqmax, lessivage)
224
225 ! Initialisation de certaines variables pour le radon et le plomb
226 ! Initialisation du traceur dans le sol (couche limite radonique)
227 trs(:, :) = 0.
228
229 open (unit=99, file='starttrac', status='old', err=999, &
230 form='formatted')
231 read(unit=99, fmt=*) (trs(i, 1), i=1, klon)
232 999 continue
233 close(unit=99)
234
235 ! Initialisation de la fraction d'aerosols lessivee
236
237 d_tr_lessi_impa(:, :, :) = 0.
238 d_tr_lessi_nucl(:, :, :) = 0.
239
240 ! Initialisation de la nature des traceurs
241
242 DO it = 1, nqmax
243 aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
244 radio(it) = .FALSE. ! par défaut pas de passage par "radiornpb"
245 clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
246 ENDDO
247
248 if (nqmax >= 3) then
249 call press_coefoz ! read input pressure levels for ozone coefficients
250 end if
251 ENDIF
252
253 ! Initialisation du traceur dans le sol (couche limite radonique)
254 if (inirnpb) THEN
255
256 radio(1)= .true.
257 radio(2)= .true.
258 clsol(1)= .true.
259 clsol(2)= .true.
260 aerosol(2) = .TRUE. ! le Pb est un aerosol
261
262 call initrrnpb(ftsol, pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, &
263 scavtr)
264 inirnpb=.false.
265 endif
266
267 ! Calcul de l'effet de la convection
268
269 if (convection) then
270 DO it=1, nqmax
271 if (iflag_con.eq.2) then
272 ! tiedke
273 CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
274 paprs, tr_seri(1, 1, it), d_tr_cv(1, 1, it))
275 else if (iflag_con.eq.3) then
276 ! KE
277 call cvltr(pdtphys, da, phi, mp, paprs, &
278 tr_seri(1, 1, it), upwd, dnwd, d_tr_cv(1, 1, it))
279 endif
280
281 DO k = 1, llm
282 DO i = 1, klon
283 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cv(i, k, it)
284 ENDDO
285 ENDDO
286 WRITE(unit=itn, fmt='(i1)') it
287 CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, &
288 'convection, tracer index = ' // itn)
289 ENDDO
290 endif
291
292 ! Calcul de l'effet des thermiques
293
294 do it=1, nqmax
295 do k=1, llm
296 do i=1, klon
297 d_tr_th(i, k, it)=0.
298 tr_seri(i, k, it)=max(tr_seri(i, k, it), 0.)
299 tr_seri(i, k, it)=min(tr_seri(i, k, it), 1.e10)
300 enddo
301 enddo
302 enddo
303
304 if (iflag_thermals > 0) then
305 nsplit=10
306 DO it=1, nqmax
307 do isplit=1, nsplit
308 ! Thermiques
309 call dqthermcell(klon, llm, pdtphys/nsplit &
310 , fm_therm, entr_therm, zmasse &
311 , tr_seri(1:klon, 1:llm, it), d_tr, ztra_th)
312
313 do k=1, llm
314 do i=1, klon
315 d_tr(i, k)=pdtphys*d_tr(i, k)/nsplit
316 d_tr_th(i, k, it)=d_tr_th(i, k, it)+d_tr(i, k)
317 tr_seri(i, k, it)=max(tr_seri(i, k, it)+d_tr(i, k), 0.)
318 enddo
319 enddo
320 enddo
321 ENDDO
322 endif
323
324 ! Calcul de l'effet de la couche limite
325
326 if (couchelimite) then
327
328 DO k = 1, llm
329 DO i = 1, klon
330 delp(i, k) = paprs(i, k)-paprs(i, k+1)
331 ENDDO
332 ENDDO
333
334 ! MAF modif pour tenir compte du cas rnpb + traceur
335 DO it=1, nqmax
336 if (clsol(it)) then
337 ! couche limite avec quantite dans le sol calculee
338 CALL cltracrn(it, pdtphys, yu1, yv1, &
339 coefh, t_seri, ftsol, pctsrf, &
340 tr_seri(1, 1, it), trs(1, it), &
341 paprs, pplay, delp, &
342 masktr(1, it), fshtr(1, it), hsoltr(it), &
343 tautr(it), vdeptr(it), &
344 rlat, &
345 d_tr_cl(1, 1, it), d_trs)
346 DO k = 1, llm
347 DO i = 1, klon
348 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
349 ENDDO
350 ENDDO
351
352 ! Traceur ds sol
353
354 DO i = 1, klon
355 trs(i, it) = trs(i, it) + d_trs(i)
356 END DO
357 else ! couche limite avec flux prescrit
358 !MAF provisoire source / traceur a creer
359 DO i=1, klon
360 source(i) = 0.0 ! pas de source, pour l'instant
361 ENDDO
362
363 CALL cltrac(pdtphys, coefh, t_seri, &
364 tr_seri(1, 1, it), source, &
365 paprs, pplay, delp, &
366 d_tr_cl(1, 1, it))
367 DO k = 1, llm
368 DO i = 1, klon
369 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
370 ENDDO
371 ENDDO
372 endif
373 ENDDO
374
375 endif ! couche limite
376
377 ! Calcul de l'effet du puits radioactif
378
379 ! MAF il faudrait faire une modification pour passer dans radiornpb
380 ! si radio=true mais pour l'instant radiornpb propre au cas rnpb
381 if (rnpb) then
382 d_tr_dec(:, :, :) = radiornpb(tr_seri, pdtphys, tautr)
383 DO it=1, nqmax
384 if (radio(it)) then
385 tr_seri(:, :, it) = tr_seri(:, :, it) + d_tr_dec(:, :, it)
386 WRITE(unit=itn, fmt='(i1)') it
387 CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, 'puits rn it='//itn)
388 endif
389 ENDDO
390 endif ! rnpb decroissance radioactive
391
392 if (nqmax >= 3) then
393 ! Ozone as a tracer:
394 if (mod(itap - 1, lmt_pas) == 0) then
395 ! Once per day, update the coefficients for ozone chemistry:
396 call regr_pr_comb_coefoz(julien)
397 end if
398 call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
399 end if
400
401 ! Calcul de l'effet de la precipitation
402
403 IF (lessivage) THEN
404 d_tr_lessi_nucl(:, :, :) = 0.
405 d_tr_lessi_impa(:, :, :) = 0.
406 flestottr(:, :, :) = 0.
407
408 ! tendance des aerosols nuclees et impactes
409
410 DO it = 1, nqmax
411 IF (aerosol(it)) THEN
412 DO k = 1, llm
413 DO i = 1, klon
414 d_tr_lessi_nucl(i, k, it) = d_tr_lessi_nucl(i, k, it) + &
415 ( 1 - frac_nucl(i, k) )*tr_seri(i, k, it)
416 d_tr_lessi_impa(i, k, it) = d_tr_lessi_impa(i, k, it) + &
417 ( 1 - frac_impa(i, k) )*tr_seri(i, k, it)
418 ENDDO
419 ENDDO
420 ENDIF
421 ENDDO
422
423 ! Mises a jour des traceurs + calcul des flux de lessivage
424 ! Mise a jour due a l'impaction et a la nucleation
425
426 DO it = 1, nqmax
427 IF (aerosol(it)) THEN
428 DO k = 1, llm
429 DO i = 1, klon
430 tr_seri(i, k, it)=tr_seri(i, k, it) &
431 *frac_impa(i, k)*frac_nucl(i, k)
432 ENDDO
433 ENDDO
434 ENDIF
435 ENDDO
436
437 ! Flux lessivage total
438
439 DO it = 1, nqmax
440 DO k = 1, llm
441 DO i = 1, klon
442 flestottr(i, k, it) = flestottr(i, k, it) - &
443 ( d_tr_lessi_nucl(i, k, it) + &
444 d_tr_lessi_impa(i, k, it) ) * &
445 ( paprs(i, k)-paprs(i, k+1) ) / &
446 (RG * pdtphys)
447 ENDDO
448 ENDDO
449 ENDDO
450 ENDIF
451
452 ! Ecriture des sorties
453 call write_histrac(lessivage, nqmax, itap, nid_tra)
454
455 if (lafin) then
456 print *, "C'est la fin de la physique."
457 open (unit=99, file='restarttrac', form='formatted')
458 do i=1, klon
459 write(unit=99, fmt=*) trs(i, 1)
460 enddo
461 PRINT *, 'Ecriture du fichier restarttrac'
462 close(99)
463 endif
464
465 contains
466
467 subroutine write_histrac(lessivage, nqmax, itap, nid_tra)
468
469 ! From phylmd/write_histrac.h, version 1.9 2006/02/21 08:08:30
470
471 use dimens_m, only: iim, jjm, llm
472 use ioipsl, only: histwrite, histsync
473 use temps, only: itau_phy
474 use advtrac_m, only: tnom
475 use comgeomphy, only: airephy
476 use dimphy, only: klon
477
478 logical, intent(in):: lessivage
479
480 integer, intent(in):: nqmax
481 ! (nombre de traceurs auxquels on applique la physique)
482
483 integer, intent(in):: itap ! number of calls to "physiq"
484 integer, intent(in):: nid_tra
485
486 ! Variables local to the procedure:
487 integer it
488 integer itau_w ! pas de temps ecriture
489 REAL zx_tmp_2d(iim, jjm+1), zx_tmp_3d(iim, jjm+1, llm)
490 logical, parameter:: ok_sync = .true.
491
492 !-----------------------------------------------------
493
494 itau_w = itau_phy + itap
495
496 CALL gr_fi_ecrit(1, klon, iim, jjm+1, pphis, zx_tmp_2d)
497 CALL histwrite(nid_tra, "phis", itau_w, zx_tmp_2d)
498
499 CALL gr_fi_ecrit(1, klon, iim, jjm+1, airephy, zx_tmp_2d)
500 CALL histwrite(nid_tra, "aire", itau_w, zx_tmp_2d)
501
502 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, zmasse, zx_tmp_3d)
503 CALL histwrite(nid_tra, "zmasse", itau_w, zx_tmp_3d)
504
505 DO it=1, nqmax
506 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, tr_seri(1, 1, it), zx_tmp_3d)
507 CALL histwrite(nid_tra, tnom(it+2), itau_w, zx_tmp_3d)
508 if (lessivage) THEN
509 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, flestottr(1, 1, it), &
510 zx_tmp_3d)
511 CALL histwrite(nid_tra, "fl"//tnom(it+2), itau_w, zx_tmp_3d)
512 endif
513
514 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_th(1, 1, it), zx_tmp_3d)
515 CALL histwrite(nid_tra, "d_tr_th_"//tnom(it+2), itau_w, zx_tmp_3d)
516 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_cv(1, 1, it), zx_tmp_3d)
517 CALL histwrite(nid_tra, "d_tr_cv_"//tnom(it+2), itau_w, zx_tmp_3d)
518 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_cl(1, 1, it), zx_tmp_3d)
519 CALL histwrite(nid_tra, "d_tr_cl_"//tnom(it+2), itau_w, zx_tmp_3d)
520 ENDDO
521
522 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pplay, zx_tmp_3d)
523 CALL histwrite(nid_tra, "pplay", itau_w, zx_tmp_3d)
524
525 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, t_seri, zx_tmp_3d)
526 CALL histwrite(nid_tra, "t", itau_w, zx_tmp_3d)
527
528 if (ok_sync) then
529 call histsync(nid_tra)
530 endif
531
532 end subroutine write_histrac
533
534 END SUBROUTINE phytrac
535
536 end module phytrac_m

  ViewVC Help
Powered by ViewVC 1.1.21