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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 35 - (show annotations)
Tue Jun 8 15:37:21 2010 UTC (13 years, 11 months ago) by guez
File size: 17211 byte(s)
Created intermediary variable for meridional wind in "calfis". Removed
unused variables.

Removed argument "firstcal" of "physiq", made it a local
variable. Removed unused argument "v" of "phytrac".

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

  ViewVC Help
Powered by ViewVC 1.1.21