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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 23 - (show annotations)
Mon Dec 14 15:25:16 2009 UTC (14 years, 5 months ago) by guez
File size: 17553 byte(s)
Split "orografi.f": one file for each procedure. Put the created files
in new directory "Orography".

Removed argument "vcov" of procedure "sortvarc". Removed arguments
"itau" and "time" of procedure "caldyn0". Removed arguments "itau",
"time" and "vcov" of procedure "sortvarc0".

Removed argument "time" of procedure "dynredem1". Removed NetCDF
variable "temps" in files "start.nc" and "restart.nc", because its
value is always 0.

Removed argument "nq" of procedures "iniadvtrac" and "leapfrog". The
number of "tracers read in "traceur.def" must now be equal to "nqmx",
or "nqmx" must equal 4 if there is no file "traceur.def". Replaced
variable "nq" by constant "nqmx" in "leapfrog".

NetCDF variable for ozone field in "coefoz.nc" must now be called
"tro3" instead of "r".

Fixed bug in "zenang".

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 logical, parameter:: ok_sync = .true.
490
491 !-----------------------------------------------------
492
493 itau_w = itau_phy + itap
494
495 CALL histwrite(nid_tra, "phis", itau_w, gr_phy_write_2d(pphis))
496 CALL histwrite(nid_tra, "aire", itau_w, gr_phy_write_2d(airephy))
497 CALL histwrite(nid_tra, "zmasse", itau_w, gr_phy_write_3d(zmasse))
498
499 DO it=1, nq_phys
500 CALL histwrite(nid_tra, tnom(it+2), itau_w, &
501 gr_phy_write_3d(tr_seri(:, :, it)))
502 if (lessivage) THEN
503 CALL histwrite(nid_tra, "fl"//tnom(it+2), itau_w, &
504 gr_phy_write_3d(flestottr(:, :, it)))
505 endif
506 CALL histwrite(nid_tra, "d_tr_th_"//tnom(it+2), itau_w, &
507 gr_phy_write_3d(d_tr_th(:, :, it)))
508 CALL histwrite(nid_tra, "d_tr_cv_"//tnom(it+2), itau_w, &
509 gr_phy_write_3d(d_tr_cv(:, :, it)))
510 CALL histwrite(nid_tra, "d_tr_cl_"//tnom(it+2), itau_w, &
511 gr_phy_write_3d(d_tr_cl(:, :, it)))
512 ENDDO
513
514 CALL histwrite(nid_tra, "pplay", itau_w, gr_phy_write_3d(pplay))
515 CALL histwrite(nid_tra, "T", itau_w, gr_phy_write_3d(t_seri))
516
517 if (ok_sync) then
518 call histsync(nid_tra)
519 endif
520
521 end subroutine write_histrac
522
523 END SUBROUTINE phytrac
524
525 end module phytrac_m

  ViewVC Help
Powered by ViewVC 1.1.21