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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 20 - (show annotations)
Wed Oct 15 16:19:57 2008 UTC (15 years, 6 months ago) by guez
File size: 17614 byte(s)
Deleted argument "presnivs" of "physiq", "ini_histhf", "ini_histhf3d",
"ini_histday", "ini_histins", "ini_histrac", "phytrac". Access it from
"comvert" instead.

Replaced calls to NetCDF Fortran 77 interface by calls to Fortran 90
interface or to NetCDF95.

Procedure "gr_phy_write_3d" now works with a variable of arbitrary
size in the second dimension.

Annotated use statements with "only" clause.

Replaced calls to NetCDF interface version 2 by calls to Fortran 90
interface in "guide.f90" and "read_reanalyse.f".

In "write_histrac", replaced calls to "gr_fi_ecrit" by calls to
"gr_phy_write_2d" and "gr_phy_write_3d".

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

  ViewVC Help
Powered by ViewVC 1.1.21