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

Annotation of /trunk/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 20 - (hide annotations)
Wed Oct 15 16:19:57 2008 UTC (15 years, 7 months ago) by guez
Original Path: trunk/libf/phylmd/phytrac.f90
File size: 17614 byte(s)
Deleted argument "presnivs" of "physiq", "ini_histhf", "ini_histhf3d",
"ini_histday", "ini_histins", "ini_histrac", "phytrac". Access it from
"comvert" instead.

Replaced calls to NetCDF Fortran 77 interface by calls to Fortran 90
interface or to NetCDF95.

Procedure "gr_phy_write_3d" now works with a variable of arbitrary
size in the second dimension.

Annotated use statements with "only" clause.

Replaced calls to NetCDF interface version 2 by calls to Fortran 90
interface in "guide.f90" and "read_reanalyse.f".

In "write_histrac", replaced calls to "gr_fi_ecrit" by calls to
"gr_phy_write_2d" and "gr_phy_write_3d".

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     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 guez 20 CALL histwrite(nid_tra, "phis", itau_w, gr_phy_write_2d(pphis))
497     CALL histwrite(nid_tra, "aire", itau_w, gr_phy_write_2d(airephy))
498     CALL histwrite(nid_tra, "zmasse", itau_w, gr_phy_write_3d(zmasse))
499 guez 3
500 guez 18 DO it=1, nq_phys
501 guez 20 CALL histwrite(nid_tra, tnom(it+2), itau_w, &
502     gr_phy_write_3d(tr_seri(:, :, it)))
503 guez 3 if (lessivage) THEN
504 guez 20 CALL histwrite(nid_tra, "fl"//tnom(it+2), itau_w, &
505     gr_phy_write_3d(flestottr(:, :, it)))
506 guez 3 endif
507 guez 20 CALL histwrite(nid_tra, "d_tr_th_"//tnom(it+2), itau_w, &
508     gr_phy_write_3d(d_tr_th(:, :, it)))
509     CALL histwrite(nid_tra, "d_tr_cv_"//tnom(it+2), itau_w, &
510     gr_phy_write_3d(d_tr_cv(:, :, it)))
511     CALL histwrite(nid_tra, "d_tr_cl_"//tnom(it+2), itau_w, &
512     gr_phy_write_3d(d_tr_cl(:, :, it)))
513 guez 3 ENDDO
514    
515 guez 20 CALL histwrite(nid_tra, "pplay", itau_w, gr_phy_write_3d(pplay))
516     CALL histwrite(nid_tra, "t", itau_w, gr_phy_write_3d(t_seri))
517 guez 3
518     if (ok_sync) then
519     call histsync(nid_tra)
520     endif
521    
522     end subroutine write_histrac
523    
524     END SUBROUTINE phytrac
525    
526     end module phytrac_m

  ViewVC Help
Powered by ViewVC 1.1.21