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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 34 - (hide annotations)
Wed Jun 2 11:01:12 2010 UTC (13 years, 11 months ago) by guez
Original Path: trunk/libf/phylmd/phytrac.f90
File size: 17236 byte(s)
Split "ini_hist.f90" into single-procedure files.

In "calfis" and "physiq", removed dummy argument "nq" since "nq" must
be equal to "nqmx".

In "calfis", renamed dummy argument "pq" to "q", same name as actual
argument in "leapfrog". Renamed local variable "zqfi" to "qx", same
name as dummy argument in "physiq".

Removed arguments "itop_con" and "ibas_con" of "phytrac", which were
not used.

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

  ViewVC Help
Powered by ViewVC 1.1.21