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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 61 - (show annotations)
Fri Apr 20 14:58:43 2012 UTC (12 years ago) by guez
File size: 17381 byte(s)
No more included file in LMDZE, not even "netcdf.inc".

Created a variable containing the list of common source files in
GNUmakefile. So we now also see clearly files that are specific to
each program.

Split module "histcom". Assembled resulting files in directory
"Histcom".

Removed aliasing in calls to "laplacien".

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, 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 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 histsync_m, 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