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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 7 - (hide annotations)
Mon Mar 31 12:24:17 2008 UTC (16 years, 1 month ago) by guez
Original Path: trunk/libf/phylmd/phytrac.f90
File size: 29267 byte(s)
This revision is not in working order. Pending some moving of files.

Important changes. In the program "etat0_lim": ozone coefficients from
Mobidic are regridded in time instead of pressure ; consequences in
"etat0". In the program "gcm", ozone coefficients from Mobidic are
read once per day only for the current day and regridded in pressure ;
consequences in "o3_chem_m", "regr_pr_coefoz", "phytrac" and
"regr_pr_comb_coefoz_m".

NetCDF95 is a library and does not export NetCDF.

New variables "nag_gl_options", "nag_fcalls_options" and
"nag_cross_options" in "nag_tools.mk".

"check_coefoz.jnl" rewritten entirely for new version of
"coefoz_LMDZ.nc".

Target "obj_etat0_lim" moved from "GNUmakefile" to "nag_rules.mk".

Added some "intent" attributes in "calfis", "clmain", "clqh",
"cltrac", "cltracrn", "cvltr", "ini_undefSTD", "moy_undefSTD",
"nflxtr", "phystokenc", "phytrac", "readsulfate", "readsulfate_preind"
and "undefSTD".

In "dynetat0", "dynredem0" and "gcm", "phis" has rank 2 instead of
1. "phis" has assumed shape in "dynredem0".

Added module containing "dynredem0". Changed some calls with NetCDF
Fortran 77 interface to calls with NetCDF95 interface.

Replaced calls to "ssum" by calls to "sum" in "inigeom".

In "make.sh", new option "-c" to change compiler.

In "aaam_bud", argument "rjour" deleted.

In "physiq": renamed some variables; deleted variable "xjour".

In "phytrac": renamed some variables; new argument "lmt_pas".

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     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 guez 7 use regr_pr_comb_coefoz_m, only: regr_pr_comb_coefoz
41 guez 3 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 guez 7 integer, intent(in):: itap ! number of calls to "physiq"
56     integer, intent(in):: lmt_pas ! number of time steps of "physics" per day
57 guez 3 integer, intent(in):: julien !jour julien, 1 <= julien <= 360
58     integer itop_con(klon)
59     integer ibas_con(klon)
60     real, intent(in):: gmtime ! heure de la journée en fraction de jour
61 guez 7 real, intent(in):: pdtphys ! pas d'integration pour la physique (s)
62 guez 3 real, intent(in):: t_seri(klon, llm) ! temperature, in K
63    
64     real tr_seri(klon, llm, nbtr)
65     ! (mass fractions of tracers, excluding water, at mid-layers)
66    
67     real u(klon, llm)
68     real v(klon, llm)
69     real sh(klon, llm) ! humidite specifique
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     real pplay(klon, llm) ! pression pour le mileu de chaque couche (en Pa)
84     real pphi(klon, llm) ! geopotentiel
85     real pphis(klon)
86     REAL, intent(in):: presnivs(llm)
87 guez 7 logical, intent(in):: firstcal ! first call to "calfis"
88 guez 3 logical, intent(in):: lafin ! fin de la physique
89    
90     integer nsplit
91     REAL pmflxr(klon, llm+1), pmflxs(klon, llm+1) !--lessivage convection
92     REAL prfl(klon, llm+1), psfl(klon, llm+1) !--lessivage large-scale
93    
94     ! convection:
95    
96     REAL pmfu(klon, llm) ! flux de masse dans le panache montant
97     REAL pmfd(klon, llm) ! flux de masse dans le panache descendant
98     REAL pen_u(klon, llm) ! flux entraine dans le panache montant
99    
100     ! thermiques:
101    
102     real fm_therm(klon, llm+1), entr_therm(klon, llm)
103    
104     REAL pde_u(klon, llm) ! flux detraine dans le panache montant
105     REAL pen_d(klon, llm) ! flux entraine dans le panache descendant
106     REAL pde_d(klon, llm) ! flux detraine dans le panache descendant
107     ! KE
108     real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
109     REAL upwd(klon, llm) ! saturated updraft mass flux
110     REAL dnwd(klon, llm) ! saturated downdraft mass flux
111    
112     ! Couche limite:
113    
114     REAL coefh(klon, llm) ! coeff melange CL
115     REAL yu1(klon) ! vents au premier niveau
116     REAL yv1(klon) ! vents au premier niveau
117    
118     ! Lessivage:
119    
120     ! pour le ON-LINE
121    
122     REAL frac_impa(klon, llm) ! fraction d'aerosols impactes
123     REAL frac_nucl(klon, llm) ! fraction d'aerosols nuclees
124    
125     ! Arguments necessaires pour les sources et puits de traceur:
126    
127     real ftsol(klon, nbsrf) ! Temperature du sol (surf)(Kelvin)
128     real pctsrf(klon, nbsrf) ! Pourcentage de sol f(nature du sol)
129    
130     real pftsol1(klon), pftsol2(klon), pftsol3(klon), pftsol4(klon)
131     real ppsrf1(klon), ppsrf2(klon), ppsrf3(klon), ppsrf4(klon)
132    
133     ! VARIABLES LOCALES TRACEURS
134    
135     ! Sources et puits des traceurs:
136    
137     ! Pour l'instant seuls les cas du rn et du pb ont ete envisages.
138    
139     REAL source(klon) ! a voir lorsque le flux est prescrit
140     !
141     ! Pour la source de radon et son reservoir de sol
142    
143     REAL, save:: trs(klon, nbtr) ! Concentration de radon dans le sol
144    
145     REAL masktr(klon, nbtr) ! Masque reservoir de sol traceur
146     ! Masque de l'echange avec la surface
147     ! (1 = reservoir) ou (possible => 1 )
148     SAVE masktr
149     REAL fshtr(klon, nbtr) ! Flux surfacique dans le reservoir de sol
150     SAVE fshtr
151     REAL hsoltr(nbtr) ! Epaisseur equivalente du reservoir de sol
152     SAVE hsoltr
153     REAL tautr(nbtr) ! Constante de decroissance radioactive
154     SAVE tautr
155     REAL vdeptr(nbtr) ! Vitesse de depot sec dans la couche Brownienne
156     SAVE vdeptr
157     REAL scavtr(nbtr) ! Coefficient de lessivage
158     SAVE scavtr
159    
160     CHARACTER itn
161     INTEGER, save:: nid_tra
162    
163     ! nature du traceur
164    
165     logical aerosol(nbtr) ! Nature du traceur
166     ! ! aerosol(it) = true => aerosol
167     ! ! aerosol(it) = false => gaz
168     logical clsol(nbtr) ! couche limite sol calculée
169     logical radio(nbtr) ! décroisssance radioactive
170     save aerosol, clsol, radio
171    
172     ! convection tiedtke
173     INTEGER i, k, it
174     REAL delp(klon, llm)
175    
176     ! Variables liees a l'ecriture de la bande histoire physique
177    
178     ! Variables locales pour effectuer les appels en serie
179    
180     REAL d_tr(klon, llm), d_trs(klon) ! tendances de traceurs
181     REAL d_tr_cl(klon, llm, nbtr) ! tendance de traceurs couche limite
182     REAL d_tr_cv(klon, llm, nbtr) ! tendance de traceurs conv pour chq traceur
183     REAL d_tr_th(klon, llm, nbtr) ! la tendance des thermiques
184     REAL d_tr_dec(klon, llm, 2) ! la tendance de la decroissance
185     ! ! radioactive du rn - > pb
186     REAL d_tr_lessi_impa(klon, llm, nbtr) ! la tendance du lessivage
187     ! ! par impaction
188     REAL d_tr_lessi_nucl(klon, llm, nbtr) ! la tendance du lessivage
189     ! ! par nucleation
190     REAL flestottr(klon, llm, nbtr) ! flux de lessivage
191     ! ! dans chaque couche
192    
193     real zmasse(klon, llm)
194     ! (column-density of mass of air in a layer, in kg m-2)
195    
196     real ztra_th(klon, llm)
197    
198     character(len=20) modname
199     character(len=80) abort_message
200     integer isplit
201    
202     ! Controls:
203     logical:: couchelimite = .true.
204     logical:: convection = .true.
205     logical:: lessivage = .true.
206     logical, save:: inirnpb
207    
208     !--------------------------------------
209    
210     modname='phytrac'
211    
212 guez 7 if (firstcal) then
213 guez 3 print *, 'phytrac: pdtphys = ', pdtphys
214     PRINT *, 'Fréquence de sortie des traceurs : ecrit_tra = ', ecrit_tra
215     if (nbtr < nqmax) then
216     abort_message='See above'
217     call abort_gcm(modname, abort_message, 1)
218     endif
219     inirnpb=rnpb
220    
221     ! Initialisation des sorties :
222     call ini_histrac(nid_tra, pdtphys, presnivs, nqmax, lessivage)
223    
224     ! Initialisation de certaines variables pour le radon et le plomb
225     ! Initialisation du traceur dans le sol (couche limite radonique)
226     trs(:, :) = 0.
227    
228     open (unit=99, file='starttrac', status='old', err=999, &
229     form='formatted')
230     read(unit=99, fmt=*) (trs(i, 1), i=1, klon)
231     999 continue
232     close(unit=99)
233    
234     ! Initialisation de la fraction d'aerosols lessivee
235    
236     d_tr_lessi_impa(:, :, :) = 0.
237     d_tr_lessi_nucl(:, :, :) = 0.
238    
239     ! Initialisation de la nature des traceurs
240    
241     DO it = 1, nqmax
242     aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
243     radio(it) = .FALSE. ! par défaut pas de passage par "radiornpb"
244     clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
245     ENDDO
246     ENDIF
247    
248     ! Initialisation du traceur dans le sol (couche limite radonique)
249     if (inirnpb) THEN
250    
251     radio(1)= .true.
252     radio(2)= .true.
253     clsol(1)= .true.
254     clsol(2)= .true.
255     aerosol(2) = .TRUE. ! le Pb est un aerosol
256    
257     call initrrnpb(ftsol, pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, &
258     scavtr)
259     inirnpb=.false.
260     endif
261    
262     do i=1, klon
263     pftsol1(i) = ftsol(i, 1)
264     pftsol2(i) = ftsol(i, 2)
265     pftsol3(i) = ftsol(i, 3)
266     pftsol4(i) = ftsol(i, 4)
267    
268     ppsrf1(i) = pctsrf(i, 1)
269     ppsrf2(i) = pctsrf(i, 2)
270     ppsrf3(i) = pctsrf(i, 3)
271     ppsrf4(i) = pctsrf(i, 4)
272    
273     enddo
274    
275     ! Calcul de l'effet de la convection
276    
277     if (convection) then
278     DO it=1, nqmax
279     if (iflag_con.eq.2) then
280     ! tiedke
281     CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
282     pplay, paprs, tr_seri(1, 1, it), d_tr_cv(1, 1, it))
283     else if (iflag_con.eq.3) then
284     ! KE
285     call cvltr(pdtphys, da, phi, mp, paprs, pplay, &
286     tr_seri(1, 1, it), upwd, dnwd, d_tr_cv(1, 1, it))
287     endif
288    
289     DO k = 1, llm
290     DO i = 1, klon
291     tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cv(i, k, it)
292     ENDDO
293     ENDDO
294     WRITE(unit=itn, fmt='(i1)') it
295     CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, &
296     'convection, tracer index = ' // itn)
297     ENDDO
298     endif
299    
300     forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
301    
302     ! Calcul de l'effet des thermiques
303    
304     do it=1, nqmax
305     do k=1, llm
306     do i=1, klon
307     d_tr_th(i, k, it)=0.
308     tr_seri(i, k, it)=max(tr_seri(i, k, it), 0.)
309     tr_seri(i, k, it)=min(tr_seri(i, k, it), 1.e10)
310     enddo
311     enddo
312     enddo
313    
314     if (iflag_thermals > 0) then
315     nsplit=10
316     DO it=1, nqmax
317     do isplit=1, nsplit
318     ! Thermiques
319     call dqthermcell(klon, llm, pdtphys/nsplit &
320     , fm_therm, entr_therm, zmasse &
321     , tr_seri(1:klon, 1:llm, it), d_tr, ztra_th)
322    
323     do k=1, llm
324     do i=1, klon
325     d_tr(i, k)=pdtphys*d_tr(i, k)/nsplit
326     d_tr_th(i, k, it)=d_tr_th(i, k, it)+d_tr(i, k)
327     tr_seri(i, k, it)=max(tr_seri(i, k, it)+d_tr(i, k), 0.)
328     enddo
329     enddo
330     enddo
331     ENDDO
332     endif
333    
334     ! Calcul de l'effet de la couche limite
335    
336     if (couchelimite) then
337    
338     DO k = 1, llm
339     DO i = 1, klon
340     delp(i, k) = paprs(i, k)-paprs(i, k+1)
341     ENDDO
342     ENDDO
343    
344     ! MAF modif pour tenir compte du cas rnpb + traceur
345     DO it=1, nqmax
346     if (clsol(it)) then
347     ! couche limite avec quantite dans le sol calculee
348     CALL cltracrn(it, pdtphys, yu1, yv1, &
349     coefh, t_seri, ftsol, pctsrf, &
350     tr_seri(1, 1, it), trs(1, it), &
351     paprs, pplay, delp, &
352     masktr(1, it), fshtr(1, it), hsoltr(it), &
353     tautr(it), vdeptr(it), &
354     rlat, &
355     d_tr_cl(1, 1, it), d_trs)
356     DO k = 1, llm
357     DO i = 1, klon
358     tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
359     ENDDO
360     ENDDO
361    
362     ! Traceur ds sol
363    
364     DO i = 1, klon
365     trs(i, it) = trs(i, it) + d_trs(i)
366     END DO
367     else ! couche limite avec flux prescrit
368     !MAF provisoire source / traceur a creer
369     DO i=1, klon
370     source(i) = 0.0 ! pas de source, pour l'instant
371     ENDDO
372    
373     CALL cltrac(pdtphys, coefh, t_seri, &
374     tr_seri(1, 1, it), source, &
375     paprs, pplay, delp, &
376     d_tr_cl(1, 1, it))
377     DO k = 1, llm
378     DO i = 1, klon
379     tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
380     ENDDO
381     ENDDO
382     endif
383     ENDDO
384    
385     endif ! couche limite
386    
387     ! Calcul de l'effet du puits radioactif
388    
389     ! MAF il faudrait faire une modification pour passer dans radiornpb
390     ! si radio=true mais pour l'instant radiornpb propre au cas rnpb
391     if (rnpb) then
392     d_tr_dec(:, :, :) = radiornpb(tr_seri, pdtphys, tautr)
393     DO it=1, nqmax
394     if (radio(it)) then
395     tr_seri(:, :, it) = tr_seri(:, :, it) + d_tr_dec(:, :, it)
396     WRITE(unit=itn, fmt='(i1)') it
397     CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, 'puits rn it='//itn)
398     endif
399     ENDDO
400     endif ! rnpb decroissance radioactive
401    
402 guez 6 if (nqmax >= 3) then
403     ! Ozone as a tracer:
404 guez 7 if (mod(itap - 1, lmt_pas) == 0) then
405     ! Once per day, update the coefficients for ozone chemistry:
406     call regr_pr_comb_coefoz(julien)
407     end if
408 guez 6 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 guez 7 call write_histrac(lessivage, nqmax, itap, nid_tra)
464 guez 3
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 guez 7 subroutine write_histrac(lessivage, nqmax, itap, nid_tra)
478 guez 3
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 guez 7 integer, intent(in):: itap ! number of calls to "physiq"
494 guez 3 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 guez 7 itau_w = itau_phy + itap
508 guez 3
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