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

Contents of /trunk/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 118 - (show annotations)
Thu Dec 18 17:30:24 2014 UTC (9 years, 4 months ago) by guez
File size: 15889 byte(s)
In file grilles_gcm.nc, renamed variable phis to orog, deleted
variable presnivs.

Removed variable bug_ozone from module clesphys.

In procedure ozonecm, moved computation of sint and cost out of the
loops on horizontal position and vertical level. Inverted the order of
the two loops. We can then move all computations from slat to aprim
out of the loop on vertical levels. Created variable slat2, following
LMDZ. Moved the limitation of column-density of ozone in cell at 1e-12
from radlwsw to ozonecm, following LMDZ.

Removed unused arguments u, albsol, rh, cldfra, rneb, diafra, cldliq,
pmflxr, pmflxs, prfl, psfl of phytrac.

In procedure yamada4, for all the arrays, replaced the dimension klon
by ngrid. At the end of the procedure, for the computation of kmn,kn,
kq and q2, changed the upper limit of the loop index from klon to ngrid.

In radlwsw, for the calculation of pozon, removed the factor
paprs(iof+i, 1)/101325, as in LMDZ. In procedure sw, removed the
factor 101325.0/PPSOL(JL), as in LMDZ.

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 cltracrn_m, only: cltracrn
33 use ctherm, only: iflag_thermals
34 use dimens_m, only: llm, nqmx
35 use dimphy, only: klon
36 use indicesol, only: nbsrf
37 use ini_histrac_m, only: ini_histrac
38 use initrrnpb_m, only: initrrnpb
39 use minmaxqfi_m, only: minmaxqfi
40 use nflxtr_m, only: nflxtr
41 use nr_util, only: assert
42 use o3_chem_m, only: o3_chem
43 use phyetat0_m, only: rlat
44 use press_coefoz_m, only: press_coefoz
45 use radiornpb_m, only: radiornpb
46 use regr_pr_comb_coefoz_m, only: regr_pr_comb_coefoz
47 use SUPHEC_M, only: rg
48
49 integer, intent(in):: itap ! number of calls to "physiq"
50 integer, intent(in):: lmt_pas ! number of time steps of "physics" per day
51 integer, intent(in):: julien !jour julien, 1 <= julien <= 360
52 real, intent(in):: gmtime ! heure de la journ\'ee en fraction de jour
53 logical, intent(in):: firstcal ! first call to "calfis"
54 logical, intent(in):: lafin ! fin de la physique
55 real, intent(in):: pdtphys ! pas d'integration pour la physique (s)
56 real, intent(in):: t_seri(klon, llm) ! temperature, in K
57
58 real, intent(in):: paprs(klon, llm+1)
59 ! (pression pour chaque inter-couche, en Pa)
60
61 real, intent(in):: pplay(klon, llm)
62 ! (pression pour le mileu de chaque couche, en Pa)
63
64 ! convection:
65
66 REAL, intent(in):: pmfu(klon, llm) ! flux de masse dans le panache montant
67
68 REAL, intent(in):: pmfd(klon, llm)
69 ! flux de masse dans le panache descendant
70
71 REAL pde_u(klon, llm) ! flux detraine dans le panache montant
72 REAL pen_d(klon, llm) ! flux entraine dans le panache descendant
73 REAL coefh(klon, llm) ! coeff melange couche limite
74
75 ! thermiques:
76 real fm_therm(klon, llm+1), entr_therm(klon, llm)
77
78 ! Couche limite:
79 REAL yu1(klon) ! vents au premier niveau
80 REAL yv1(klon) ! vents au premier niveau
81
82 ! Arguments n\'ecessaires pour les sources et puits de traceur :
83 real, intent(in):: ftsol(klon, nbsrf) ! Temperature du sol (surf)(Kelvin)
84 real pctsrf(klon, nbsrf) ! Pourcentage de sol f(nature du sol)
85
86 ! Lessivage pour le on-line
87 REAL frac_impa(klon, llm) ! fraction d'aerosols impactes
88 REAL frac_nucl(klon, llm) ! fraction d'aerosols nuclees
89
90 real, intent(in):: pphis(klon)
91
92 ! Kerry Emanuel
93 real, intent(in):: da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
94 REAL, intent(in):: upwd(klon, llm) ! saturated updraft mass flux
95 REAL, intent(in):: dnwd(klon, llm) ! saturated downdraft mass flux
96
97 real, intent(inout):: tr_seri(:, :, :) ! (klon, llm, nqmx - 2)
98 ! (mass fractions of tracers, excluding water, at mid-layers)
99
100 real, intent(in):: zmasse(:, :) ! (klon, llm)
101 ! (column-density of mass of air in a cell, in kg m-2)
102
103 ! Variables local to the procedure:
104
105 integer nsplit
106
107 ! TRACEURS
108
109 ! Sources et puits des traceurs:
110
111 ! Pour l'instant seuls les cas du rn et du pb ont ete envisages.
112
113 REAL source(klon) ! a voir lorsque le flux est prescrit
114 !
115 ! Pour la source de radon et son reservoir de sol
116
117 REAL, save:: trs(klon, nqmx - 2) ! Concentration de radon dans le sol
118
119 REAL masktr(klon, nqmx - 2) ! Masque reservoir de sol traceur
120 ! Masque de l'echange avec la surface
121 ! (1 = reservoir) ou (possible => 1)
122 SAVE masktr
123 REAL fshtr(klon, nqmx - 2) ! Flux surfacique dans le reservoir de sol
124 SAVE fshtr
125 REAL hsoltr(nqmx - 2) ! Epaisseur equivalente du reservoir de sol
126 SAVE hsoltr
127 REAL tautr(nqmx - 2) ! Constante de decroissance radioactive
128 SAVE tautr
129 REAL vdeptr(nqmx - 2) ! Vitesse de depot sec dans la couche Brownienne
130 SAVE vdeptr
131 REAL scavtr(nqmx - 2) ! Coefficient de lessivage
132 SAVE scavtr
133
134 CHARACTER itn
135 INTEGER, save:: nid_tra
136
137 ! nature du traceur
138
139 logical aerosol(nqmx - 2) ! Nature du traceur
140 ! ! aerosol(it) = true => aerosol
141 ! ! aerosol(it) = false => gaz
142 logical clsol(nqmx - 2) ! couche limite sol calcul\'ee
143 logical radio(nqmx - 2) ! d\'ecroisssance radioactive
144 save aerosol, clsol, radio
145
146 ! convection tiedtke
147 INTEGER i, k, it
148 REAL delp(klon, llm)
149
150 ! Variables liees a l'ecriture de la bande histoire physique
151
152 ! Variables locales pour effectuer les appels en serie
153
154 REAL d_tr(klon, llm), d_trs(klon) ! tendances de traceurs
155 REAL d_tr_cl(klon, llm, nqmx - 2) ! tendance de traceurs couche limite
156 REAL d_tr_cv(klon, llm, nqmx - 2) ! tendance de traceurs conv pour chq traceur
157 REAL d_tr_th(klon, llm, nqmx - 2) ! la tendance des thermiques
158 REAL d_tr_dec(klon, llm, 2) ! la tendance de la decroissance
159 ! ! radioactive du rn - > pb
160 REAL d_tr_lessi_impa(klon, llm, nqmx - 2) ! la tendance du lessivage
161 ! ! par impaction
162 REAL d_tr_lessi_nucl(klon, llm, nqmx - 2) ! la tendance du lessivage
163 ! ! par nucleation
164 REAL flestottr(klon, llm, nqmx - 2) ! flux de lessivage
165 ! ! dans chaque couche
166
167 real ztra_th(klon, llm)
168 integer isplit
169
170 ! Controls:
171 logical:: couchelimite = .true.
172 logical:: convection = .true.
173 logical:: lessivage = .true.
174 logical, save:: inirnpb
175
176 !--------------------------------------
177
178 call assert(shape(zmasse) == (/klon, llm/), "phytrac zmasse")
179 call assert(shape(tr_seri) == (/klon, llm, nqmx - 2/), "phytrac tr_seri")
180
181 if (firstcal) then
182 print *, 'phytrac: pdtphys = ', pdtphys
183 PRINT *, 'Frequency of tracer output: ecrit_tra = ', ecrit_tra
184 inirnpb = .true.
185
186 ! Initialisation des sorties :
187 call ini_histrac(nid_tra, pdtphys, nqmx - 2, lessivage)
188
189 ! Initialisation de certaines variables pour le radon et le plomb
190 ! Initialisation du traceur dans le sol (couche limite radonique)
191 trs(:, :) = 0.
192
193 open (unit=99, file='starttrac', status='old', err=999, &
194 form='formatted')
195 read(unit=99, fmt=*) (trs(i, 1), i=1, klon)
196 999 continue
197 close(unit=99)
198
199 ! Initialisation de la fraction d'aerosols lessivee
200
201 d_tr_lessi_impa = 0.
202 d_tr_lessi_nucl = 0.
203
204 ! Initialisation de la nature des traceurs
205
206 DO it = 1, nqmx - 2
207 aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
208 radio(it) = .FALSE. ! par d\'efaut pas de passage par "radiornpb"
209 clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
210 ENDDO
211
212 if (nqmx >= 5) then
213 call press_coefoz ! read input pressure levels for ozone coefficients
214 end if
215 ENDIF
216
217 if (inirnpb) THEN
218 ! Initialisation du traceur dans le sol (couche limite radonique)
219 radio(1)= .true.
220 radio(2)= .true.
221 clsol(1)= .true.
222 clsol(2)= .true.
223 aerosol(2) = .TRUE. ! le Pb est un aerosol
224 call initrrnpb(pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, scavtr)
225 inirnpb=.false.
226 endif
227
228 if (convection) then
229 ! Calcul de l'effet de la convection
230 DO it=1, nqmx - 2
231 if (iflag_con == 2) then
232 ! Tiedke
233 CALL nflxtr(pdtphys, pmfu, pmfd, pde_u, pen_d, paprs, &
234 tr_seri(:, :, it), d_tr_cv(:, :, it))
235 else if (iflag_con == 3) then
236 ! Emanuel
237 call cvltr(pdtphys, da, phi, mp, paprs, tr_seri(1, 1, it), upwd, &
238 dnwd, d_tr_cv(1, 1, it))
239 endif
240
241 DO k = 1, llm
242 DO i = 1, klon
243 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cv(i, k, it)
244 ENDDO
245 ENDDO
246 WRITE(unit=itn, fmt='(i1)') it
247 CALL minmaxqfi(tr_seri(:, :, it), 0., 1e33, &
248 'convection, tracer index = ' // itn)
249 ENDDO
250 endif
251
252 ! Calcul de l'effet des thermiques
253
254 do it=1, nqmx - 2
255 do k=1, llm
256 do i=1, klon
257 d_tr_th(i, k, it)=0.
258 tr_seri(i, k, it)=max(tr_seri(i, k, it), 0.)
259 tr_seri(i, k, it)=min(tr_seri(i, k, it), 1e10)
260 enddo
261 enddo
262 enddo
263
264 if (iflag_thermals > 0) then
265 nsplit=10
266 DO it=1, nqmx - 2
267 do isplit=1, nsplit
268 ! Thermiques
269 call dqthermcell(klon, llm, pdtphys/nsplit &
270 , fm_therm, entr_therm, zmasse &
271 , tr_seri(1:klon, 1:llm, it), d_tr, ztra_th)
272
273 do k=1, llm
274 do i=1, klon
275 d_tr(i, k)=pdtphys*d_tr(i, k)/nsplit
276 d_tr_th(i, k, it)=d_tr_th(i, k, it)+d_tr(i, k)
277 tr_seri(i, k, it)=max(tr_seri(i, k, it)+d_tr(i, k), 0.)
278 enddo
279 enddo
280 enddo
281 ENDDO
282 endif
283
284 ! Calcul de l'effet de la couche limite
285
286 if (couchelimite) then
287 DO k = 1, llm
288 DO i = 1, klon
289 delp(i, k) = paprs(i, k)-paprs(i, k+1)
290 ENDDO
291 ENDDO
292
293 ! MAF modif pour tenir compte du cas traceur
294 DO it=1, nqmx - 2
295 if (clsol(it)) then
296 ! couche limite avec quantite dans le sol calculee
297 CALL cltracrn(it, pdtphys, yu1, yv1, &
298 coefh, t_seri, ftsol, pctsrf, &
299 tr_seri(:, :, it), trs(1, it), &
300 paprs, pplay, delp, &
301 masktr(1, it), fshtr(1, it), hsoltr(it), &
302 tautr(it), vdeptr(it), &
303 rlat, &
304 d_tr_cl(1, 1, it), d_trs)
305 DO k = 1, llm
306 DO i = 1, klon
307 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
308 ENDDO
309 ENDDO
310
311 ! Traceur ds sol
312
313 DO i = 1, klon
314 trs(i, it) = trs(i, it) + d_trs(i)
315 END DO
316 else ! couche limite avec flux prescrit
317 !MAF provisoire source / traceur a creer
318 DO i=1, klon
319 source(i) = 0.0 ! pas de source, pour l'instant
320 ENDDO
321
322 CALL cltrac(pdtphys, coefh, t_seri, &
323 tr_seri(1, 1, it), source, &
324 paprs, pplay, delp, &
325 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