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

Annotation of /trunk/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 23 - (hide annotations)
Mon Dec 14 15:25:16 2009 UTC (14 years, 5 months ago) by guez
Original Path: trunk/libf/phylmd/phytrac.f90
File size: 17553 byte(s)
Split "orografi.f": one file for each procedure. Put the created files
in new directory "Orography".

Removed argument "vcov" of procedure "sortvarc". Removed arguments
"itau" and "time" of procedure "caldyn0". Removed arguments "itau",
"time" and "vcov" of procedure "sortvarc0".

Removed argument "time" of procedure "dynredem1". Removed NetCDF
variable "temps" in files "start.nc" and "restart.nc", because its
value is always 0.

Removed argument "nq" of procedures "iniadvtrac" and "leapfrog". The
number of "tracers read in "traceur.def" must now be equal to "nqmx",
or "nqmx" must equal 4 if there is no file "traceur.def". Replaced
variable "nq" by constant "nqmx" in "leapfrog".

NetCDF variable for ozone field in "coefoz.nc" must now be called
"tro3" instead of "r".

Fixed bug in "zenang".

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

  ViewVC Help
Powered by ViewVC 1.1.21