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

Annotation of /trunk/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6 - (hide annotations)
Tue Mar 4 14:00:42 2008 UTC (16 years, 2 months ago) by guez
Original Path: trunk/libf/phylmd/phytrac.f90
File size: 29099 byte(s)
Removed test on coefoz_LMDZ in gcm.sh.
Added test on nbtr in ini_histhf3d and write_histhf3d.
Added test on nqmax in phytrac.
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 guez 6 if (nqmax >= 3) then
247     ! Get the parameters for ozone chemistry:
248     call read_coefoz
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     do i=1, klon
267     pftsol1(i) = ftsol(i, 1)
268     pftsol2(i) = ftsol(i, 2)
269     pftsol3(i) = ftsol(i, 3)
270     pftsol4(i) = ftsol(i, 4)
271    
272     ppsrf1(i) = pctsrf(i, 1)
273     ppsrf2(i) = pctsrf(i, 2)
274     ppsrf3(i) = pctsrf(i, 3)
275     ppsrf4(i) = pctsrf(i, 4)
276    
277     enddo
278    
279     ! Calcul de l'effet de la convection
280    
281     if (convection) then
282     DO it=1, nqmax
283     if (iflag_con.eq.2) then
284     ! tiedke
285     CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
286     pplay, paprs, tr_seri(1, 1, it), d_tr_cv(1, 1, it))
287     else if (iflag_con.eq.3) then
288     ! KE
289     call cvltr(pdtphys, da, phi, mp, paprs, pplay, &
290     tr_seri(1, 1, it), upwd, dnwd, d_tr_cv(1, 1, it))
291     endif
292    
293     DO k = 1, llm
294     DO i = 1, klon
295     tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cv(i, k, it)
296     ENDDO
297     ENDDO
298     WRITE(unit=itn, fmt='(i1)') it
299     CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, &
300     'convection, tracer index = ' // itn)
301     ENDDO
302     endif
303    
304     forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
305    
306     ! Calcul de l'effet des thermiques
307    
308     do it=1, nqmax
309     do k=1, llm
310     do i=1, klon
311     d_tr_th(i, k, it)=0.
312     tr_seri(i, k, it)=max(tr_seri(i, k, it), 0.)
313     tr_seri(i, k, it)=min(tr_seri(i, k, it), 1.e10)
314     enddo
315     enddo
316     enddo
317    
318     if (iflag_thermals > 0) then
319     nsplit=10
320     DO it=1, nqmax
321     do isplit=1, nsplit
322     ! Thermiques
323     call dqthermcell(klon, llm, pdtphys/nsplit &
324     , fm_therm, entr_therm, zmasse &
325     , tr_seri(1:klon, 1:llm, it), d_tr, ztra_th)
326    
327     do k=1, llm
328     do i=1, klon
329     d_tr(i, k)=pdtphys*d_tr(i, k)/nsplit
330     d_tr_th(i, k, it)=d_tr_th(i, k, it)+d_tr(i, k)
331     tr_seri(i, k, it)=max(tr_seri(i, k, it)+d_tr(i, k), 0.)
332     enddo
333     enddo
334     enddo
335     ENDDO
336     endif
337    
338     ! Calcul de l'effet de la couche limite
339    
340     if (couchelimite) then
341    
342     DO k = 1, llm
343     DO i = 1, klon
344     delp(i, k) = paprs(i, k)-paprs(i, k+1)
345     ENDDO
346     ENDDO
347    
348     ! MAF modif pour tenir compte du cas rnpb + traceur
349     DO it=1, nqmax
350     if (clsol(it)) then
351     ! couche limite avec quantite dans le sol calculee
352     CALL cltracrn(it, pdtphys, yu1, yv1, &
353     coefh, t_seri, ftsol, pctsrf, &
354     tr_seri(1, 1, it), trs(1, it), &
355     paprs, pplay, delp, &
356     masktr(1, it), fshtr(1, it), hsoltr(it), &
357     tautr(it), vdeptr(it), &
358     rlat, &
359     d_tr_cl(1, 1, it), d_trs)
360     DO k = 1, llm
361     DO i = 1, klon
362     tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
363     ENDDO
364     ENDDO
365    
366     ! Traceur ds sol
367    
368     DO i = 1, klon
369     trs(i, it) = trs(i, it) + d_trs(i)
370     END DO
371     else ! couche limite avec flux prescrit
372     !MAF provisoire source / traceur a creer
373     DO i=1, klon
374     source(i) = 0.0 ! pas de source, pour l'instant
375     ENDDO
376    
377     CALL cltrac(pdtphys, coefh, t_seri, &
378     tr_seri(1, 1, it), source, &
379     paprs, pplay, delp, &
380     d_tr_cl(1, 1, it))
381     DO k = 1, llm
382     DO i = 1, klon
383     tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
384     ENDDO
385     ENDDO
386     endif
387     ENDDO
388    
389     endif ! couche limite
390    
391     ! Calcul de l'effet du puits radioactif
392    
393     ! MAF il faudrait faire une modification pour passer dans radiornpb
394     ! si radio=true mais pour l'instant radiornpb propre au cas rnpb
395     if (rnpb) then
396     d_tr_dec(:, :, :) = radiornpb(tr_seri, pdtphys, tautr)
397     DO it=1, nqmax
398     if (radio(it)) then
399     tr_seri(:, :, it) = tr_seri(:, :, it) + d_tr_dec(:, :, it)
400     WRITE(unit=itn, fmt='(i1)') it
401     CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, 'puits rn it='//itn)
402     endif
403     ENDDO
404     endif ! rnpb decroissance radioactive
405    
406 guez 6 if (nqmax >= 3) then
407     ! Ozone as a tracer:
408     call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
409     end if
410 guez 3
411     ! Calcul de l'effet de la precipitation
412    
413     IF (lessivage) THEN
414     d_tr_lessi_nucl(:, :, :) = 0.
415     d_tr_lessi_impa(:, :, :) = 0.
416     flestottr(:, :, :) = 0.
417    
418     ! tendance des aerosols nuclees et impactes
419    
420     DO it = 1, nqmax
421     IF (aerosol(it)) THEN
422     DO k = 1, llm
423     DO i = 1, klon
424     d_tr_lessi_nucl(i, k, it) = d_tr_lessi_nucl(i, k, it) + &
425     ( 1 - frac_nucl(i, k) )*tr_seri(i, k, it)
426     d_tr_lessi_impa(i, k, it) = d_tr_lessi_impa(i, k, it) + &
427     ( 1 - frac_impa(i, k) )*tr_seri(i, k, it)
428     ENDDO
429     ENDDO
430     ENDIF
431     ENDDO
432    
433     ! Mises a jour des traceurs + calcul des flux de lessivage
434     ! Mise a jour due a l'impaction et a la nucleation
435    
436     DO it = 1, nqmax
437     IF (aerosol(it)) THEN
438     DO k = 1, llm
439     DO i = 1, klon
440     tr_seri(i, k, it)=tr_seri(i, k, it) &
441     *frac_impa(i, k)*frac_nucl(i, k)
442     ENDDO
443     ENDDO
444     ENDIF
445     ENDDO
446    
447     ! Flux lessivage total
448    
449     DO it = 1, nqmax
450     DO k = 1, llm
451     DO i = 1, klon
452     flestottr(i, k, it) = flestottr(i, k, it) - &
453     ( d_tr_lessi_nucl(i, k, it) + &
454     d_tr_lessi_impa(i, k, it) ) * &
455     ( paprs(i, k)-paprs(i, k+1) ) / &
456     (RG * pdtphys)
457     ENDDO
458     ENDDO
459     ENDDO
460     ENDIF
461    
462     ! Ecriture des sorties
463     call write_histrac(lessivage, nqmax, nstep, nid_tra)
464    
465     if (lafin) then
466     print *, "C'est la fin de la physique."
467     open (unit=99, file='restarttrac', form='formatted')
468     do i=1, klon
469     write(unit=99, fmt=*) trs(i, 1)
470     enddo
471     PRINT *, 'Ecriture du fichier restarttrac'
472     close(99)
473     endif
474    
475     contains
476    
477     subroutine write_histrac(lessivage, nqmax, nstep, nid_tra)
478    
479     ! From phylmd/write_histrac.h, version 1.9 2006/02/21 08:08:30
480    
481     use dimens_m, only: iim, jjm, llm
482     use ioipsl, only: histwrite, histsync
483     use temps, only: itau_phy
484     use advtrac_m, only: tnom
485     use comgeomphy, only: airephy
486     use dimphy, only: klon
487    
488     logical, intent(in):: lessivage
489    
490     integer, intent(in):: nqmax
491     ! (nombre de traceurs auxquels on applique la physique)
492    
493     integer, intent(in):: nstep ! appel physique
494     integer, intent(in):: nid_tra
495    
496     ! Variables local to the procedure:
497     INTEGER ndex2d(iim*(jjm+1)), ndex3d(iim*(jjm+1)*llm)
498     integer it
499     integer itau_w ! pas de temps ecriture
500     REAL zx_tmp_2d(iim, jjm+1), zx_tmp_3d(iim, jjm+1, llm)
501     logical, parameter:: ok_sync = .true.
502    
503     !-----------------------------------------------------
504    
505     ndex2d = 0
506     ndex3d = 0
507     itau_w = itau_phy + nstep
508    
509     CALL gr_fi_ecrit(1, klon, iim, jjm+1, pphis, zx_tmp_2d)
510     CALL histwrite(nid_tra, "phis", itau_w, zx_tmp_2d, iim*(jjm+1), ndex2d)
511    
512     CALL gr_fi_ecrit(1, klon, iim, jjm+1, airephy, zx_tmp_2d)
513     CALL histwrite(nid_tra, "aire", itau_w, zx_tmp_2d, iim*(jjm+1), ndex2d)
514    
515     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, zmasse, zx_tmp_3d)
516     CALL histwrite(nid_tra, "zmasse", itau_w, zx_tmp_3d, iim*(jjm+1)*llm, &
517     ndex3d)
518    
519     DO it=1, nqmax
520     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, tr_seri(1, 1, it), zx_tmp_3d)
521     CALL histwrite(nid_tra, tnom(it+2), itau_w, zx_tmp_3d, &
522     iim*(jjm+1)*llm, ndex3d)
523     if (lessivage) THEN
524     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, flestottr(1, 1, it), &
525     zx_tmp_3d)
526     CALL histwrite(nid_tra, "fl"//tnom(it+2), itau_w, zx_tmp_3d, &
527     iim*(jjm+1)*llm, ndex3d)
528     endif
529    
530     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_th(1, 1, it), zx_tmp_3d)
531     CALL histwrite(nid_tra, "d_tr_th_"//tnom(it+2), itau_w, zx_tmp_3d, &
532     iim*(jjm+1)*llm, ndex3d)
533     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_cv(1, 1, it), zx_tmp_3d)
534     CALL histwrite(nid_tra, "d_tr_cv_"//tnom(it+2), itau_w, zx_tmp_3d, &
535     iim*(jjm+1)*llm, ndex3d)
536     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, d_tr_cl(1, 1, it), zx_tmp_3d)
537     CALL histwrite(nid_tra, "d_tr_cl_"//tnom(it+2), itau_w, zx_tmp_3d, &
538     iim*(jjm+1)*llm, ndex3d)
539     ENDDO
540    
541     CALL gr_fi_ecrit(1, klon, iim, jjm+1, yu1, zx_tmp_2d)
542     CALL histwrite(nid_tra, "pyu1", itau_w, zx_tmp_2d, &
543     iim*(jjm+1), ndex2d)
544    
545     CALL gr_fi_ecrit(1, klon, iim, jjm+1, yv1, zx_tmp_2d)
546     CALL histwrite(nid_tra, "pyv1", itau_w, zx_tmp_2d, &
547     iim*(jjm+1), ndex2d)
548    
549     CALL gr_fi_ecrit(1, klon, iim, jjm+1, pftsol1, zx_tmp_2d)
550     CALL histwrite(nid_tra, "ftsol1", itau_w, zx_tmp_2d, &
551     iim*(jjm+1), ndex2d)
552    
553     CALL gr_fi_ecrit(1, klon, iim, jjm+1, pftsol2, zx_tmp_2d)
554     CALL histwrite(nid_tra, "ftsol2", itau_w, zx_tmp_2d, &
555     iim*(jjm+1), ndex2d)
556    
557     CALL gr_fi_ecrit(1, klon, iim, jjm+1, pftsol3, zx_tmp_2d)
558     CALL histwrite(nid_tra, "ftsol3", itau_w, zx_tmp_2d, &
559     iim*(jjm+1), ndex2d)
560    
561     CALL gr_fi_ecrit(1, klon, iim, jjm+1, pftsol4, zx_tmp_2d)
562     CALL histwrite(nid_tra, "ftsol4", itau_w, zx_tmp_2d, &
563     iim*(jjm+1), ndex2d)
564    
565     CALL gr_fi_ecrit(1, klon, iim, jjm+1, ppsrf1, zx_tmp_2d)
566     CALL histwrite(nid_tra, "psrf1", itau_w, zx_tmp_2d, &
567     iim*(jjm+1), ndex2d)
568    
569     CALL gr_fi_ecrit(1, klon, iim, jjm+1, ppsrf2, zx_tmp_2d)
570     CALL histwrite(nid_tra, "psrf2", itau_w, zx_tmp_2d, &
571     iim*(jjm+1), ndex2d)
572    
573     CALL gr_fi_ecrit(1, klon, iim, jjm+1, ppsrf3, zx_tmp_2d)
574     CALL histwrite(nid_tra, "psrf3", itau_w, zx_tmp_2d, &
575     iim*(jjm+1), ndex2d)
576    
577     CALL gr_fi_ecrit(1, klon, iim, jjm+1, ppsrf4, zx_tmp_2d)
578     CALL histwrite(nid_tra, "psrf4", itau_w, zx_tmp_2d, &
579     iim*(jjm+1), ndex2d)
580     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pplay, zx_tmp_3d)
581     CALL histwrite(nid_tra, "pplay", itau_w, zx_tmp_3d, &
582     iim*(jjm+1)*llm, ndex3d)
583    
584     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, t_seri, zx_tmp_3d)
585     CALL histwrite(nid_tra, "t", itau_w, zx_tmp_3d, &
586     iim*(jjm+1)*llm, ndex3d)
587     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pmfu, zx_tmp_3d)
588     CALL histwrite(nid_tra, "mfu", itau_w, zx_tmp_3d, &
589     iim*(jjm+1)*llm, ndex3d)
590     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pmfd, zx_tmp_3d)
591     CALL histwrite(nid_tra, "mfd", itau_w, zx_tmp_3d, &
592     iim*(jjm+1)*llm, ndex3d)
593     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pen_u, zx_tmp_3d)
594     CALL histwrite(nid_tra, "en_u", itau_w, zx_tmp_3d, &
595     iim*(jjm+1)*llm, ndex3d)
596     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pen_d, zx_tmp_3d)
597     CALL histwrite(nid_tra, "en_d", itau_w, zx_tmp_3d, &
598     iim*(jjm+1)*llm, ndex3d)
599     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pde_d, zx_tmp_3d)
600     CALL histwrite(nid_tra, "de_d", itau_w, zx_tmp_3d, &
601     iim*(jjm+1)*llm, ndex3d)
602     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, pde_u, zx_tmp_3d)
603     CALL histwrite(nid_tra, "de_u", itau_w, zx_tmp_3d, &
604     iim*(jjm+1)*llm, ndex3d)
605     CALL gr_fi_ecrit(llm, klon, iim, jjm+1, coefh, zx_tmp_3d)
606     CALL histwrite(nid_tra, "coefh", itau_w, zx_tmp_3d, &
607     iim*(jjm+1)*llm, ndex3d)
608    
609     ! abder
610    
611     if (ok_sync) then
612     call histsync(nid_tra)
613     endif
614    
615     end subroutine write_histrac
616    
617     END SUBROUTINE phytrac
618    
619     !*************************************************
620    
621     subroutine ini_histrac(nid_tra, pdtphys, presnivs, nqmax, lessivage)
622    
623     ! From phylmd/ini_histrac.h, version 1.10 2006/02/21 08:08:30
624    
625     use dimens_m, only: iim, jjm, llm
626     use ioipsl, only: ymds2ju, histbeg_totreg, histvert, histdef, histend
627     use temps, only: annee_ref, day_ref, itau_phy
628     use advtrac_m, only: niadv, tnom, ttext
629     use dimphy, only: klon
630     use clesphys, only: ecrit_tra
631     use grid_change, only: gr_phy_write
632     use phyetat0_m, only: rlon, rlat
633    
634     INTEGER, intent(out):: nid_tra
635     real, intent(in):: pdtphys ! pas d'integration pour la physique (s)
636     REAL, intent(in):: presnivs(:)
637    
638     integer, intent(in):: nqmax
639     ! (nombre de traceurs auxquels on applique la physique)
640    
641     logical, intent(in):: lessivage
642    
643     ! Variables local to the procedure:
644    
645     REAL zjulian
646     REAL zx_lat(iim, jjm+1)
647     INTEGER nhori, nvert
648     REAL zsto, zout
649     integer it, iq, iiq
650    
651     !---------------------------------------------------------
652    
653     CALL ymds2ju(annee_ref, month=1, day=day_ref, sec=0.0, julian=zjulian)
654     zx_lat(:, :) = gr_phy_write(rlat)
655     CALL histbeg_totreg("histrac", iim, rlon(2:iim+1), jjm+1, zx_lat(1, :), &
656     1, iim, 1, jjm+1, itau_phy, zjulian, pdtphys, nhori, nid_tra)
657     CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb", llm, &
658     presnivs, nvert)
659    
660     zsto = pdtphys
661     zout = pdtphys * REAL(ecrit_tra)
662    
663     CALL histdef(nid_tra, "phis", "Surface geop. height", "-", &
664     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
665     "once", zsto, zout)
666     CALL histdef(nid_tra, "aire", "Grid area", "-", &
667     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
668     "once", zsto, zout)
669     CALL histdef(nid_tra, "zmasse", "column density of air in cell", &
670     "kg m-2", iim, jjm + 1, nhori, llm, 1, llm, nvert, 32, "ave(X)", &
671     zsto, zout)
672    
673     DO it=1, nqmax
674     ! champ 2D
675     iq=it+2
676     iiq=niadv(iq)
677     CALL histdef(nid_tra, tnom(iq), ttext(iiq), "U/kga", iim, jjm+1, &
678     nhori, llm, 1, llm, nvert, 32, "ave(X)", zsto, zout)
679     if (lessivage) THEN
680     CALL histdef(nid_tra, "fl"//tnom(iq), "Flux "//ttext(iiq), &
681     "U/m2/s", iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
682     "ave(X)", zsto, zout)
683     endif
684    
685     !---Ajout Olivia
686     CALL histdef(nid_tra, "d_tr_th_"//tnom(iq), &
687     "tendance thermique"// ttext(iiq), "?", &
688     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
689     "ave(X)", zsto, zout)
690     CALL histdef(nid_tra, "d_tr_cv_"//tnom(iq), &
691     "tendance convection"// ttext(iiq), "?", &
692     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
693     "ave(X)", zsto, zout)
694     CALL histdef(nid_tra, "d_tr_cl_"//tnom(iq), &
695     "tendance couche limite"// ttext(iiq), "?", &
696     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
697     "ave(X)", zsto, zout)
698     !---fin Olivia
699    
700     ENDDO
701    
702     CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-", &
703     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
704     "inst(X)", zout, zout)
705    
706     CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-", &
707     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
708     "inst(X)", zout, zout)
709     CALL histdef(nid_tra, "psrf1", "nature sol", "-", &
710     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
711     "inst(X)", zout, zout)
712     CALL histdef(nid_tra, "psrf2", "nature sol", "-", &
713     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
714     "inst(X)", zout, zout)
715     CALL histdef(nid_tra, "psrf3", "nature sol", "-", &
716     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
717     "inst(X)", zout, zout)
718     CALL histdef(nid_tra, "psrf4", "nature sol", "-", &
719     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
720     "inst(X)", zout, zout)
721     CALL histdef(nid_tra, "ftsol1", "temper sol", "-", &
722     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
723     "inst(X)", zout, zout)
724     CALL histdef(nid_tra, "ftsol2", "temper sol", "-", &
725     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
726     "inst(X)", zout, zout)
727     CALL histdef(nid_tra, "ftsol3", "temper sol", "-", &
728     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
729     "inst(X)", zout, zout)
730     CALL histdef(nid_tra, "ftsol4", "temper sol", "-", &
731     iim, jjm+1, nhori, 1, 1, 1, -99, 32, &
732     "inst(X)", zout, zout)
733     CALL histdef(nid_tra, "pplay", "flux u mont", "-", &
734     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
735     "inst(X)", zout, zout)
736     CALL histdef(nid_tra, "t", "flux u mont", "-", &
737     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
738     "inst(X)", zout, zout)
739     CALL histdef(nid_tra, "mfu", "flux u mont", "-", &
740     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
741     "ave(X)", zsto, zout)
742     CALL histdef(nid_tra, "mfd", "flux u decen", "-", &
743     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
744     "ave(X)", zsto, zout)
745     CALL histdef(nid_tra, "en_u", "flux u mont", "-", &
746     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
747     "ave(X)", zsto, zout)
748     CALL histdef(nid_tra, "en_d", "flux u mont", "-", &
749     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
750     "ave(X)", zsto, zout)
751     CALL histdef(nid_tra, "de_d", "flux u mont", "-", &
752     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
753     "ave(X)", zsto, zout)
754     CALL histdef(nid_tra, "de_u", "flux u decen", "-", &
755     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
756     "ave(X)", zsto, zout)
757     CALL histdef(nid_tra, "coefh", "turbulent coef", "-", &
758     iim, jjm+1, nhori, llm, 1, llm, nvert, 32, &
759     "ave(X)", zsto, zout)
760    
761     CALL histend(nid_tra)
762    
763     end subroutine ini_histrac
764    
765     !*************************************************
766    
767     function radiornpb(tr_seri, pdtphys, tautr)
768    
769     ! From phylmd/radiornpb.F, v 1.2 2005/05/25 13:10:09
770    
771     ! Auteurs: AA + CG (LGGE/CNRS) Date 24-06-94
772     ! Objet: Decroissance radioactive d'un traceur dans l'atmosphere
773     !G 24 06 94 : Pour un traceur, le radon
774     !G 16 12 94 : Plus un 2eme traceur, le 210Pb. Le radon decroit en plomb.
775    
776     ! Le pas de temps "pdtphys" est supposé beaucoup plus petit que la
777     ! constante de temps de décroissance.
778    
779     use dimens_m, only: llm
780     use dimphy, only: klon, nbtr
781     use nrutil, only: assert
782    
783     IMPLICIT none
784    
785     REAL, intent(in):: tr_seri(:, :, :), pdtphys, tautr(:)
786     real radiornpb(klon, llm, 2)
787    
788     ! Variable local to the procedure:
789     INTEGER it
790    
791     !-----------------------------------------------
792    
793     call assert(shape(tr_seri) == (/klon, llm, nbtr/), "radiornpb tr_seri")
794     call assert(size(tautr) == nbtr, "radiornpb tautr")
795    
796     DO it = 1, 2
797     IF (tautr(it) > 0.) THEN
798     radiornpb(:, :, it) = - tr_seri(:, :, it) * pdtphys / tautr(it)
799     ELSE
800     radiornpb(:, :, it) = 0.
801     END IF
802     END DO
803    
804     !G161294 : Cas particulier radon 1 => plomb 2
805     radiornpb(:, :, 2) = radiornpb(:, :, 2) - radiornpb(:, :, 1)
806    
807     END function radiornpb
808    
809     !*************************************************
810    
811     SUBROUTINE minmaxqfi(zq, qmin, qmax, comment)
812    
813     ! From phylmd/minmaxqfi.F, version 1.1.1.1 2004/05/19 12:53:09
814    
815     use dimens_m, only: llm
816     use dimphy, only: klon
817     use nrutil, only: assert
818    
819     IMPLICIT none
820    
821     real, intent(in):: zq(:, :), qmin, qmax
822     CHARACTER(len=*), intent(in):: comment
823    
824     ! Variables local to the procedure:
825    
826     INTEGER jadrs(klon), jbad, k, i
827    
828     !---------------------------------
829    
830     call assert(shape(zq) == (/klon, llm/), "minmaxqfi")
831    
832     DO k = 1, llm
833     jbad = 0
834     DO i = 1, klon
835     IF (zq(i, k) > qmax .OR. zq(i, k) < qmin) THEN
836     jbad = jbad + 1
837     jadrs(jbad) = i
838     ENDIF
839     ENDDO
840     IF (jbad > 0) THEN
841     PRINT *, comment
842     DO i = 1, jbad
843     PRINT *, "zq(", jadrs(i), ", ", k, ") = ", zq(jadrs(i), k)
844     ENDDO
845     ENDIF
846     ENDDO
847    
848     end SUBROUTINE minmaxqfi
849    
850     end module phytrac_m

  ViewVC Help
Powered by ViewVC 1.1.21