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

Annotation of /trunk/Sources/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 17 - (hide annotations)
Tue Aug 5 13:31:32 2008 UTC (15 years, 9 months ago) by guez
Original Path: trunk/libf/phylmd/phytrac.f90
File size: 17956 byte(s)
Created rule for "compare_sampl_*" files in
"Documentation/Manuel_LMDZE.texfol/Graphiques/GNUmakefile".

Extracted "qcheck", "radiornpb", "minmaxqfi" into separate files.

Read pressure coordinate of ozone coefficients once per run instead of
every day.

Added some "intent" attributes.

Added argument "nq" to "ini_histday". Replaced calls to "gr_fi_ecrit"
by calls to "gr_phy_write_2d". "Sigma_O3_Royer" is written to
"histday.nc" only if "nq >= 4". Moved "ini_histrac" to module
"ini_hist".

Compute "zmasse" in "physiq", pass it to "phytrac".

Removed computations of "pftsol*" and "ppsrf*" in "phytrac".

Do not use variable "rg" from module "YOMCST" in "TLIFT".

1 guez 3 module phytrac_m
2    
3     ! This module is clean: no C preprocessor directive, no include line.
4    
5     IMPLICIT none
6    
7     private
8     public phytrac
9    
10     contains
11    
12 guez 7 SUBROUTINE phytrac(rnpb, itap, lmt_pas, julien, gmtime, firstcal, lafin, &
13     nqmax, pdtphys, u, v, t_seri, paprs, pplay, pmfu, pmfd, pen_u, &
14 guez 3 pde_u, pen_d, pde_d, coefh, fm_therm, entr_therm, yu1, yv1, ftsol, &
15     pctsrf, frac_impa, frac_nucl, presnivs, pphis, &
16 guez 12 pphi, albsol, rh, cldfra, rneb, diafra, cldliq, itop_con, &
17 guez 17 ibas_con, pmflxr, pmflxs, prfl, psfl, da, phi, mp, upwd, dnwd, &
18     tr_seri, zmasse)
19 guez 3
20     ! From phylmd/phytrac.F, version 1.15 2006/02/21 08:08:30
21    
22     ! Authors : Frédéric Hourdin, Abderrahmane Idelkadi, Marie-Alice
23     ! Foujols, Olivia
24     ! Objet : moniteur général des tendances des traceurs
25    
26     ! Remarques :
27     ! 1/ L'appel de "phytrac" se fait avec "nq-2" donc nous avons bien
28     ! les vrais traceurs (en nombre "nbtr", sans la vapeur d'eau ni l'eau
29     ! liquide) dans "phytrac".
30     ! 2/ Le choix du radon et du plomb se fait juste avec un "data"
31     ! (peu propre).
32     ! Pourrait-on avoir une variable qui indiquerait le type de traceur ?
33    
34 guez 17 use dimens_m, only: llm
35 guez 3 use indicesol, only: nbsrf
36     use dimphy, only: klon, nbtr
37 guez 12 use clesphys, only: ecrit_tra
38     use clesphys2, only: iflag_con
39 guez 3 use abort_gcm_m, only: abort_gcm
40     use YOMCST, only: rg
41     use ctherm, only: iflag_thermals
42 guez 7 use regr_pr_comb_coefoz_m, only: regr_pr_comb_coefoz
43 guez 3 use phyetat0_m, only: rlat
44     use o3_chem_m, only: o3_chem
45 guez 17 use ini_hist, only: ini_histrac
46     use radiornpb_m, only: radiornpb
47     use minmaxqfi_m, only: minmaxqfi
48     use numer_rec, only: assert
49     use press_coefoz_m, only: press_coefoz
50 guez 3
51     ! Arguments:
52    
53     ! EN ENTREE:
54    
55     ! divers:
56    
57     logical, intent(in):: rnpb
58    
59     integer, intent(in):: nqmax
60     ! (nombre de traceurs auxquels on applique la physique)
61    
62 guez 7 integer, intent(in):: itap ! number of calls to "physiq"
63     integer, intent(in):: lmt_pas ! number of time steps of "physics" per day
64 guez 3 integer, intent(in):: julien !jour julien, 1 <= julien <= 360
65     integer itop_con(klon)
66     integer ibas_con(klon)
67     real, intent(in):: gmtime ! heure de la journée en fraction de jour
68 guez 7 real, intent(in):: pdtphys ! pas d'integration pour la physique (s)
69 guez 3 real, intent(in):: t_seri(klon, llm) ! temperature, in K
70    
71 guez 17 real, intent(inout):: tr_seri(klon, llm, nbtr)
72 guez 3 ! (mass fractions of tracers, excluding water, at mid-layers)
73    
74     real u(klon, llm)
75     real v(klon, llm)
76     real rh(klon, llm) ! humidite relative
77     real cldliq(klon, llm) ! eau liquide nuageuse
78     real cldfra(klon, llm) ! fraction nuageuse (tous les nuages)
79    
80     real diafra(klon, llm)
81     ! (fraction nuageuse (convection ou stratus artificiels))
82    
83     real rneb(klon, llm) ! fraction nuageuse (grande echelle)
84     real albsol(klon) ! albedo surface
85    
86     real, intent(in):: paprs(klon, llm+1)
87     ! (pression pour chaque inter-couche, en Pa)
88    
89 guez 10 real, intent(in):: pplay(klon, llm)
90     ! (pression pour le mileu de chaque couche, en Pa)
91    
92 guez 3 real pphi(klon, llm) ! geopotentiel
93     real pphis(klon)
94     REAL, intent(in):: presnivs(llm)
95 guez 7 logical, intent(in):: firstcal ! first call to "calfis"
96 guez 3 logical, intent(in):: lafin ! fin de la physique
97    
98     REAL pmflxr(klon, llm+1), pmflxs(klon, llm+1) !--lessivage convection
99     REAL prfl(klon, llm+1), psfl(klon, llm+1) !--lessivage large-scale
100    
101     ! convection:
102     REAL pmfu(klon, llm) ! flux de masse dans le panache montant
103     REAL pmfd(klon, llm) ! flux de masse dans le panache descendant
104     REAL pen_u(klon, llm) ! flux entraine dans le panache montant
105    
106     ! thermiques:
107    
108     real fm_therm(klon, llm+1), entr_therm(klon, llm)
109    
110     REAL pde_u(klon, llm) ! flux detraine dans le panache montant
111     REAL pen_d(klon, llm) ! flux entraine dans le panache descendant
112     REAL pde_d(klon, llm) ! flux detraine dans le panache descendant
113     ! KE
114     real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
115     REAL upwd(klon, llm) ! saturated updraft mass flux
116     REAL dnwd(klon, llm) ! saturated downdraft mass flux
117    
118     ! Couche limite:
119    
120     REAL coefh(klon, llm) ! coeff melange CL
121     REAL yu1(klon) ! vents au premier niveau
122     REAL yv1(klon) ! vents au premier niveau
123    
124     ! Lessivage:
125    
126     ! pour le ON-LINE
127    
128     REAL frac_impa(klon, llm) ! fraction d'aerosols impactes
129     REAL frac_nucl(klon, llm) ! fraction d'aerosols nuclees
130    
131     ! Arguments necessaires pour les sources et puits de traceur:
132    
133     real ftsol(klon, nbsrf) ! Temperature du sol (surf)(Kelvin)
134     real pctsrf(klon, nbsrf) ! Pourcentage de sol f(nature du sol)
135    
136 guez 17 real, intent(in):: zmasse(:, :) ! (klon, llm)
137     ! (column-density of mass of air in a cell, in kg m-2)
138 guez 3
139 guez 17 ! Variables local to the procedure:
140 guez 3
141 guez 17 integer nsplit
142    
143     ! TRACEURS
144    
145 guez 3 ! Sources et puits des traceurs:
146    
147     ! Pour l'instant seuls les cas du rn et du pb ont ete envisages.
148    
149     REAL source(klon) ! a voir lorsque le flux est prescrit
150     !
151     ! Pour la source de radon et son reservoir de sol
152    
153     REAL, save:: trs(klon, nbtr) ! Concentration de radon dans le sol
154    
155     REAL masktr(klon, nbtr) ! Masque reservoir de sol traceur
156     ! Masque de l'echange avec la surface
157     ! (1 = reservoir) ou (possible => 1 )
158     SAVE masktr
159     REAL fshtr(klon, nbtr) ! Flux surfacique dans le reservoir de sol
160     SAVE fshtr
161     REAL hsoltr(nbtr) ! Epaisseur equivalente du reservoir de sol
162     SAVE hsoltr
163     REAL tautr(nbtr) ! Constante de decroissance radioactive
164     SAVE tautr
165     REAL vdeptr(nbtr) ! Vitesse de depot sec dans la couche Brownienne
166     SAVE vdeptr
167     REAL scavtr(nbtr) ! Coefficient de lessivage
168     SAVE scavtr
169    
170     CHARACTER itn
171     INTEGER, save:: nid_tra
172    
173     ! nature du traceur
174    
175     logical aerosol(nbtr) ! Nature du traceur
176     ! ! aerosol(it) = true => aerosol
177     ! ! aerosol(it) = false => gaz
178     logical clsol(nbtr) ! couche limite sol calculée
179     logical radio(nbtr) ! décroisssance radioactive
180     save aerosol, clsol, radio
181    
182     ! convection tiedtke
183     INTEGER i, k, it
184     REAL delp(klon, llm)
185    
186     ! Variables liees a l'ecriture de la bande histoire physique
187    
188     ! Variables locales pour effectuer les appels en serie
189    
190     REAL d_tr(klon, llm), d_trs(klon) ! tendances de traceurs
191     REAL d_tr_cl(klon, llm, nbtr) ! tendance de traceurs couche limite
192     REAL d_tr_cv(klon, llm, nbtr) ! tendance de traceurs conv pour chq traceur
193     REAL d_tr_th(klon, llm, nbtr) ! la tendance des thermiques
194     REAL d_tr_dec(klon, llm, 2) ! la tendance de la decroissance
195     ! ! radioactive du rn - > pb
196     REAL d_tr_lessi_impa(klon, llm, nbtr) ! la tendance du lessivage
197     ! ! par impaction
198     REAL d_tr_lessi_nucl(klon, llm, nbtr) ! la tendance du lessivage
199     ! ! par nucleation
200     REAL flestottr(klon, llm, nbtr) ! flux de lessivage
201     ! ! dans chaque couche
202    
203     real ztra_th(klon, llm)
204     integer isplit
205    
206     ! Controls:
207     logical:: couchelimite = .true.
208     logical:: convection = .true.
209     logical:: lessivage = .true.
210     logical, save:: inirnpb
211    
212     !--------------------------------------
213    
214 guez 17 call assert(shape(zmasse) == (/klon, llm/), "phytrac")
215 guez 3
216 guez 7 if (firstcal) then
217 guez 3 print *, 'phytrac: pdtphys = ', pdtphys
218     PRINT *, 'Fréquence de sortie des traceurs : ecrit_tra = ', ecrit_tra
219 guez 17 if (nbtr < nqmax) call abort_gcm('phytrac', 'See above', 1)
220 guez 3 inirnpb=rnpb
221    
222     ! Initialisation des sorties :
223     call ini_histrac(nid_tra, pdtphys, presnivs, nqmax, lessivage)
224    
225     ! Initialisation de certaines variables pour le radon et le plomb
226     ! Initialisation du traceur dans le sol (couche limite radonique)
227     trs(:, :) = 0.
228    
229     open (unit=99, file='starttrac', status='old', err=999, &
230     form='formatted')
231     read(unit=99, fmt=*) (trs(i, 1), i=1, klon)
232     999 continue
233     close(unit=99)
234    
235     ! Initialisation de la fraction d'aerosols lessivee
236    
237     d_tr_lessi_impa(:, :, :) = 0.
238     d_tr_lessi_nucl(:, :, :) = 0.
239    
240     ! Initialisation de la nature des traceurs
241    
242     DO it = 1, nqmax
243     aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
244     radio(it) = .FALSE. ! par défaut pas de passage par "radiornpb"
245     clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
246     ENDDO
247 guez 17
248     if (nqmax >= 3) then
249     call press_coefoz ! read input pressure levels for ozone coefficients
250     end if
251 guez 3 ENDIF
252    
253     ! Initialisation du traceur dans le sol (couche limite radonique)
254     if (inirnpb) THEN
255    
256     radio(1)= .true.
257     radio(2)= .true.
258     clsol(1)= .true.
259     clsol(2)= .true.
260     aerosol(2) = .TRUE. ! le Pb est un aerosol
261    
262     call initrrnpb(ftsol, pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, &
263     scavtr)
264     inirnpb=.false.
265     endif
266    
267     ! Calcul de l'effet de la convection
268    
269     if (convection) then
270     DO it=1, nqmax
271     if (iflag_con.eq.2) then
272     ! tiedke
273     CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
274 guez 10 paprs, tr_seri(1, 1, it), d_tr_cv(1, 1, it))
275 guez 3 else if (iflag_con.eq.3) then
276     ! KE
277 guez 10 call cvltr(pdtphys, da, phi, mp, paprs, &
278 guez 3 tr_seri(1, 1, it), upwd, dnwd, d_tr_cv(1, 1, it))
279     endif
280    
281     DO k = 1, llm
282     DO i = 1, klon
283     tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cv(i, k, it)
284     ENDDO
285     ENDDO
286     WRITE(unit=itn, fmt='(i1)') it
287     CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, &
288     'convection, tracer index = ' // itn)
289     ENDDO
290     endif
291    
292     ! Calcul de l'effet des thermiques
293    
294     do it=1, nqmax
295     do k=1, llm
296     do i=1, klon
297     d_tr_th(i, k, it)=0.
298     tr_seri(i, k, it)=max(tr_seri(i, k, it), 0.)
299     tr_seri(i, k, it)=min(tr_seri(i, k, it), 1.e10)
300     enddo
301     enddo
302     enddo
303    
304     if (iflag_thermals > 0) then
305     nsplit=10
306     DO it=1, nqmax
307     do isplit=1, nsplit
308     ! Thermiques
309     call dqthermcell(klon, llm, pdtphys/nsplit &
310     , fm_therm, entr_therm, zmasse &
311     , tr_seri(1:klon, 1:llm, it), d_tr, ztra_th)
312    
313     do k=1, llm
314     do i=1, klon
315     d_tr(i, k)=pdtphys*d_tr(i, k)/nsplit
316     d_tr_th(i, k, it)=d_tr_th(i, k, it)+d_tr(i, k)
317     tr_seri(i, k, it)=max(tr_seri(i, k, it)+d_tr(i, k), 0.)
318     enddo
319     enddo
320     enddo
321     ENDDO
322     endif
323    
324     ! Calcul de l'effet de la couche limite
325    
326     if (couchelimite) then
327    
328     DO k = 1, llm
329     DO i = 1, klon
330     delp(i, k) = paprs(i, k)-paprs(i, k+1)
331     ENDDO
332     ENDDO
333    
334     ! MAF modif pour tenir compte du cas rnpb + traceur
335     DO it=1, nqmax
336     if (clsol(it)) then
337     ! couche limite avec quantite dans le sol calculee
338     CALL cltracrn(it, pdtphys, yu1, yv1, &
339     coefh, t_seri, ftsol, pctsrf, &
340     tr_seri(1, 1, it), trs(1, it), &
341     paprs, pplay, delp, &
342     masktr(1, it), fshtr(1, it), hsoltr(it), &
343     tautr(it), vdeptr(it), &
344     rlat, &
345     d_tr_cl(1, 1, it), d_trs)
346     DO k = 1, llm
347     DO i = 1, klon
348     tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
349     ENDDO
350     ENDDO
351    
352     ! Traceur ds sol
353    
354     DO i = 1, klon
355     trs(i, it) = trs(i, it) + d_trs(i)
356     END DO
357     else ! couche limite avec flux prescrit
358     !MAF provisoire source / traceur a creer
359     DO i=1, klon
360     source(i) = 0.0 ! pas de source, pour l'instant
361     ENDDO
362    
363     CALL cltrac(pdtphys, coefh, t_seri, &
364     tr_seri(1, 1, it), source, &
365     paprs, pplay, delp, &
366     d_tr_cl(1, 1, it))
367     DO k = 1, llm
368     DO i = 1, klon
369     tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
370     ENDDO
371     ENDDO
372     endif
373     ENDDO
374    
375     endif ! couche limite
376    
377     ! Calcul de l'effet du puits radioactif
378    
379     ! MAF il faudrait faire une modification pour passer dans radiornpb
380     ! si radio=true mais pour l'instant radiornpb propre au cas rnpb
381     if (rnpb) then
382     d_tr_dec(:, :, :) = radiornpb(tr_seri, pdtphys, tautr)
383     DO it=1, nqmax
384     if (radio(it)) then
385     tr_seri(:, :, it) = tr_seri(:, :, it) + d_tr_dec(:, :, it)
386     WRITE(unit=itn, fmt='(i1)') it
387     CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, 'puits rn it='//itn)
388     endif
389     ENDDO
390     endif ! rnpb decroissance radioactive
391    
392 guez 6 if (nqmax >= 3) then
393     ! Ozone as a tracer:
394 guez 7 if (mod(itap - 1, lmt_pas) == 0) then
395     ! Once per day, update the coefficients for ozone chemistry:
396     call regr_pr_comb_coefoz(julien)
397     end if
398 guez 6 call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
399     end if
400 guez 3
401     ! Calcul de l'effet de la precipitation
402    
403     IF (lessivage) THEN
404     d_tr_lessi_nucl(:, :, :) = 0.
405     d_tr_lessi_impa(:, :, :) = 0.
406     flestottr(:, :, :) = 0.
407    
408     ! tendance des aerosols nuclees et impactes
409    
410     DO it = 1, nqmax
411     IF (aerosol(it)) THEN
412     DO k = 1, llm
413     DO i = 1, klon
414     d_tr_lessi_nucl(i, k, it) = d_tr_lessi_nucl(i, k, it) + &
415     ( 1 - frac_nucl(i, k) )*tr_seri(i, k, it)
416     d_tr_lessi_impa(i, k, it) = d_tr_lessi_impa(i, k, it) + &
417     ( 1 - frac_impa(i, k) )*tr_seri(i, k, it)
418     ENDDO
419     ENDDO
420     ENDIF
421     ENDDO
422    
423     ! Mises a jour des traceurs + calcul des flux de lessivage
424     ! Mise a jour due a l'impaction et a la nucleation
425    
426     DO it = 1, nqmax
427     IF (aerosol(it)) THEN
428     DO k = 1, llm
429     DO i = 1, klon
430     tr_seri(i, k, it)=tr_seri(i, k, it) &
431     *frac_impa(i, k)*frac_nucl(i, k)
432     ENDDO
433     ENDDO
434     ENDIF
435     ENDDO
436    
437     ! Flux lessivage total
438    
439     DO it = 1, nqmax
440     DO k = 1, llm
441     DO i = 1, klon
442     flestottr(i, k, it) = flestottr(i, k, it) - &
443     ( d_tr_lessi_nucl(i, k, it) + &
444     d_tr_lessi_impa(i, k, it) ) * &
445     ( paprs(i, k)-paprs(i, k+1) ) / &
446     (RG * pdtphys)
447     ENDDO
448     ENDDO
449     ENDDO
450     ENDIF
451    
452     ! Ecriture des sorties
453 guez 7 call write_histrac(lessivage, nqmax, itap, nid_tra)
454 guez 3
455     if (lafin) then
456     print *, "C'est la fin de la physique."
457     open (unit=99, file='restarttrac', form='formatted')
458     do i=1, klon
459     write(unit=99, fmt=*) trs(i, 1)
460     enddo
461     PRINT *, 'Ecriture du fichier restarttrac'
462     close(99)
463     endif
464    
465     contains
466    
467 guez 7 subroutine write_histrac(lessivage, nqmax, itap, nid_tra)
468 guez 3
469     ! From phylmd/write_histrac.h, version 1.9 2006/02/21 08:08:30
470    
471     use dimens_m, only: iim, jjm, llm
472     use ioipsl, only: histwrite, histsync
473     use temps, only: itau_phy
474     use advtrac_m, only: tnom
475     use comgeomphy, only: airephy
476     use dimphy, only: klon
477    
478     logical, intent(in):: lessivage
479    
480     integer, intent(in):: nqmax
481     ! (nombre de traceurs auxquels on applique la physique)
482    
483 guez 7 integer, intent(in):: itap ! number of calls to "physiq"
484 guez 3 integer, intent(in):: nid_tra
485    
486     ! Variables local to the procedure:
487     integer it
488     integer itau_w ! pas de temps ecriture
489     REAL zx_tmp_2d(iim, jjm+1), zx_tmp_3d(iim, jjm+1, llm)
490     logical, parameter:: ok_sync = .true.
491    
492     !-----------------------------------------------------
493    
494 guez 7 itau_w = itau_phy + itap
495 guez 3
496     CALL gr_fi_ecrit(1, klon, iim, jjm+1, pphis, zx_tmp_2d)
497 guez 15 CALL histwrite(nid_tra, "phis", itau_w, zx_tmp_2d)
498 guez 3
499     CALL gr_fi_ecrit(1, klon, iim, jjm+1, airephy, zx_tmp_2d)
500 guez 15 CALL histwrite(nid_tra, "aire", itau_w, zx_tmp_2d)
501 guez 3
502     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, zmasse, zx_tmp_3d)
503 guez 15 CALL histwrite(nid_tra, "zmasse", itau_w, zx_tmp_3d)
504 guez 3
505     DO it=1, nqmax
506     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, tr_seri(1, 1, it), zx_tmp_3d)
507 guez 15 CALL histwrite(nid_tra, tnom(it+2), itau_w, zx_tmp_3d)
508 guez 3 if (lessivage) THEN
509     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, flestottr(1, 1, it), &
510     zx_tmp_3d)
511 guez 15 CALL histwrite(nid_tra, "fl"//tnom(it+2), itau_w, zx_tmp_3d)
512 guez 3 endif
513    
514     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_th(1, 1, it), zx_tmp_3d)
515 guez 15 CALL histwrite(nid_tra, "d_tr_th_"//tnom(it+2), itau_w, zx_tmp_3d)
516 guez 3 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_cv(1, 1, it), zx_tmp_3d)
517 guez 15 CALL histwrite(nid_tra, "d_tr_cv_"//tnom(it+2), itau_w, zx_tmp_3d)
518 guez 3 CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_cl(1, 1, it), zx_tmp_3d)
519 guez 15 CALL histwrite(nid_tra, "d_tr_cl_"//tnom(it+2), itau_w, zx_tmp_3d)
520 guez 3 ENDDO
521    
522     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pplay, zx_tmp_3d)
523 guez 15 CALL histwrite(nid_tra, "pplay", itau_w, zx_tmp_3d)
524 guez 3
525     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, t_seri, zx_tmp_3d)
526 guez 15 CALL histwrite(nid_tra, "t", itau_w, zx_tmp_3d)
527 guez 3
528     if (ok_sync) then
529     call histsync(nid_tra)
530     endif
531    
532     end subroutine write_histrac
533    
534     END SUBROUTINE phytrac
535    
536     end module phytrac_m

  ViewVC Help
Powered by ViewVC 1.1.21