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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 120 - (show annotations)
Tue Jan 13 14:56:15 2015 UTC (9 years, 4 months ago) by guez
Original Path: trunk/phylmd/phytrac.f
File size: 15909 byte(s)
In procedure fxhyp, removed the possibility to set scal180 to
false. The useful lower bound of fhyp and xxpr is not 0. It does not
make sense to give the save attribute to is2 since fxhyp is only
called one per run. Bug fix: is2 could be used without being
defined. The bug did not appear because is2 had the save attribute so
it was initialized at 0.

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

  ViewVC Help
Powered by ViewVC 1.1.21