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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 47 - (show annotations)
Fri Jul 1 15:00:48 2011 UTC (12 years, 10 months ago) by guez
File size: 17364 byte(s)
Split "thermcell.f" and "cv3_routines.f".
Removed copies of files that are now in "L_util".
Moved "mva9" and "diagetpq" to their own files.
Unified variable names across procedures.

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 ds phytrac
28 ! - stockage des moyennes des champs necessaires 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 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 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 ! KE
99 real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
100 REAL upwd(klon, llm) ! saturated updraft mass flux
101 REAL 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 ! Initialisation du traceur dans le sol (couche limite radonique)
240 if (inirnpb) THEN
241
242 radio(1)= .true.
243 radio(2)= .true.
244 clsol(1)= .true.
245 clsol(2)= .true.
246 aerosol(2) = .TRUE. ! le Pb est un aerosol
247
248 call initrrnpb(ftsol, pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, &
249 scavtr)
250 inirnpb=.false.
251 endif
252
253 ! Calcul de l'effet de la convection
254
255 if (convection) then
256 DO it=1, nq_phys
257 if (iflag_con.eq.2) then
258 ! tiedke
259 CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
260 paprs, tr_seri(1, 1, it), d_tr_cv(1, 1, it))
261 else if (iflag_con.eq.3) then
262 ! KE
263 call cvltr(pdtphys, da, phi, mp, paprs, &
264 tr_seri(1, 1, it), upwd, dnwd, d_tr_cv(1, 1, it))
265 endif
266
267 DO k = 1, llm
268 DO i = 1, klon
269 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cv(i, k, it)
270 ENDDO
271 ENDDO
272 WRITE(unit=itn, fmt='(i1)') it
273 CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, &
274 'convection, tracer index = ' // itn)
275 ENDDO
276 endif
277
278 ! Calcul de l'effet des thermiques
279
280 do it=1, nq_phys
281 do k=1, llm
282 do i=1, klon
283 d_tr_th(i, k, it)=0.
284 tr_seri(i, k, it)=max(tr_seri(i, k, it), 0.)
285 tr_seri(i, k, it)=min(tr_seri(i, k, it), 1.e10)
286 enddo
287 enddo
288 enddo
289
290 if (iflag_thermals > 0) then
291 nsplit=10
292 DO it=1, nq_phys
293 do isplit=1, nsplit
294 ! Thermiques
295 call dqthermcell(klon, llm, pdtphys/nsplit &
296 , fm_therm, entr_therm, zmasse &
297 , tr_seri(1:klon, 1:llm, it), d_tr, ztra_th)
298
299 do k=1, llm
300 do i=1, klon
301 d_tr(i, k)=pdtphys*d_tr(i, k)/nsplit
302 d_tr_th(i, k, it)=d_tr_th(i, k, it)+d_tr(i, k)
303 tr_seri(i, k, it)=max(tr_seri(i, k, it)+d_tr(i, k), 0.)
304 enddo
305 enddo
306 enddo
307 ENDDO
308 endif
309
310 ! Calcul de l'effet de la couche limite
311
312 if (couchelimite) then
313
314 DO k = 1, llm
315 DO i = 1, klon
316 delp(i, k) = paprs(i, k)-paprs(i, k+1)
317 ENDDO
318 ENDDO
319
320 ! MAF modif pour tenir compte du cas rnpb + traceur
321 DO it=1, nq_phys
322 if (clsol(it)) then
323 ! couche limite avec quantite dans le sol calculee
324 CALL cltracrn(it, pdtphys, yu1, yv1, &
325 coefh, t_seri, ftsol, pctsrf, &
326 tr_seri(1, 1, it), trs(1, it), &
327 paprs, pplay, delp, &
328 masktr(1, it), fshtr(1, it), hsoltr(it), &
329 tautr(it), vdeptr(it), &
330 rlat, &
331 d_tr_cl(1, 1, it), d_trs)
332 DO k = 1, llm
333 DO i = 1, klon
334 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
335 ENDDO
336 ENDDO
337
338 ! Traceur ds sol
339
340 DO i = 1, klon
341 trs(i, it) = trs(i, it) + d_trs(i)
342 END DO
343 else ! couche limite avec flux prescrit
344 !MAF provisoire source / traceur a creer
345 DO i=1, klon
346 source(i) = 0.0 ! pas de source, pour l'instant
347 ENDDO
348
349 CALL cltrac(pdtphys, coefh, t_seri, &
350 tr_seri(1, 1, it), source, &
351 paprs, pplay, delp, &
352 d_tr_cl(1, 1, it))
353 DO k = 1, llm
354 DO i = 1, klon
355 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
356 ENDDO
357 ENDDO
358 endif
359 ENDDO
360
361 endif ! couche limite
362
363 ! Calcul de l'effet du puits radioactif
364
365 ! MAF il faudrait faire une modification pour passer dans radiornpb
366 ! si radio=true mais pour l'instant radiornpb propre au cas rnpb
367 if (rnpb) then
368 d_tr_dec(:, :, :) = radiornpb(tr_seri, pdtphys, tautr)
369 DO it=1, nq_phys
370 if (radio(it)) then
371 tr_seri(:, :, it) = tr_seri(:, :, it) + d_tr_dec(:, :, it)
372 WRITE(unit=itn, fmt='(i1)') it
373 CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, 'puits rn it='//itn)
374 endif
375 ENDDO
376 endif ! rnpb decroissance radioactive
377
378 if (nq_phys >= 3) then
379 ! Ozone as a tracer:
380 if (mod(itap - 1, lmt_pas) == 0) then
381 ! Once per day, update the coefficients for ozone chemistry:
382 call regr_pr_comb_coefoz(julien)
383 end if
384 call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
385 end if
386
387 ! Calcul de l'effet de la precipitation
388
389 IF (lessivage) THEN
390 d_tr_lessi_nucl(:, :, :) = 0.
391 d_tr_lessi_impa(:, :, :) = 0.
392 flestottr(:, :, :) = 0.
393
394 ! tendance des aerosols nuclees et impactes
395
396 DO it = 1, nq_phys
397 IF (aerosol(it)) THEN
398 DO k = 1, llm
399 DO i = 1, klon
400 d_tr_lessi_nucl(i, k, it) = d_tr_lessi_nucl(i, k, it) + &
401 ( 1 - frac_nucl(i, k) )*tr_seri(i, k, it)
402 d_tr_lessi_impa(i, k, it) = d_tr_lessi_impa(i, k, it) + &
403 ( 1 - frac_impa(i, k) )*tr_seri(i, k, it)
404 ENDDO
405 ENDDO
406 ENDIF
407 ENDDO
408
409 ! Mises a jour des traceurs + calcul des flux de lessivage
410 ! Mise a jour due a l'impaction et a la nucleation
411
412 DO it = 1, nq_phys
413 IF (aerosol(it)) THEN
414 DO k = 1, llm
415 DO i = 1, klon
416 tr_seri(i, k, it)=tr_seri(i, k, it) &
417 *frac_impa(i, k)*frac_nucl(i, k)
418 ENDDO
419 ENDDO
420 ENDIF
421 ENDDO
422
423 ! Flux lessivage total
424
425 DO it = 1, nq_phys
426 DO k = 1, llm
427 DO i = 1, klon
428 flestottr(i, k, it) = flestottr(i, k, it) - &
429 ( d_tr_lessi_nucl(i, k, it) + &
430 d_tr_lessi_impa(i, k, it) ) * &
431 ( paprs(i, k)-paprs(i, k+1) ) / &
432 (RG * pdtphys)
433 ENDDO
434 ENDDO
435 ENDDO
436 ENDIF
437
438 ! Ecriture des sorties
439 call write_histrac(lessivage, nq_phys, itap, nid_tra)
440
441 if (lafin) then
442 print *, "C'est la fin de la physique."
443 open (unit=99, file='restarttrac', form='formatted')
444 do i=1, klon
445 write(unit=99, fmt=*) trs(i, 1)
446 enddo
447 PRINT *, 'Ecriture du fichier restarttrac'
448 close(99)
449 endif
450
451 contains
452
453 subroutine write_histrac(lessivage, nq_phys, itap, nid_tra)
454
455 ! From phylmd/write_histrac.h, version 1.9 2006/02/21 08:08:30
456
457 use dimens_m, only: iim, jjm, llm
458 use histcom, only: histsync
459 use histwrite_m, only: histwrite
460 use temps, only: itau_phy
461 use iniadvtrac_m, only: tnom
462 use comgeomphy, only: airephy
463 use dimphy, only: klon
464 use grid_change, only: gr_phy_write_2d
465 use gr_phy_write_3d_m, only: gr_phy_write_3d
466
467 logical, intent(in):: lessivage
468
469 integer, intent(in):: nq_phys
470 ! (nombre de traceurs auxquels on applique la physique)
471
472 integer, intent(in):: itap ! number of calls to "physiq"
473 integer, intent(in):: nid_tra
474
475 ! Variables local to the procedure:
476 integer it
477 integer itau_w ! pas de temps ecriture
478 logical, parameter:: ok_sync = .true.
479
480 !-----------------------------------------------------
481
482 itau_w = itau_phy + itap
483
484 CALL histwrite(nid_tra, "phis", itau_w, gr_phy_write_2d(pphis))
485 CALL histwrite(nid_tra, "aire", itau_w, gr_phy_write_2d(airephy))
486 CALL histwrite(nid_tra, "zmasse", itau_w, gr_phy_write_3d(zmasse))
487
488 DO it=1, nq_phys
489 CALL histwrite(nid_tra, tnom(it+2), itau_w, &
490 gr_phy_write_3d(tr_seri(:, :, it)))
491 if (lessivage) THEN
492 CALL histwrite(nid_tra, "fl"//tnom(it+2), itau_w, &
493 gr_phy_write_3d(flestottr(:, :, it)))
494 endif
495 CALL histwrite(nid_tra, "d_tr_th_"//tnom(it+2), itau_w, &
496 gr_phy_write_3d(d_tr_th(:, :, it)))
497 CALL histwrite(nid_tra, "d_tr_cv_"//tnom(it+2), itau_w, &
498 gr_phy_write_3d(d_tr_cv(:, :, it)))
499 CALL histwrite(nid_tra, "d_tr_cl_"//tnom(it+2), itau_w, &
500 gr_phy_write_3d(d_tr_cl(:, :, it)))
501 ENDDO
502
503 CALL histwrite(nid_tra, "pplay", itau_w, gr_phy_write_3d(pplay))
504 CALL histwrite(nid_tra, "T", itau_w, gr_phy_write_3d(t_seri))
505
506 if (ok_sync) then
507 call histsync(nid_tra)
508 endif
509
510 end subroutine write_histrac
511
512 END SUBROUTINE phytrac
513
514 end module phytrac_m

  ViewVC Help
Powered by ViewVC 1.1.21