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

Annotation of /trunk/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21