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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 18 - (show annotations)
Thu Aug 7 12:29:13 2008 UTC (15 years, 9 months ago) by guez
File size: 18094 byte(s)
In module "regr_pr", rewrote scanning of horizontal positions as a
single set of loops, using a mask.

Added some "intent" attributes.

In "dynredem0", replaced calls to Fortran 77 interface of NetCDF by
calls to NetCDF95. Removed calls to "nf_redef", regrouped all writing
operations. In "dynredem1", replaced some calls to Fortran 77
interface of NetCDF by calls to Fortran 90 interface.

Renamed variable "nqmax" to "nq_phys".

In "physiq", if "nq >= 5" then "wo" is computed from the
parameterization of "Cariolle".

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, 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):: nq_phys
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 zmasse")
215 call assert(shape(tr_seri) == (/klon, llm, nbtr/), "phytrac tr_seri")
216
217 if (firstcal) then
218 print *, 'phytrac: pdtphys = ', pdtphys
219 PRINT *, 'Fréquence de sortie des traceurs : ecrit_tra = ', ecrit_tra
220 if (nbtr < nq_phys) call abort_gcm('phytrac', 'nbtr < nq_phys', 1)
221 inirnpb=rnpb
222
223 ! Initialisation des sorties :
224 call ini_histrac(nid_tra, pdtphys, presnivs, nq_phys, lessivage)
225
226 ! Initialisation de certaines variables pour le radon et le plomb
227 ! Initialisation du traceur dans le sol (couche limite radonique)
228 trs(:, :) = 0.
229
230 open (unit=99, file='starttrac', status='old', err=999, &
231 form='formatted')
232 read(unit=99, fmt=*) (trs(i, 1), i=1, klon)
233 999 continue
234 close(unit=99)
235
236 ! Initialisation de la fraction d'aerosols lessivee
237
238 d_tr_lessi_impa(:, :, :) = 0.
239 d_tr_lessi_nucl(:, :, :) = 0.
240
241 ! Initialisation de la nature des traceurs
242
243 DO it = 1, nq_phys
244 aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
245 radio(it) = .FALSE. ! par défaut pas de passage par "radiornpb"
246 clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
247 ENDDO
248
249 if (nq_phys >= 3) then
250 call press_coefoz ! read input pressure levels for ozone coefficients
251 end if
252 ENDIF
253
254 ! Initialisation du traceur dans le sol (couche limite radonique)
255 if (inirnpb) THEN
256
257 radio(1)= .true.
258 radio(2)= .true.
259 clsol(1)= .true.
260 clsol(2)= .true.
261 aerosol(2) = .TRUE. ! le Pb est un aerosol
262
263 call initrrnpb(ftsol, pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, &
264 scavtr)
265 inirnpb=.false.
266 endif
267
268 ! Calcul de l'effet de la convection
269
270 if (convection) then
271 DO it=1, nq_phys
272 if (iflag_con.eq.2) then
273 ! tiedke
274 CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
275 paprs, tr_seri(1, 1, it), d_tr_cv(1, 1, it))
276 else if (iflag_con.eq.3) then
277 ! KE
278 call cvltr(pdtphys, da, phi, mp, paprs, &
279 tr_seri(1, 1, it), upwd, dnwd, d_tr_cv(1, 1, it))
280 endif
281
282 DO k = 1, llm
283 DO i = 1, klon
284 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cv(i, k, it)
285 ENDDO
286 ENDDO
287 WRITE(unit=itn, fmt='(i1)') it
288 CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, &
289 'convection, tracer index = ' // itn)
290 ENDDO
291 endif
292
293 ! Calcul de l'effet des thermiques
294
295 do it=1, nq_phys
296 do k=1, llm
297 do i=1, klon
298 d_tr_th(i, k, it)=0.
299 tr_seri(i, k, it)=max(tr_seri(i, k, it), 0.)
300 tr_seri(i, k, it)=min(tr_seri(i, k, it), 1.e10)
301 enddo
302 enddo
303 enddo
304
305 if (iflag_thermals > 0) then
306 nsplit=10
307 DO it=1, nq_phys
308 do isplit=1, nsplit
309 ! Thermiques
310 call dqthermcell(klon, llm, pdtphys/nsplit &
311 , fm_therm, entr_therm, zmasse &
312 , tr_seri(1:klon, 1:llm, it), d_tr, ztra_th)
313
314 do k=1, llm
315 do i=1, klon
316 d_tr(i, k)=pdtphys*d_tr(i, k)/nsplit
317 d_tr_th(i, k, it)=d_tr_th(i, k, it)+d_tr(i, k)
318 tr_seri(i, k, it)=max(tr_seri(i, k, it)+d_tr(i, k), 0.)
319 enddo
320 enddo
321 enddo
322 ENDDO
323 endif
324
325 ! Calcul de l'effet de la couche limite
326
327 if (couchelimite) then
328
329 DO k = 1, llm
330 DO i = 1, klon
331 delp(i, k) = paprs(i, k)-paprs(i, k+1)
332 ENDDO
333 ENDDO
334
335 ! MAF modif pour tenir compte du cas rnpb + traceur
336 DO it=1, nq_phys
337 if (clsol(it)) then
338 ! couche limite avec quantite dans le sol calculee
339 CALL cltracrn(it, pdtphys, yu1, yv1, &
340 coefh, t_seri, ftsol, pctsrf, &
341 tr_seri(1, 1, it), trs(1, it), &
342 paprs, pplay, delp, &
343 masktr(1, it), fshtr(1, it), hsoltr(it), &
344 tautr(it), vdeptr(it), &
345 rlat, &
346 d_tr_cl(1, 1, it), d_trs)
347 DO k = 1, llm
348 DO i = 1, klon
349 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
350 ENDDO
351 ENDDO
352
353 ! Traceur ds sol
354
355 DO i = 1, klon
356 trs(i, it) = trs(i, it) + d_trs(i)
357 END DO
358 else ! couche limite avec flux prescrit
359 !MAF provisoire source / traceur a creer
360 DO i=1, klon
361 source(i) = 0.0 ! pas de source, pour l'instant
362 ENDDO
363
364 CALL cltrac(pdtphys, coefh, t_seri, &
365 tr_seri(1, 1, it), source, &
366 paprs, pplay, delp, &
367 d_tr_cl(1, 1, it))
368 DO k = 1, llm
369 DO i = 1, klon
370 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
371 ENDDO
372 ENDDO
373 endif
374 ENDDO
375
376 endif ! couche limite
377
378 ! Calcul de l'effet du puits radioactif
379
380 ! MAF il faudrait faire une modification pour passer dans radiornpb
381 ! si radio=true mais pour l'instant radiornpb propre au cas rnpb
382 if (rnpb) then
383 d_tr_dec(:, :, :) = radiornpb(tr_seri, pdtphys, tautr)
384 DO it=1, nq_phys
385 if (radio(it)) then
386 tr_seri(:, :, it) = tr_seri(:, :, it) + d_tr_dec(:, :, it)
387 WRITE(unit=itn, fmt='(i1)') it
388 CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, 'puits rn it='//itn)
389 endif
390 ENDDO
391 endif ! rnpb decroissance radioactive
392
393 if (nq_phys >= 3) then
394 ! Ozone as a tracer:
395 if (mod(itap - 1, lmt_pas) == 0) then
396 ! Once per day, update the coefficients for ozone chemistry:
397 call regr_pr_comb_coefoz(julien)
398 end if
399 call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
400 end if
401
402 ! Calcul de l'effet de la precipitation
403
404 IF (lessivage) THEN
405 d_tr_lessi_nucl(:, :, :) = 0.
406 d_tr_lessi_impa(:, :, :) = 0.
407 flestottr(:, :, :) = 0.
408
409 ! tendance des aerosols nuclees et impactes
410
411 DO it = 1, nq_phys
412 IF (aerosol(it)) THEN
413 DO k = 1, llm
414 DO i = 1, klon
415 d_tr_lessi_nucl(i, k, it) = d_tr_lessi_nucl(i, k, it) + &
416 ( 1 - frac_nucl(i, k) )*tr_seri(i, k, it)
417 d_tr_lessi_impa(i, k, it) = d_tr_lessi_impa(i, k, it) + &
418 ( 1 - frac_impa(i, k) )*tr_seri(i, k, it)
419 ENDDO
420 ENDDO
421 ENDIF
422 ENDDO
423
424 ! Mises a jour des traceurs + calcul des flux de lessivage
425 ! Mise a jour due a l'impaction et a la nucleation
426
427 DO it = 1, nq_phys
428 IF (aerosol(it)) THEN
429 DO k = 1, llm
430 DO i = 1, klon
431 tr_seri(i, k, it)=tr_seri(i, k, it) &
432 *frac_impa(i, k)*frac_nucl(i, k)
433 ENDDO
434 ENDDO
435 ENDIF
436 ENDDO
437
438 ! Flux lessivage total
439
440 DO it = 1, nq_phys
441 DO k = 1, llm
442 DO i = 1, klon
443 flestottr(i, k, it) = flestottr(i, k, it) - &
444 ( d_tr_lessi_nucl(i, k, it) + &
445 d_tr_lessi_impa(i, k, it) ) * &
446 ( paprs(i, k)-paprs(i, k+1) ) / &
447 (RG * pdtphys)
448 ENDDO
449 ENDDO
450 ENDDO
451 ENDIF
452
453 ! Ecriture des sorties
454 call write_histrac(lessivage, nq_phys, itap, nid_tra)
455
456 if (lafin) then
457 print *, "C'est la fin de la physique."
458 open (unit=99, file='restarttrac', form='formatted')
459 do i=1, klon
460 write(unit=99, fmt=*) trs(i, 1)
461 enddo
462 PRINT *, 'Ecriture du fichier restarttrac'
463 close(99)
464 endif
465
466 contains
467
468 subroutine write_histrac(lessivage, nq_phys, itap, nid_tra)
469
470 ! From phylmd/write_histrac.h, version 1.9 2006/02/21 08:08:30
471
472 use dimens_m, only: iim, jjm, llm
473 use ioipsl, only: histwrite, histsync
474 use temps, only: itau_phy
475 use iniadvtrac_m, only: tnom
476 use comgeomphy, only: airephy
477 use dimphy, only: klon
478
479 logical, intent(in):: lessivage
480
481 integer, intent(in):: nq_phys
482 ! (nombre de traceurs auxquels on applique la physique)
483
484 integer, intent(in):: itap ! number of calls to "physiq"
485 integer, intent(in):: nid_tra
486
487 ! Variables local to the procedure:
488 integer it
489 integer itau_w ! pas de temps ecriture
490 REAL zx_tmp_2d(iim, jjm+1), zx_tmp_3d(iim, jjm+1, llm)
491 logical, parameter:: ok_sync = .true.
492
493 !-----------------------------------------------------
494
495 itau_w = itau_phy + itap
496
497 CALL gr_fi_ecrit(1, klon, iim, jjm+1, pphis, zx_tmp_2d)
498 CALL histwrite(nid_tra, "phis", itau_w, zx_tmp_2d)
499
500 CALL gr_fi_ecrit(1, klon, iim, jjm+1, airephy, zx_tmp_2d)
501 CALL histwrite(nid_tra, "aire", itau_w, zx_tmp_2d)
502
503 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, zmasse, zx_tmp_3d)
504 CALL histwrite(nid_tra, "zmasse", itau_w, zx_tmp_3d)
505
506 DO it=1, nq_phys
507 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, tr_seri(1, 1, it), zx_tmp_3d)
508 CALL histwrite(nid_tra, tnom(it+2), itau_w, zx_tmp_3d)
509 if (lessivage) THEN
510 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, flestottr(1, 1, it), &
511 zx_tmp_3d)
512 CALL histwrite(nid_tra, "fl"//tnom(it+2), itau_w, zx_tmp_3d)
513 endif
514
515 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_th(1, 1, it), zx_tmp_3d)
516 CALL histwrite(nid_tra, "d_tr_th_"//tnom(it+2), itau_w, zx_tmp_3d)
517 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_cv(1, 1, it), zx_tmp_3d)
518 CALL histwrite(nid_tra, "d_tr_cv_"//tnom(it+2), itau_w, zx_tmp_3d)
519 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_cl(1, 1, it), zx_tmp_3d)
520 CALL histwrite(nid_tra, "d_tr_cl_"//tnom(it+2), itau_w, zx_tmp_3d)
521 ENDDO
522
523 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pplay, zx_tmp_3d)
524 CALL histwrite(nid_tra, "pplay", itau_w, zx_tmp_3d)
525
526 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, t_seri, zx_tmp_3d)
527 CALL histwrite(nid_tra, "t", itau_w, zx_tmp_3d)
528
529 if (ok_sync) then
530 call histsync(nid_tra)
531 endif
532
533 end subroutine write_histrac
534
535 END SUBROUTINE phytrac
536
537 end module phytrac_m

  ViewVC Help
Powered by ViewVC 1.1.21