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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 47 - (hide annotations)
Fri Jul 1 15:00:48 2011 UTC (12 years, 10 months ago) by guez
Original Path: trunk/libf/phylmd/phytrac.f90
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 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 47 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 guez 3
16 guez 47 ! From phylmd/phytrac.F, version 1.15 2006/02/21 08:08:30 (SVN revision 679)
17 guez 3
18 guez 18 ! Authors: Frédéric Hourdin, Abderrahmane Idelkadi, Marie-Alice
19 guez 3 ! Foujols, Olivia
20     ! Objet : moniteur général des tendances des traceurs
21    
22 guez 34 ! L'appel de "phytrac" se fait avec "nqmx-2" donc nous avons bien
23 guez 3 ! les vrais traceurs (en nombre "nbtr", sans la vapeur d'eau ni l'eau
24     ! liquide) dans "phytrac".
25    
26 guez 47 ! Modifications pour les traceurs :
27     ! - uniformisation des parametrisations ds phytrac
28     ! - stockage des moyennes des champs necessaires en mode traceur off-line
29    
30 guez 17 use dimens_m, only: llm
31 guez 3 use indicesol, only: nbsrf
32     use dimphy, only: klon, nbtr
33 guez 12 use clesphys, only: ecrit_tra
34     use clesphys2, only: iflag_con
35 guez 3 use abort_gcm_m, only: abort_gcm
36 guez 38 use SUPHEC_M, only: rg
37 guez 3 use ctherm, only: iflag_thermals
38 guez 7 use regr_pr_comb_coefoz_m, only: regr_pr_comb_coefoz
39 guez 3 use phyetat0_m, only: rlat
40     use o3_chem_m, only: o3_chem
41 guez 34 use ini_histrac_m, only: ini_histrac
42 guez 17 use radiornpb_m, only: radiornpb
43     use minmaxqfi_m, only: minmaxqfi
44 guez 36 use nr_util, only: assert
45 guez 17 use press_coefoz_m, only: press_coefoz
46 guez 3
47     logical, intent(in):: rnpb
48    
49 guez 18 integer, intent(in):: nq_phys
50 guez 3 ! (nombre de traceurs auxquels on applique la physique)
51    
52 guez 7 integer, intent(in):: itap ! number of calls to "physiq"
53     integer, intent(in):: lmt_pas ! number of time steps of "physics" per day
54 guez 3 integer, intent(in):: julien !jour julien, 1 <= julien <= 360
55     real, intent(in):: gmtime ! heure de la journée en fraction de jour
56 guez 7 real, intent(in):: pdtphys ! pas d'integration pour la physique (s)
57 guez 3 real, intent(in):: t_seri(klon, llm) ! temperature, in K
58    
59 guez 18 real, intent(inout):: tr_seri(:, :, :) ! (klon, llm, nbtr)
60 guez 3 ! (mass fractions of tracers, excluding water, at mid-layers)
61    
62 guez 47 real, intent(in):: u(klon, llm)
63 guez 3 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 guez 10 real, intent(in):: pplay(klon, llm)
77     ! (pression pour le mileu de chaque couche, en Pa)
78    
79 guez 3 real pphis(klon)
80 guez 7 logical, intent(in):: firstcal ! first call to "calfis"
81 guez 3 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 guez 17 real, intent(in):: zmasse(:, :) ! (klon, llm)
122     ! (column-density of mass of air in a cell, in kg m-2)
123 guez 3
124 guez 17 ! Variables local to the procedure:
125 guez 3
126 guez 17 integer nsplit
127    
128     ! TRACEURS
129    
130 guez 3 ! 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 guez 18 call assert(shape(zmasse) == (/klon, llm/), "phytrac zmasse")
200     call assert(shape(tr_seri) == (/klon, llm, nbtr/), "phytrac tr_seri")
201 guez 3
202 guez 7 if (firstcal) then
203 guez 3 print *, 'phytrac: pdtphys = ', pdtphys
204     PRINT *, 'Fréquence de sortie des traceurs : ecrit_tra = ', ecrit_tra
205 guez 18 if (nbtr < nq_phys) call abort_gcm('phytrac', 'nbtr < nq_phys', 1)
206 guez 3 inirnpb=rnpb
207    
208     ! Initialisation des sorties :
209 guez 20 call ini_histrac(nid_tra, pdtphys, nq_phys, lessivage)
210 guez 3
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 guez 18 DO it = 1, nq_phys
229 guez 3 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 guez 17
234 guez 18 if (nq_phys >= 3) then
235 guez 17 call press_coefoz ! read input pressure levels for ozone coefficients
236     end if
237 guez 3 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 guez 18 DO it=1, nq_phys
257 guez 3 if (iflag_con.eq.2) then
258     ! tiedke
259     CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
260 guez 10 paprs, tr_seri(1, 1, it), d_tr_cv(1, 1, it))
261 guez 3 else if (iflag_con.eq.3) then
262     ! KE
263 guez 10 call cvltr(pdtphys, da, phi, mp, paprs, &
264 guez 3 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 guez 18 do it=1, nq_phys
281 guez 3 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 guez 18 DO it=1, nq_phys
293 guez 3 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 guez 18 DO it=1, nq_phys
322 guez 3 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 guez 18 DO it=1, nq_phys
370 guez 3 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 guez 18 if (nq_phys >= 3) then
379 guez 6 ! Ozone as a tracer:
380 guez 7 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 guez 6 call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
385     end if
386 guez 3
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 guez 18 DO it = 1, nq_phys
397 guez 3 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 guez 18 DO it = 1, nq_phys
413 guez 3 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 guez 18 DO it = 1, nq_phys
426 guez 3 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 guez 18 call write_histrac(lessivage, nq_phys, itap, nid_tra)
440 guez 3
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 guez 18 subroutine write_histrac(lessivage, nq_phys, itap, nid_tra)
454 guez 3
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 guez 30 use histcom, only: histsync
459     use histwrite_m, only: histwrite
460 guez 3 use temps, only: itau_phy
461 guez 18 use iniadvtrac_m, only: tnom
462 guez 3 use comgeomphy, only: airephy
463     use dimphy, only: klon
464 guez 32 use grid_change, only: gr_phy_write_2d
465     use gr_phy_write_3d_m, only: gr_phy_write_3d
466 guez 3
467     logical, intent(in):: lessivage
468    
469 guez 18 integer, intent(in):: nq_phys
470 guez 3 ! (nombre de traceurs auxquels on applique la physique)
471    
472 guez 7 integer, intent(in):: itap ! number of calls to "physiq"
473 guez 3 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 guez 7 itau_w = itau_phy + itap
483 guez 3
484 guez 20 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 guez 3
488 guez 18 DO it=1, nq_phys
489 guez 20 CALL histwrite(nid_tra, tnom(it+2), itau_w, &
490     gr_phy_write_3d(tr_seri(:, :, it)))
491 guez 3 if (lessivage) THEN
492 guez 20 CALL histwrite(nid_tra, "fl"//tnom(it+2), itau_w, &
493     gr_phy_write_3d(flestottr(:, :, it)))
494 guez 3 endif
495 guez 20 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 guez 3 ENDDO
502    
503 guez 20 CALL histwrite(nid_tra, "pplay", itau_w, gr_phy_write_3d(pplay))
504 guez 22 CALL histwrite(nid_tra, "T", itau_w, gr_phy_write_3d(t_seri))
505 guez 3
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