/[lmdze]/trunk/libf/phylmd/phytrac.f90
ViewVC logotype

Contents of /trunk/libf/phylmd/phytrac.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 62 - (show annotations)
Thu Jul 26 14:37:37 2012 UTC (11 years, 9 months ago) by guez
File size: 17433 byte(s)
Changed handling of compiler in compilation system.

Removed the prefix letters "y", "p", "t" or "z" in some names of variables.

Replaced calls to NetCDF by calls to NetCDF95.

Extracted "ioget_calendar" procedures from "calendar.f90" into a
separate file.

Extracted to a separate file, "mathop2.f90", procedures that were not
part of the generic interface "mathop" in "mathop.f90".

Removed computation of "dq" in "bilan_dyn", which was not used.

In "iniadvtrac", removed schemes 20 Slopes and 30 Prather. Was not
compatible with declarations of array sizes.

In "clcdrag", "ustarhb", "vdif_kcay", "yamada4" and "coefkz", changed
the size of some arrays from "klon" to "knon".

Removed possible call to "conema3" in "physiq".

Removed unused argument "cd" in "yamada".

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

  ViewVC Help
Powered by ViewVC 1.1.21