/[lmdze]/trunk/Sources/phylmd/phytrac.f
ViewVC logotype

Contents of /trunk/Sources/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21