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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 13 - (hide annotations)
Fri Jul 25 19:59:34 2008 UTC (15 years, 9 months ago) by guez
Original Path: trunk/libf/phylmd/phytrac.f90
File size: 29249 byte(s)
-- Minor change of behaviour:

"etat0" does not compute "rugsrel" nor "radpas". Deleted arguments
"radpas" and "rugsrel" of "phyredem". Deleted argument "rugsrel" of
"phyetat0". "startphy.nc" does not contain the variable "RUGSREL". In
"physiq", "rugoro" is set to 0 if not "ok_orodr". The whole program
"etat0_lim" does not use "clesphys2".

-- Minor modification of input/output:

Created subroutine "read_clesphys2". Variables of "clesphys2" are read
in "read_clesphys2" instead of "conf_gcm". "printflag" does not print
variables of "clesphys2".

-- Should not change any result at run time:

References to module "numer_rec" instead of individual modules of
"Numer_rec_Lionel".

Deleted argument "clesphy0" of "calfis", "physiq", "conf_gcm",
"leapfrog", "phyetat0". Deleted variable "clesphy0" in
"gcm". "phyetat0" does not modify variables of "clesphys2".

The program unit "gcm" does not modify "itau_phy".

Added some "intent" attributes.

"regr11_lint" does not call "polint".

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     nqmax, 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     pctsrf, frac_impa, frac_nucl, presnivs, pphis, &
16 guez 12 pphi, albsol, rh, cldfra, rneb, diafra, cldliq, itop_con, &
17 guez 3 ibas_con, pmflxr, pmflxs, prfl, psfl, da, phi, mp, upwd, dnwd, tr_seri)
18    
19     ! From phylmd/phytrac.F, version 1.15 2006/02/21 08:08:30
20    
21     ! Authors : Frédéric Hourdin, Abderrahmane Idelkadi, Marie-Alice
22     ! 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     use dimens_m, only: iim, jjm, llm
34     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    
45     ! Arguments:
46    
47     ! EN ENTREE:
48    
49     ! divers:
50    
51     logical, intent(in):: rnpb
52    
53     integer, intent(in):: nqmax
54     ! (nombre de traceurs auxquels on applique la physique)
55    
56 guez 7 integer, intent(in):: itap ! number of calls to "physiq"
57     integer, intent(in):: lmt_pas ! number of time steps of "physics" per day
58 guez 3 integer, intent(in):: julien !jour julien, 1 <= julien <= 360
59     integer itop_con(klon)
60     integer ibas_con(klon)
61     real, intent(in):: gmtime ! heure de la journée en fraction de jour
62 guez 7 real, intent(in):: pdtphys ! pas d'integration pour la physique (s)
63 guez 3 real, intent(in):: t_seri(klon, llm) ! temperature, in K
64    
65     real tr_seri(klon, llm, nbtr)
66     ! (mass fractions of tracers, excluding water, at mid-layers)
67    
68     real u(klon, llm)
69     real v(klon, llm)
70     real rh(klon, llm) ! humidite relative
71     real cldliq(klon, llm) ! eau liquide nuageuse
72     real cldfra(klon, llm) ! fraction nuageuse (tous les nuages)
73    
74     real diafra(klon, llm)
75     ! (fraction nuageuse (convection ou stratus artificiels))
76    
77     real rneb(klon, llm) ! fraction nuageuse (grande echelle)
78     real albsol(klon) ! albedo surface
79    
80     real, intent(in):: paprs(klon, llm+1)
81     ! (pression pour chaque inter-couche, en Pa)
82    
83 guez 10 real, intent(in):: pplay(klon, llm)
84     ! (pression pour le mileu de chaque couche, en Pa)
85    
86 guez 3 real pphi(klon, llm) ! geopotentiel
87     real pphis(klon)
88     REAL, intent(in):: presnivs(llm)
89 guez 7 logical, intent(in):: firstcal ! first call to "calfis"
90 guez 3 logical, intent(in):: lafin ! fin de la physique
91    
92     integer nsplit
93     REAL pmflxr(klon, llm+1), pmflxs(klon, llm+1) !--lessivage convection
94     REAL prfl(klon, llm+1), psfl(klon, llm+1) !--lessivage large-scale
95    
96     ! convection:
97    
98     REAL pmfu(klon, llm) ! flux de masse dans le panache montant
99     REAL pmfd(klon, llm) ! flux de masse dans le panache descendant
100     REAL pen_u(klon, llm) ! flux entraine dans le panache montant
101    
102     ! thermiques:
103    
104     real fm_therm(klon, llm+1), entr_therm(klon, llm)
105    
106     REAL pde_u(klon, llm) ! flux detraine dans le panache montant
107     REAL pen_d(klon, llm) ! flux entraine dans le panache descendant
108     REAL pde_d(klon, llm) ! flux detraine dans le panache descendant
109     ! KE
110     real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
111     REAL upwd(klon, llm) ! saturated updraft mass flux
112     REAL dnwd(klon, llm) ! saturated downdraft mass flux
113    
114     ! Couche limite:
115    
116     REAL coefh(klon, llm) ! coeff melange CL
117     REAL yu1(klon) ! vents au premier niveau
118     REAL yv1(klon) ! vents au premier niveau
119    
120     ! Lessivage:
121    
122     ! pour le ON-LINE
123    
124     REAL frac_impa(klon, llm) ! fraction d'aerosols impactes
125     REAL frac_nucl(klon, llm) ! fraction d'aerosols nuclees
126    
127     ! Arguments necessaires pour les sources et puits de traceur:
128    
129     real ftsol(klon, nbsrf) ! Temperature du sol (surf)(Kelvin)
130     real pctsrf(klon, nbsrf) ! Pourcentage de sol f(nature du sol)
131    
132     real pftsol1(klon), pftsol2(klon), pftsol3(klon), pftsol4(klon)
133     real ppsrf1(klon), ppsrf2(klon), ppsrf3(klon), ppsrf4(klon)
134    
135     ! VARIABLES LOCALES TRACEURS
136    
137     ! Sources et puits des traceurs:
138    
139     ! Pour l'instant seuls les cas du rn et du pb ont ete envisages.
140    
141     REAL source(klon) ! a voir lorsque le flux est prescrit
142     !
143     ! Pour la source de radon et son reservoir de sol
144    
145     REAL, save:: trs(klon, nbtr) ! Concentration de radon dans le sol
146    
147     REAL masktr(klon, nbtr) ! Masque reservoir de sol traceur
148     ! Masque de l'echange avec la surface
149     ! (1 = reservoir) ou (possible => 1 )
150     SAVE masktr
151     REAL fshtr(klon, nbtr) ! Flux surfacique dans le reservoir de sol
152     SAVE fshtr
153     REAL hsoltr(nbtr) ! Epaisseur equivalente du reservoir de sol
154     SAVE hsoltr
155     REAL tautr(nbtr) ! Constante de decroissance radioactive
156     SAVE tautr
157     REAL vdeptr(nbtr) ! Vitesse de depot sec dans la couche Brownienne
158     SAVE vdeptr
159     REAL scavtr(nbtr) ! Coefficient de lessivage
160     SAVE scavtr
161    
162     CHARACTER itn
163     INTEGER, save:: nid_tra
164    
165     ! nature du traceur
166    
167     logical aerosol(nbtr) ! Nature du traceur
168     ! ! aerosol(it) = true => aerosol
169     ! ! aerosol(it) = false => gaz
170     logical clsol(nbtr) ! couche limite sol calculée
171     logical radio(nbtr) ! décroisssance radioactive
172     save aerosol, clsol, radio
173    
174     ! convection tiedtke
175     INTEGER i, k, it
176     REAL delp(klon, llm)
177    
178     ! Variables liees a l'ecriture de la bande histoire physique
179    
180     ! Variables locales pour effectuer les appels en serie
181    
182     REAL d_tr(klon, llm), d_trs(klon) ! tendances de traceurs
183     REAL d_tr_cl(klon, llm, nbtr) ! tendance de traceurs couche limite
184     REAL d_tr_cv(klon, llm, nbtr) ! tendance de traceurs conv pour chq traceur
185     REAL d_tr_th(klon, llm, nbtr) ! la tendance des thermiques
186     REAL d_tr_dec(klon, llm, 2) ! la tendance de la decroissance
187     ! ! radioactive du rn - > pb
188     REAL d_tr_lessi_impa(klon, llm, nbtr) ! la tendance du lessivage
189     ! ! par impaction
190     REAL d_tr_lessi_nucl(klon, llm, nbtr) ! la tendance du lessivage
191     ! ! par nucleation
192     REAL flestottr(klon, llm, nbtr) ! flux de lessivage
193     ! ! dans chaque couche
194    
195     real zmasse(klon, llm)
196     ! (column-density of mass of air in a layer, in kg m-2)
197    
198     real ztra_th(klon, llm)
199    
200     character(len=20) modname
201     character(len=80) abort_message
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     modname='phytrac'
213    
214 guez 7 if (firstcal) then
215 guez 3 print *, 'phytrac: pdtphys = ', pdtphys
216     PRINT *, 'Fréquence de sortie des traceurs : ecrit_tra = ', ecrit_tra
217     if (nbtr < nqmax) then
218     abort_message='See above'
219     call abort_gcm(modname, abort_message, 1)
220     endif
221     inirnpb=rnpb
222    
223     ! Initialisation des sorties :
224     call ini_histrac(nid_tra, pdtphys, presnivs, nqmax, lessivage)
225    
226     ! Initialisation de certaines variables pour le radon et le plomb
227     ! Initialisation du traceur dans le sol (couche limite radonique)
228     trs(:, :) = 0.
229    
230     open (unit=99, file='starttrac', status='old', err=999, &
231     form='formatted')
232     read(unit=99, fmt=*) (trs(i, 1), i=1, klon)
233     999 continue
234     close(unit=99)
235    
236     ! Initialisation de la fraction d'aerosols lessivee
237    
238     d_tr_lessi_impa(:, :, :) = 0.
239     d_tr_lessi_nucl(:, :, :) = 0.
240    
241     ! Initialisation de la nature des traceurs
242    
243     DO it = 1, nqmax
244     aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
245     radio(it) = .FALSE. ! par défaut pas de passage par "radiornpb"
246     clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
247     ENDDO
248     ENDIF
249    
250     ! Initialisation du traceur dans le sol (couche limite radonique)
251     if (inirnpb) THEN
252    
253     radio(1)= .true.
254     radio(2)= .true.
255     clsol(1)= .true.
256     clsol(2)= .true.
257     aerosol(2) = .TRUE. ! le Pb est un aerosol
258    
259     call initrrnpb(ftsol, pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, &
260     scavtr)
261     inirnpb=.false.
262     endif
263    
264     do i=1, klon
265     pftsol1(i) = ftsol(i, 1)
266     pftsol2(i) = ftsol(i, 2)
267     pftsol3(i) = ftsol(i, 3)
268     pftsol4(i) = ftsol(i, 4)
269    
270     ppsrf1(i) = pctsrf(i, 1)
271     ppsrf2(i) = pctsrf(i, 2)
272     ppsrf3(i) = pctsrf(i, 3)
273     ppsrf4(i) = pctsrf(i, 4)
274    
275     enddo
276    
277     ! Calcul de l'effet de la convection
278    
279     if (convection) then
280     DO it=1, nqmax
281     if (iflag_con.eq.2) then
282     ! tiedke
283     CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
284 guez 10 paprs, tr_seri(1, 1, it), d_tr_cv(1, 1, it))
285 guez 3 else if (iflag_con.eq.3) then
286     ! KE
287 guez 10 call cvltr(pdtphys, da, phi, mp, paprs, &
288 guez 3 tr_seri(1, 1, it), upwd, dnwd, d_tr_cv(1, 1, it))
289     endif
290    
291     DO k = 1, llm
292     DO i = 1, klon
293     tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cv(i, k, it)
294     ENDDO
295     ENDDO
296     WRITE(unit=itn, fmt='(i1)') it
297     CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, &
298     'convection, tracer index = ' // itn)
299     ENDDO
300     endif
301    
302     forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
303    
304     ! Calcul de l'effet des thermiques
305    
306     do it=1, nqmax
307     do k=1, llm
308     do i=1, klon
309     d_tr_th(i, k, it)=0.
310     tr_seri(i, k, it)=max(tr_seri(i, k, it), 0.)
311     tr_seri(i, k, it)=min(tr_seri(i, k, it), 1.e10)
312     enddo
313     enddo
314     enddo
315    
316     if (iflag_thermals > 0) then
317     nsplit=10
318     DO it=1, nqmax
319     do isplit=1, nsplit
320     ! Thermiques
321     call dqthermcell(klon, llm, pdtphys/nsplit &
322     , fm_therm, entr_therm, zmasse &
323     , tr_seri(1:klon, 1:llm, it), d_tr, ztra_th)
324    
325     do k=1, llm
326     do i=1, klon
327     d_tr(i, k)=pdtphys*d_tr(i, k)/nsplit
328     d_tr_th(i, k, it)=d_tr_th(i, k, it)+d_tr(i, k)
329     tr_seri(i, k, it)=max(tr_seri(i, k, it)+d_tr(i, k), 0.)
330     enddo
331     enddo
332     enddo
333     ENDDO
334     endif
335    
336     ! Calcul de l'effet de la couche limite
337    
338     if (couchelimite) then
339    
340     DO k = 1, llm
341     DO i = 1, klon
342     delp(i, k) = paprs(i, k)-paprs(i, k+1)
343     ENDDO
344     ENDDO
345    
346     ! MAF modif pour tenir compte du cas rnpb + traceur
347     DO it=1, nqmax
348     if (clsol(it)) then
349     ! couche limite avec quantite dans le sol calculee
350     CALL cltracrn(it, pdtphys, yu1, yv1, &
351     coefh, t_seri, ftsol, pctsrf, &
352     tr_seri(1, 1, it), trs(1, it), &
353     paprs, pplay, delp, &
354     masktr(1, it), fshtr(1, it), hsoltr(it), &
355     tautr(it), vdeptr(it), &
356     rlat, &
357     d_tr_cl(1, 1, it), d_trs)
358     DO k = 1, llm
359     DO i = 1, klon
360     tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
361     ENDDO
362     ENDDO
363    
364     ! Traceur ds sol
365    
366     DO i = 1, klon
367     trs(i, it) = trs(i, it) + d_trs(i)
368     END DO
369     else ! couche limite avec flux prescrit
370     !MAF provisoire source / traceur a creer
371     DO i=1, klon
372     source(i) = 0.0 ! pas de source, pour l'instant
373     ENDDO
374    
375     CALL cltrac(pdtphys, coefh, t_seri, &
376     tr_seri(1, 1, it), source, &
377     paprs, pplay, delp, &
378     d_tr_cl(1, 1, it))
379     DO k = 1, llm
380     DO i = 1, klon
381     tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
382     ENDDO
383     ENDDO
384     endif
385     ENDDO
386    
387     endif ! couche limite
388    
389     ! Calcul de l'effet du puits radioactif
390    
391     ! MAF il faudrait faire une modification pour passer dans radiornpb
392     ! si radio=true mais pour l'instant radiornpb propre au cas rnpb
393     if (rnpb) then
394     d_tr_dec(:, :, :) = radiornpb(tr_seri, pdtphys, tautr)
395     DO it=1, nqmax
396     if (radio(it)) then
397     tr_seri(:, :, it) = tr_seri(:, :, it) + d_tr_dec(:, :, it)
398     WRITE(unit=itn, fmt='(i1)') it
399     CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, 'puits rn it='//itn)
400     endif
401     ENDDO
402     endif ! rnpb decroissance radioactive
403    
404 guez 6 if (nqmax >= 3) then
405     ! Ozone as a tracer:
406 guez 7 if (mod(itap - 1, lmt_pas) == 0) then
407     ! Once per day, update the coefficients for ozone chemistry:
408     call regr_pr_comb_coefoz(julien)
409     end if
410 guez 6 call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
411     end if
412 guez 3
413     ! Calcul de l'effet de la precipitation
414    
415     IF (lessivage) THEN
416     d_tr_lessi_nucl(:, :, :) = 0.
417     d_tr_lessi_impa(:, :, :) = 0.
418     flestottr(:, :, :) = 0.
419    
420     ! tendance des aerosols nuclees et impactes
421    
422     DO it = 1, nqmax
423     IF (aerosol(it)) THEN
424     DO k = 1, llm
425     DO i = 1, klon
426     d_tr_lessi_nucl(i, k, it) = d_tr_lessi_nucl(i, k, it) + &
427     ( 1 - frac_nucl(i, k) )*tr_seri(i, k, it)
428     d_tr_lessi_impa(i, k, it) = d_tr_lessi_impa(i, k, it) + &
429     ( 1 - frac_impa(i, k) )*tr_seri(i, k, it)
430     ENDDO
431     ENDDO
432     ENDIF
433     ENDDO
434    
435     ! Mises a jour des traceurs + calcul des flux de lessivage
436     ! Mise a jour due a l'impaction et a la nucleation
437    
438     DO it = 1, nqmax
439     IF (aerosol(it)) THEN
440     DO k = 1, llm
441     DO i = 1, klon
442     tr_seri(i, k, it)=tr_seri(i, k, it) &
443     *frac_impa(i, k)*frac_nucl(i, k)
444     ENDDO
445     ENDDO
446     ENDIF
447     ENDDO
448    
449     ! Flux lessivage total
450    
451     DO it = 1, nqmax
452     DO k = 1, llm
453     DO i = 1, klon
454     flestottr(i, k, it) = flestottr(i, k, it) - &
455     ( d_tr_lessi_nucl(i, k, it) + &
456     d_tr_lessi_impa(i, k, it) ) * &
457     ( paprs(i, k)-paprs(i, k+1) ) / &
458     (RG * pdtphys)
459     ENDDO
460     ENDDO
461     ENDDO
462     ENDIF
463    
464     ! Ecriture des sorties
465 guez 7 call write_histrac(lessivage, nqmax, itap, nid_tra)
466 guez 3
467     if (lafin) then
468     print *, "C'est la fin de la physique."
469     open (unit=99, file='restarttrac', form='formatted')
470     do i=1, klon
471     write(unit=99, fmt=*) trs(i, 1)
472     enddo
473     PRINT *, 'Ecriture du fichier restarttrac'
474     close(99)
475     endif
476    
477     contains
478    
479 guez 7 subroutine write_histrac(lessivage, nqmax, itap, nid_tra)
480 guez 3
481     ! From phylmd/write_histrac.h, version 1.9 2006/02/21 08:08:30
482    
483     use dimens_m, only: iim, jjm, llm
484     use ioipsl, only: histwrite, histsync
485     use temps, only: itau_phy
486     use advtrac_m, only: tnom
487     use comgeomphy, only: airephy
488     use dimphy, only: klon
489    
490     logical, intent(in):: lessivage
491    
492     integer, intent(in):: nqmax
493     ! (nombre de traceurs auxquels on applique la physique)
494    
495 guez 7 integer, intent(in):: itap ! number of calls to "physiq"
496 guez 3 integer, intent(in):: nid_tra
497    
498     ! Variables local to the procedure:
499     INTEGER ndex2d(iim*(jjm+1)), ndex3d(iim*(jjm+1)*llm)
500     integer it
501     integer itau_w ! pas de temps ecriture
502     REAL zx_tmp_2d(iim, jjm+1), zx_tmp_3d(iim, jjm+1, llm)
503     logical, parameter:: ok_sync = .true.
504    
505     !-----------------------------------------------------
506    
507     ndex2d = 0
508     ndex3d = 0
509 guez 7 itau_w = itau_phy + itap
510 guez 3
511     CALL gr_fi_ecrit(1, klon, iim, jjm+1, pphis, zx_tmp_2d)
512     CALL histwrite(nid_tra, "phis", itau_w, zx_tmp_2d, iim*(jjm+1), ndex2d)
513    
514     CALL gr_fi_ecrit(1, klon, iim, jjm+1, airephy, zx_tmp_2d)
515     CALL histwrite(nid_tra, "aire", itau_w, zx_tmp_2d, iim*(jjm+1), ndex2d)
516    
517     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, zmasse, zx_tmp_3d)
518     CALL histwrite(nid_tra, "zmasse", itau_w, zx_tmp_3d, iim*(jjm+1)*llm, &
519     ndex3d)
520    
521     DO it=1, nqmax
522     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, tr_seri(1, 1, it), zx_tmp_3d)
523     CALL histwrite(nid_tra, tnom(it+2), itau_w, zx_tmp_3d, &
524     iim*(jjm+1)*llm, ndex3d)
525     if (lessivage) THEN
526     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, flestottr(1, 1, it), &
527     zx_tmp_3d)
528     CALL histwrite(nid_tra, "fl"//tnom(it+2), itau_w, zx_tmp_3d, &
529     iim*(jjm+1)*llm, ndex3d)
530     endif
531    
532     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_th(1, 1, it), zx_tmp_3d)
533     CALL histwrite(nid_tra, "d_tr_th_"//tnom(it+2), itau_w, zx_tmp_3d, &
534     iim*(jjm+1)*llm, ndex3d)
535     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_cv(1, 1, it), zx_tmp_3d)
536     CALL histwrite(nid_tra, "d_tr_cv_"//tnom(it+2), itau_w, zx_tmp_3d, &
537     iim*(jjm+1)*llm, ndex3d)
538     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_cl(1, 1, it), zx_tmp_3d)
539     CALL histwrite(nid_tra, "d_tr_cl_"//tnom(it+2), itau_w, zx_tmp_3d, &
540     iim*(jjm+1)*llm, ndex3d)
541     ENDDO
542    
543     CALL gr_fi_ecrit(1, klon, iim, jjm+1, yu1, zx_tmp_2d)
544     CALL histwrite(nid_tra, "pyu1", itau_w, zx_tmp_2d, &
545     iim*(jjm+1), ndex2d)
546    
547     CALL gr_fi_ecrit(1, klon, iim, jjm+1, yv1, zx_tmp_2d)
548     CALL histwrite(nid_tra, "pyv1", itau_w, zx_tmp_2d, &
549     iim*(jjm+1), ndex2d)
550    
551     CALL gr_fi_ecrit(1, klon, iim, jjm+1, pftsol1, zx_tmp_2d)
552     CALL histwrite(nid_tra, "ftsol1", itau_w, zx_tmp_2d, &
553     iim*(jjm+1), ndex2d)
554    
555     CALL gr_fi_ecrit(1, klon, iim, jjm+1, pftsol2, zx_tmp_2d)
556     CALL histwrite(nid_tra, "ftsol2", itau_w, zx_tmp_2d, &
557     iim*(jjm+1), ndex2d)
558    
559     CALL gr_fi_ecrit(1, klon, iim, jjm+1, pftsol3, zx_tmp_2d)
560     CALL histwrite(nid_tra, "ftsol3", itau_w, zx_tmp_2d, &
561     iim*(jjm+1), ndex2d)
562    
563     CALL gr_fi_ecrit(1, klon, iim, jjm+1, pftsol4, zx_tmp_2d)
564     CALL histwrite(nid_tra, "ftsol4", itau_w, zx_tmp_2d, &
565     iim*(jjm+1), ndex2d)
566    
567     CALL gr_fi_ecrit(1, klon, iim, jjm+1, ppsrf1, zx_tmp_2d)
568     CALL histwrite(nid_tra, "psrf1", itau_w, zx_tmp_2d, &
569     iim*(jjm+1), ndex2d)
570    
571     CALL gr_fi_ecrit(1, klon, iim, jjm+1, ppsrf2, zx_tmp_2d)
572     CALL histwrite(nid_tra, "psrf2", itau_w, zx_tmp_2d, &
573     iim*(jjm+1), ndex2d)
574    
575     CALL gr_fi_ecrit(1, klon, iim, jjm+1, ppsrf3, zx_tmp_2d)
576     CALL histwrite(nid_tra, "psrf3", itau_w, zx_tmp_2d, &
577     iim*(jjm+1), ndex2d)
578    
579     CALL gr_fi_ecrit(1, klon, iim, jjm+1, ppsrf4, zx_tmp_2d)
580     CALL histwrite(nid_tra, "psrf4", itau_w, zx_tmp_2d, &
581     iim*(jjm+1), ndex2d)
582     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pplay, zx_tmp_3d)
583     CALL histwrite(nid_tra, "pplay", itau_w, zx_tmp_3d, &
584     iim*(jjm+1)*llm, ndex3d)
585    
586     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, t_seri, zx_tmp_3d)
587     CALL histwrite(nid_tra, "t", itau_w, zx_tmp_3d, &
588     iim*(jjm+1)*llm, ndex3d)
589     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pmfu, zx_tmp_3d)
590     CALL histwrite(nid_tra, "mfu", itau_w, zx_tmp_3d, &
591     iim*(jjm+1)*llm, ndex3d)
592     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pmfd, zx_tmp_3d)
593     CALL histwrite(nid_tra, "mfd", itau_w, zx_tmp_3d, &
594     iim*(jjm+1)*llm, ndex3d)
595     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pen_u, zx_tmp_3d)
596     CALL histwrite(nid_tra, "en_u", itau_w, zx_tmp_3d, &
597     iim*(jjm+1)*llm, ndex3d)
598     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pen_d, zx_tmp_3d)
599     CALL histwrite(nid_tra, "en_d", itau_w, zx_tmp_3d, &
600     iim*(jjm+1)*llm, ndex3d)
601     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pde_d, zx_tmp_3d)
602     CALL histwrite(nid_tra, "de_d", itau_w, zx_tmp_3d, &
603     iim*(jjm+1)*llm, ndex3d)
604     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pde_u, zx_tmp_3d)
605     CALL histwrite(nid_tra, "de_u", itau_w, zx_tmp_3d, &
606     iim*(jjm+1)*llm, ndex3d)
607     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, coefh, zx_tmp_3d)
608     CALL histwrite(nid_tra, "coefh", itau_w, zx_tmp_3d, &
609     iim*(jjm+1)*llm, ndex3d)
610    
611     ! abder
612    
613     if (ok_sync) then
614     call histsync(nid_tra)
615     endif
616    
617     end subroutine write_histrac
618    
619     END SUBROUTINE phytrac
620    
621     !*************************************************
622    
623     subroutine ini_histrac(nid_tra, pdtphys, presnivs, nqmax, lessivage)
624    
625     ! From phylmd/ini_histrac.h, version 1.10 2006/02/21 08:08:30
626    
627     use dimens_m, only: iim, jjm, llm
628     use ioipsl, only: ymds2ju, histbeg_totreg, histvert, histdef, histend
629     use temps, only: annee_ref, day_ref, itau_phy
630     use advtrac_m, only: niadv, tnom, ttext
631     use dimphy, only: klon
632     use clesphys, only: ecrit_tra
633     use grid_change, only: gr_phy_write
634     use phyetat0_m, only: rlon, rlat
635    
636     INTEGER, intent(out):: nid_tra
637     real, intent(in):: pdtphys ! pas d'integration pour la physique (s)
638     REAL, intent(in):: presnivs(:)
639    
640     integer, intent(in):: nqmax
641     ! (nombre de traceurs auxquels on applique la physique)
642    
643     logical, intent(in):: lessivage
644    
645     ! Variables local to the procedure:
646    
647     REAL zjulian
648     REAL zx_lat(iim, jjm+1)
649     INTEGER nhori, nvert
650     REAL zsto, zout
651     integer it, iq, iiq
652    
653     !---------------------------------------------------------
654    
655     CALL ymds2ju(annee_ref, month=1, day=day_ref, sec=0.0, julian=zjulian)
656     zx_lat(:, :) = gr_phy_write(rlat)
657     CALL histbeg_totreg("histrac", iim, rlon(2:iim+1), jjm+1, zx_lat(1, :), &
658     1, iim, 1, jjm+1, itau_phy, zjulian, pdtphys, nhori, nid_tra)
659     CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb", llm, &
660     presnivs, nvert)
661    
662     zsto = pdtphys
663     zout = pdtphys * REAL(ecrit_tra)
664    
665     CALL histdef(nid_tra, "phis", "Surface geop. height", "-", &
666     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
667     "once", zsto, zout)
668     CALL histdef(nid_tra, "aire", "Grid area", "-", &
669     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
670     "once", zsto, zout)
671     CALL histdef(nid_tra, "zmasse", "column density of air in cell", &
672     "kg m-2", iim, jjm + 1, nhori, llm, 1, llm, nvert, 32, "ave(X)", &
673     zsto, zout)
674    
675     DO it=1, nqmax
676     ! champ 2D
677     iq=it+2
678     iiq=niadv(iq)
679     CALL histdef(nid_tra, tnom(iq), ttext(iiq), "U/kga", iim, jjm+1, &
680     nhori, llm, 1, llm, nvert, 32, "ave(X)", zsto, zout)
681     if (lessivage) THEN
682     CALL histdef(nid_tra, "fl"//tnom(iq), "Flux "//ttext(iiq), &
683     "U/m2/s", iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
684     "ave(X)", zsto, zout)
685     endif
686    
687     !---Ajout Olivia
688     CALL histdef(nid_tra, "d_tr_th_"//tnom(iq), &
689     "tendance thermique"// ttext(iiq), "?", &
690     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
691     "ave(X)", zsto, zout)
692     CALL histdef(nid_tra, "d_tr_cv_"//tnom(iq), &
693     "tendance convection"// ttext(iiq), "?", &
694     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
695     "ave(X)", zsto, zout)
696     CALL histdef(nid_tra, "d_tr_cl_"//tnom(iq), &
697     "tendance couche limite"// ttext(iiq), "?", &
698     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
699     "ave(X)", zsto, zout)
700     !---fin Olivia
701    
702     ENDDO
703    
704     CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-", &
705     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
706     "inst(X)", zout, zout)
707    
708     CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-", &
709     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
710     "inst(X)", zout, zout)
711     CALL histdef(nid_tra, "psrf1", "nature sol", "-", &
712     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
713     "inst(X)", zout, zout)
714     CALL histdef(nid_tra, "psrf2", "nature sol", "-", &
715     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
716     "inst(X)", zout, zout)
717     CALL histdef(nid_tra, "psrf3", "nature sol", "-", &
718     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
719     "inst(X)", zout, zout)
720     CALL histdef(nid_tra, "psrf4", "nature sol", "-", &
721     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
722     "inst(X)", zout, zout)
723     CALL histdef(nid_tra, "ftsol1", "temper sol", "-", &
724     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
725     "inst(X)", zout, zout)
726     CALL histdef(nid_tra, "ftsol2", "temper sol", "-", &
727     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
728     "inst(X)", zout, zout)
729     CALL histdef(nid_tra, "ftsol3", "temper sol", "-", &
730     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
731     "inst(X)", zout, zout)
732     CALL histdef(nid_tra, "ftsol4", "temper sol", "-", &
733     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
734     "inst(X)", zout, zout)
735     CALL histdef(nid_tra, "pplay", "flux u mont", "-", &
736     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
737     "inst(X)", zout, zout)
738     CALL histdef(nid_tra, "t", "flux u mont", "-", &
739     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
740     "inst(X)", zout, zout)
741     CALL histdef(nid_tra, "mfu", "flux u mont", "-", &
742     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
743     "ave(X)", zsto, zout)
744     CALL histdef(nid_tra, "mfd", "flux u decen", "-", &
745     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
746     "ave(X)", zsto, zout)
747     CALL histdef(nid_tra, "en_u", "flux u mont", "-", &
748     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
749     "ave(X)", zsto, zout)
750     CALL histdef(nid_tra, "en_d", "flux u mont", "-", &
751     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
752     "ave(X)", zsto, zout)
753     CALL histdef(nid_tra, "de_d", "flux u mont", "-", &
754     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
755     "ave(X)", zsto, zout)
756     CALL histdef(nid_tra, "de_u", "flux u decen", "-", &
757     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
758     "ave(X)", zsto, zout)
759     CALL histdef(nid_tra, "coefh", "turbulent coef", "-", &
760     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
761     "ave(X)", zsto, zout)
762    
763     CALL histend(nid_tra)
764    
765     end subroutine ini_histrac
766    
767     !*************************************************
768    
769     function radiornpb(tr_seri, pdtphys, tautr)
770    
771     ! From phylmd/radiornpb.F, v 1.2 2005/05/25 13:10:09
772    
773     ! Auteurs: AA + CG (LGGE/CNRS) Date 24-06-94
774     ! Objet: Decroissance radioactive d'un traceur dans l'atmosphere
775     !G 24 06 94 : Pour un traceur, le radon
776     !G 16 12 94 : Plus un 2eme traceur, le 210Pb. Le radon decroit en plomb.
777    
778     ! Le pas de temps "pdtphys" est supposé beaucoup plus petit que la
779     ! constante de temps de décroissance.
780    
781     use dimens_m, only: llm
782     use dimphy, only: klon, nbtr
783 guez 13 use numer_rec, only: assert
784 guez 3
785     IMPLICIT none
786    
787     REAL, intent(in):: tr_seri(:, :, :), pdtphys, tautr(:)
788     real radiornpb(klon, llm, 2)
789    
790     ! Variable local to the procedure:
791     INTEGER it
792    
793     !-----------------------------------------------
794    
795     call assert(shape(tr_seri) == (/klon, llm, nbtr/), "radiornpb tr_seri")
796     call assert(size(tautr) == nbtr, "radiornpb tautr")
797    
798     DO it = 1, 2
799     IF (tautr(it) > 0.) THEN
800     radiornpb(:, :, it) = - tr_seri(:, :, it) * pdtphys / tautr(it)
801     ELSE
802     radiornpb(:, :, it) = 0.
803     END IF
804     END DO
805    
806     !G161294 : Cas particulier radon 1 => plomb 2
807     radiornpb(:, :, 2) = radiornpb(:, :, 2) - radiornpb(:, :, 1)
808    
809     END function radiornpb
810    
811     !*************************************************
812    
813     SUBROUTINE minmaxqfi(zq, qmin, qmax, comment)
814    
815     ! From phylmd/minmaxqfi.F, version 1.1.1.1 2004/05/19 12:53:09
816    
817     use dimens_m, only: llm
818     use dimphy, only: klon
819 guez 13 use numer_rec, only: assert
820 guez 3
821     IMPLICIT none
822    
823     real, intent(in):: zq(:, :), qmin, qmax
824     CHARACTER(len=*), intent(in):: comment
825    
826     ! Variables local to the procedure:
827    
828     INTEGER jadrs(klon), jbad, k, i
829    
830     !---------------------------------
831    
832     call assert(shape(zq) == (/klon, llm/), "minmaxqfi")
833    
834     DO k = 1, llm
835     jbad = 0
836     DO i = 1, klon
837     IF (zq(i, k) > qmax .OR. zq(i, k) < qmin) THEN
838     jbad = jbad + 1
839     jadrs(jbad) = i
840     ENDIF
841     ENDDO
842     IF (jbad > 0) THEN
843     PRINT *, comment
844     DO i = 1, jbad
845     PRINT *, "zq(", jadrs(i), ", ", k, ") = ", zq(jadrs(i), k)
846     ENDDO
847     ENDIF
848     ENDDO
849    
850     end SUBROUTINE minmaxqfi
851    
852     end module phytrac_m

  ViewVC Help
Powered by ViewVC 1.1.21