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

Contents of /trunk/libf/phylmd/phytrac.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 34 - (show annotations)
Wed Jun 2 11:01:12 2010 UTC (13 years, 11 months ago) by guez
File size: 17236 byte(s)
Split "ini_hist.f90" into single-procedure files.

In "calfis" and "physiq", removed dummy argument "nq" since "nq" must
be equal to "nqmx".

In "calfis", renamed dummy argument "pq" to "q", same name as actual
argument in "leapfrog". Renamed local variable "zqfi" to "qx", same
name as dummy argument in "physiq".

Removed arguments "itop_con" and "ibas_con" of "phytrac", which were
not used.

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

  ViewVC Help
Powered by ViewVC 1.1.21