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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 36 - (hide annotations)
Thu Dec 2 17:11:04 2010 UTC (13 years, 5 months ago) by guez
File size: 17209 byte(s)
Now using the library "NR_util".

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

  ViewVC Help
Powered by ViewVC 1.1.21