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

Contents of /trunk/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 72 - (show annotations)
Tue Jul 23 13:00:07 2013 UTC (10 years, 9 months ago) by guez
Original Path: trunk/libf/phylmd/phytrac.f90
File size: 17453 byte(s)
NaN to signalling NaN in gfortran_debug.mk.

Removed unused procedures in getincom and getincom2. In procedure
conf_interface, replaced call to getincom by new namelist. Moved
procedure conf_interface into module interface_surf.

Added variables sig1 and w01 to startphy.nc and restartphy.nc, for
procedure cv_driver. Renamed (ema_)?work1 and (ema_)?work2 to sig1 and
w01 in concvl and physiq.

Deleted unused arguments of clmain, clqh and intersurf_hq, among which
(y)?sollwdown. Following LMDZ, in physiq, read sollw instead of
sollwdown from startphy.nc, write sollw instead of sollwdown to
restartphy.nc.

In procedure sw, initialized zfs[ud][pn]a[di], for runs where ok_ade
and ok_aie are false. (Following LMDZ.)

Added dimension klev to startphy.nc and restartphy.nc, and deleted
dimension horizon_vertical. Made t_ancien and q_ancien two-dimensional
NetCDF variables. Bug fix: in phyetat0, define ratqs, clwcon and
rnebcon for vertical levels >=2.

Bug fix: set mfg, p[de]n_[ud] to 0. when iflag_con >= 3. (Following LMDZ.)

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

  ViewVC Help
Powered by ViewVC 1.1.21