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

Contents of /trunk/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 104 - (show annotations)
Thu Sep 4 10:05:52 2014 UTC (9 years, 8 months ago) by guez
File size: 16502 byte(s)
Removed procedure sortvarc0. Called sortvarc with an additional
argument resetvarc instead. (Following LMDZ.) Moved current time
computations and some printing statements from sortvarc to
caldyn. Could then remove arguments itau and time_0 of sortvarc, and
could remove "use dynetat0". Better to keep "dynetat0.f" as a gcm-only
file.

Moved some variables from module ener to module sortvarc.

Split file "mathelp.f" into single-procedure files.

Removed unused argument nadv of adaptdt. Removed dimension arguments
of bernoui.

Removed unused argument nisurf of interfoce_lim. Changed the size of
argument lmt_sst of interfoce_lim from klon to knon. Removed case when
newlmt is false.

dynredem1 is called only once in each run, either ce0l or gcm. So
variable nb in call to nf95_put_var was always 1. Removed variable nb.

Removed dimension arguments of calcul_fluxs. Removed unused arguments
precip_rain, precip_snow, snow of calcul_fluxs. Changed the size of
all the arrays in calcul_fluxs from klon to knon.

Removed dimension arguments of fonte_neige. Changed the size of all
the arrays in fonte_neige from klon to knon.

Changed the size of arguments tsurf and tsurf_new of interfsurf_hq
from klon to knon. Changed the size of argument ptsrf of soil from
klon to knon.

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

  ViewVC Help
Powered by ViewVC 1.1.21