source: LMDZ.3.3/trunk/libf/phylmd/phytrac.F @ 344

Last change on this file since 344 was 344, checked in by lmdz, 22 years ago

Inclusion des modifs de D. Hauglustaine pour la version 1 de INCA
LF
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 39.4 KB
Line 
1c
2c $Header$
3c
4      SUBROUTINE phytrac (rnpb,
5     I                    nstep,
6     I                    julien,
7     I                    gmtime,
8     I                    debutphy,
9     I                    lafin,
10     I                    nqmax,
11     I                    nlon,
12     I                    nlev,
13     I                    pdtphys,
14     I                    u,
15     I                    v,
16     I                    t_seri,
17     I                    paprs,
18     I                    pplay,
19     I                    pmfu,
20     I                    pmfd,
21     I                    pen_u,
22     I                    pde_u,
23     I                    pen_d,
24     I                    pde_d,
25     I                    coefh,
26     I                    yu1,
27     I                    yv1,
28     I                    ftsol,
29     I                    pctsrf,
30     I                    xlat,
31     I                    frac_impa,
32     I                    frac_nucl,
33     I                    xlon,
34     I                    paire,
35     I                    pphis,
36     I                    pphi,
37     I                    albsol,
38     I                    sh,
39     I                    rh,
40     I                    cldfra,
41     I                    rneb,
42     I                    diafra,
43     I                    cldliq,
44     I                    itop_con,
45     I                    ibas_con,
46     I                    pmflxr,
47     I                    pmflxs,
48     I                    prfl,
49     I                    psfl,
50     I                    flxmass_w,
51     O                    tr_seri)
52      USE ioipsl
53
54#ifdef INCA
55      USE sflx
56      USE chem_tracnm
57      USE species_names
58      USE chem_mods
59      USE pht_tables, ONLY : jrates
60      USE transport_controls, ONLY : conv_flg, pbl_flg
61      USE airplane_src, ONLY : ptrop
62      USE lightning, ONLY : prod_light
63#ifdef INCA_AER
64      USE AEROSOL_DIAG
65#endif
66#endif
67
68      IMPLICIT none
69c======================================================================
70c Auteur(s) FH
71c Objet: Moniteur general des tendances traceurs
72c
73cAA Remarques en vrac:
74cAA--------------------
75cAA 1/ le call phytrac se fait avec nqmax-2 donc nous avons bien
76cAA les vrais traceurs (nbtr) dans phytrac (pas la vapeur ni eau liquide)
77cAA 2/ Le choix du radon et du pb se fait juste avec un data
78cAA    (peu propre). Peut-etre pourrait-on prevoir dans l'avenir
79cAA    une variable "type de traceur"
80c======================================================================
81#include "YOMCST.h"
82#include "dimensions.h"
83#include "paramet.h"
84#include "dimphy.h"
85#include "indicesol.h"
86#include "comvert.h"
87#include "temps.h"
88#include "control.h"
89c======================================================================
90
91c Arguments:
92c
93c   EN ENTREE:
94c   ==========
95c
96c   divers:
97c   -------
98c
99      integer nlon   ! nombre de points horizontaux
100      integer nlev   ! nombre de couches verticales
101      integer nqmax  ! nombre de traceurs auxquels on applique la physique
102      integer nstep  ! appel physique
103      integer julien !jour julien
104      integer itop_con(nlon)
105      integer ibas_con(nlon)
106      real gmtime
107      real pdtphys  ! pas d'integration pour la physique (seconde)
108      real t_seri(nlon,nlev) ! temperature
109      real u(klon,klev)
110      real v(klon,klev)
111      real sh(nlon,nlev)     ! humidite specifique
112      real rh(nlon,nlev)     ! humidite relative
113      real cldliq(nlon,nlev) ! eau liquide nuageuse
114      real cldfra(nlon,nlev) ! fraction nuageuse (tous les nuages)
115      real diafra(nlon,nlev) ! fraction nuageuse (convection ou stratus artificiels)
116      real rneb(nlon,nlev)   ! fraction nuageuse (grande echelle)
117      real albsol(nlon)      ! albedo surface
118      real tr_seri(nlon,nlev,nbtr) ! traceur 
119      real paprs(nlon,nlev+1)  ! pression pour chaque inter-couche (en Pa)
120      real ps(nlon)  ! pression surface
121      real pplay(nlon,nlev)  ! pression pour le mileu de chaque couche (en Pa)
122      real pphi(nlon,klev) ! geopotentiel
123      real znivsig(klev) ! niveaux sigma
124      real paire(klon)
125      real pphis(klon)
126      REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)   !--lessivage convection
127      REAL prfl(klon,klev+1),   psfl(klon,klev+1)     !--lessivage large-scale
128      REAL flxmass_w(klon,klev)
129      logical debutphy       ! le flag de l'initialisation de la physique
130      logical lafin          ! le flag de la fin de la physique
131
132      integer ll
133c
134cAA Rem : nbtr : nombre de vrais traceurs est defini dans dimphy.h
135c
136c   convection:
137c   -----------
138c
139      REAL pmfu(nlon,nlev)  ! flux de masse dans le panache montant
140      REAL pmfd(nlon,nlev)  ! flux de masse dans le panache descendant
141      REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
142      REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
143      REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
144      REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
145c
146c   Couche limite:
147c   --------------
148c
149      REAL coefh(nlon,nlev) ! coeff melange CL
150      REAL yu1(nlon)        ! vents au premier niveau
151      REAL yv1(nlon)        ! vents au premier niveau
152      REAL xlat(nlon)       ! latitudes pour chaque point
153      REAL xlon(nlon)       ! longitudes pour chaque point
154
155c
156c   Lessivage:
157c   ----------
158c
159c pour le ON-LINE
160c
161      REAL frac_impa(nlon,nlev)  ! fraction d'aerosols impactes
162      REAL frac_nucl(nlon,nlev)  ! fraction d'aerosols nuclees
163c
164cAA
165cAA Arguments necessaires pour les sources et puits de traceur:
166cAA ----------------
167cAA
168      real ftsol(nlon,nbsrf)  ! Temperature du sol (surf)(Kelvin)
169      real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol)
170c abder
171      real pftsol1(nlon),pftsol2(nlon),pftsol3(nlon),pftsol4(nlon)
172      real ppsrf1(nlon),ppsrf2(nlon),ppsrf3(nlon),ppsrf4(nlon)
173c fin
174cAA ----------------------------
175cAA  VARIABLES LOCALES TRACEURS
176cAA ----------------------------
177cAA
178cAA Sources et puits des traceurs:
179cAA ------------------------------
180cAA
181cAA Pour l'instant seuls les cas du rn et du pb ont ete envisages.
182
183      REAL source(klon)       ! a voir lorsque le flux est prescrit
184      REAL topuits(klev,nbtr)  ! a voir
185cAA
186cAA Pour la source de radon et son reservoir de sol
187cAA ................................................
188 
189      REAL trs(klon,nbtr)    ! Conc. radon ds le sol
190      SAVE trs
191
192      REAL masktr(klon,nbtr) ! Masque reservoir de sol traceur
193c                            Masque de l'echange avec la surface
194c                           (1 = reservoir) ou (possible => 1 )
195      SAVE masktr
196      REAL fshtr(klon,nbtr)  ! Flux surfacique dans le reservoir de sol
197      SAVE fshtr
198      REAL hsoltr(nbtr)      ! Epaisseur equivalente du reservoir de sol
199      SAVE hsoltr
200      REAL tautr(nbtr)       ! Constante de decroissance radioactive
201      SAVE tautr
202      REAL vdeptr(nbtr)      ! Vitesse de depot sec dans la couche Brownienne
203      SAVE vdeptr
204      REAL scavtr(nbtr)      ! Coefficient de lessivage
205      SAVE scavtr
206cAA
207      CHARACTER*2 itn
208C maf ioipsl
209      CHARACTER*2 str2
210      INTEGER nhori, nvert, nverta, nvertb, nverts
211      REAL zsto, zout, zjulian
212      INTEGER nid_tra
213      SAVE nid_tra
214c     REAL x(klon,klev,nbtr+2) ! traceurs
215      INTEGER ndex(1)
216      REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
217      REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
218      INTEGER itra
219      SAVE itra                   ! compteur pour la physique
220C
221C Variables liees a l'ecriture de la bande histoire : phytrac.nc
222c
223      INTEGER ecrit_tra
224      SAVE ecrit_tra   
225      logical ok_sync
226      parameter (ok_sync = .true.)
227C
228C nature du traceur
229c
230      logical aerosol(nbtr)  ! Nature du traceur
231c                            ! aerosol(it) = true  => aerosol
232c                            ! aerosol(it) = false => gaz
233c                            ! nat_trac(it) = 1. aerosolc
234      logical clsol(nbtr)    ! clsol(it) = true => CL sol calculee
235      logical radio(nbtr)    ! radio(it)=true => decroisssance radioactive
236      save aerosol,clsol,radio
237C
238c======================================================================
239c
240c Declaration des procedures appelees
241c
242c--modif convection tiedtke
243      INTEGER i, k, it,itap
244        save itap
245      REAL delp(klon,klev)
246c--end modif
247c
248c Variables liees a l'ecriture de la bande histoire physique
249c
250c Variables locales pour effectuer les appels en serie
251c----------------------------------------------------
252c
253      REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
254      REAL d_tr_cl(klon,klev) ! tendance de traceurs  couche limite
255      REAL d_tr_cv(klon,klev) ! tendance de traceurs  convection
256      REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
257c                                   ! radioactive du rn - > pb
258      REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
259c                                          ! par impaction
260      REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
261c                                          ! par nucleation
262      REAL fluxrn(klon,klev)
263      REAL fluxpb(klon,klev)
264      REAL pbimpa(klon,klev)
265      REAL pbnucl(klon,klev)
266      REAL rn(klon,klev)
267      REAL pb(klon,klev)
268      REAL flestottr(klon,klev,nbtr) ! flux de lessivage
269c                                    ! dans chaque couche
270
271C
272      character*20 modname
273      character*80 abort_message
274c
275c   Controles
276c-------------
277      logical first,couchelimite,convection,lessivage,sorties,
278     s        rnpb,inirnpb
279      save first,couchelimite,convection,lessivage,sorties,
280     s     inirnpb
281      data first,couchelimite,convection,lessivage,sorties
282     s     /.true.,.true.,.true.,.true.,.true./
283
284#ifdef INCA
285      INTEGER           :: ncsec
286      INTEGER           :: prt_flag_ts(nbtr)=(/1,1,1
287#ifdef INCA_CH4
288     .                                              ,0,0,1,1,1,1,1,
289     .                                         0,1,0,0,0,0,0,1,0,0,
290     .                                         0,1,1,1,1,0,1,1,1,0,
291     .                                         1,1,1,1,1,1,1,1,1,1,
292     .                                         1,0,0
293#endif
294#ifdef INCA_AER
295     .                                              ,1,1,1,1,1,1
296#endif
297     .                                         /)
298      REAL, PARAMETER   :: dry_mass = 28.966
299      REAL, POINTER     :: hbuf(:)
300      REAL, ALLOCATABLE :: obuf(:)
301      REAL              :: calday
302      REAL              :: pdel(klon,klev)
303      REAL              :: dummy(klon,klev) = 0.
304#endif
305
306c
307c======================================================================
308c         ecrit_tra = NINT(86400./pdtphys *1.0)  ! tous les jours
309         modname='phytrac'
310         ecrit_tra = NINT(86400./pdtphys *ecritphy)   
311!DH      print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
312         ps(:) = paprs(:,1)
313
314c        print*,'DANS PHYTRAC debutphy=',debutphy
315
316         if (debutphy) then
317
318          print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
319          ecrit_tra = NINT(86400./pdtphys/2.) ! tous les 12H
320c         ecrit_tra = NINT(86400./pdtphys) ! tous les 24H
321
322         if(nbtr.lt.nqmax) then
323           print*,'NQMAX=',nqmax
324           print*,'NBTR=',nbtr
325           abort_message='See above'
326           call abort_gcm(modname,abort_message,1)
327         endif
328
329         inirnpb=rnpb
330!DH      PRINT*, 'La frequence de sortie traceurs est  ', ecrit_tra
331
332C         
333         CALL ymds2ju(anne_ini, 1, 1, 0.0, zjulian)
334         itra = NINT(FLOAT(day_ini)*86400./pdtphys)
335         itap = NINT(FLOAT(day_ini)*86400./pdtphys)
336!        itra=0
337!        itap=0
338!        zjulian = zjulian + day_ini
339 
340         CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlon,zx_lon)
341         DO i = 1, iim
342            zx_lon(i,1) = xlon(i+1)
343            zx_lon(i,jjm+1) = xlon(i+1)
344         ENDDO
345         DO ll=1,klev
346            znivsig(ll)=float(ll)
347         ENDDO
348         CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlat,zx_lat)
349         CALL histbeg("histrac", iim,zx_lon, jjm+1,zx_lat,
350     .                 1,iim,1,jjm+1, itra, zjulian, pdtphys,
351     .                 nhori, nid_tra)
352
353         call histvert(nid_tra, "presnivs", "presnivs", "mb",
354     .                 klev, presnivs, nvert)
355#ifdef INCA
356!        call histvert(nid_tra, "ap", "Hybrid A parameter", "-",
357!    .                 klev+1, ap, nverta)
358!        call histvert(nid_tra, "bp", "Hybrid B parameter", "-",
359!    .                 klev+1, bp, nvertb)
360#endif
361
362         zsto = pdtphys
363         zout = pdtphys * FLOAT(ecrit_tra)
364c
365         CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
366     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
367     .                "once",  zsto,zout)
368c
369         CALL histdef(nid_tra, "aire", "Grid area", "-",
370     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
371     .                "once",  zsto,zout)
372
373#ifdef INCA
374         CALL histdef(nid_tra, "ps", "Surface pressure", "Pa",
375     .                iim,jjm+1,nhori, 1,1,1,-99, 32,
376     .                "ave(X)", zsto,zout)
377
378         CALL histdef(nid_tra, "ptrop", "Tropopause pressure", "Pa",
379     .                iim,jjm+1,nhori, 1,1,1,-99, 32,
380     .                "ave(X)", zsto,zout)
381
382         CALL histdef(nid_tra, "temp", "Air temperature", "K",
383     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
384     .                "ave(X)", zsto,zout)
385
386         CALL histdef(nid_tra, "u", "zonal wind component", "m/s",
387     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
388     .                "ave(X)", zsto,zout)
389
390         CALL histdef(nid_tra, "v", "zonal wind component", "m/s",
391     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
392     .                "ave(X)", zsto,zout)
393
394         CALL histdef(nid_tra, "h2o", "Specific Humidity", "MMR",
395     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
396     .                "ave(X)", zsto,zout)
397
398         CALL histdef(nid_tra, "pmid", "Pressure", "Pa",
399     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
400     .                "ave(X)", zsto,zout)
401
402         CALL histdef(nid_tra, "pdel", "Delta Pressure", "Pa",
403     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
404     .                "ave(X)", zsto,zout)
405#ifdef INCA_CH4
406#ifdef INCAINFO
407         DO it=1, phtcnt
408         WRITE(str2,'(i2.2)') it
409         CALL histdef(nid_tra, "j"//str2,"j"//str2, "CM-3 S-1",
410     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
411     .                "ave(X)", zsto,zout)
412         ENDDO
413         DO it=1, hetcnt
414         WRITE(str2,'(i2.2)') it
415         CALL histdef(nid_tra, "w"//str2,"w"//str2, "S-1",
416     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
417     .                "ave(X)", zsto,zout)
418         ENDDO
419
420         DO it=1, extcnt
421         WRITE(str2,'(i2.2)') it
422        CALL histdef(nid_tra, "ext"//str2,"ext"//str2, "CM-3 S-1",
423     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
424     .                "ave(X)", zsto,zout)
425         ENDDO
426
427         DO it=1, nfs
428         WRITE(str2,'(i2.2)') it
429         CALL histdef(nid_tra, "INV"//str2, "INV"//str2, "CM-3",
430     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
431     .                "ave(X)", zsto,zout)
432         ENDDO
433
434#else
435         CALL histdef(nid_tra, "jO3","jO3", "CM-3 S-1",
436     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
437     .                "ave(X)", zsto,zout)
438         CALL histdef(nid_tra, "jNO2","jNO2", "CM-3 S-1",
439     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
440     .                "ave(X)", zsto,zout)
441         CALL histdef(nid_tra, "jH2O2","jH2O2", "CM-3 S-1",
442     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
443     .                "ave(X)", zsto,zout)
444         CALL histdef(nid_tra, "wHNO3","wHNO3", "S-1",
445     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
446     .                "ave(X)", zsto,zout)
447         CALL histdef(nid_tra, "kN2O5", "kN2O5","CM-3 S-1",
448     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
449     .                "ave(X)", zsto,zout)
450         CALL histdef(nid_tra, "LghtNO","LghtNO", "CM-3 S-1",
451     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
452     .                "ave(X)", zsto,zout)
453#endif
454
455         DO it=1, grpcnt
456         CALL histdef(nid_tra, grpsym(it), grpsym(it), "VMR",
457     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
458     .                "ave(X)", zsto,zout)
459         ENDDO
460#endif
461#endif
462
463#ifdef INCA_AER
464C   3d FIELDS
465
466        CALL histdef(nid_tra, "scavcoef_st","scavcoef_st", "S-1",
467     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
468     .                "ave(X)", zsto,zout)
469        CALL histdef(nid_tra, "scavcoef_cv","scavcoef_cv", "S-1",
470     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
471     .                "ave(X)", zsto,zout)
472#endif
473 
474         DO it=1,nqmax
475         IF (it.LE.99) THEN
476         WRITE(str2,'(i2.2)') it
477C champ 2D
478
479#ifdef INCA
480         IF ( prt_flag_ts(it) == 0 ) CYCLE
481
482         CALL histdef(nid_tra, "Emi_"//solsym(it), "Emi_"//solsym(it),
483     .           "kg/m2/s", iim,jjm+1,nhori, 1,1,1, -99, 32,
484     .           "ave(X)", zsto,zout)
485         CALL histdef(nid_tra, "Dep_"//solsym(it), "Dep_"//solsym(it),
486     .           "cm/s", iim,jjm+1,nhori, 1,1,1, -99, 32,
487     .           "ave(X)", zsto,zout)
488#ifdef INCA_AER
489        IF (solsym(it) .eq. 'CIDUSTM'.or.solsym(it) .eq. 'CIN'
490     .  .or.solsym(it) .eq. 'CSSSM'  .or.solsym(it) .eq. 'CSN'
491     .  .or.solsym(it) .eq. 'ASSSM'  .or.solsym(it) .eq. 'ASN' ) THEN
492         CALL histdef(nid_tra, "Sed_"//solsym(it), "Sed_"//solsym(it),
493     .           "kg/m2", iim,jjm+1,nhori, 1,1,1, -99, 32,
494     .           "ave(X)", zsto,zout)
495        endif
496        IF (solsym(it) .eq. 'CIDUSTM') THEN
497         CALL histdef(nid_tra, "OD_"//solsym(it), "OD_"//solsym(it),
498     .           "opt. depth", iim,jjm+1,nhori, 1,1,1, -99, 32,
499     .           "ave(X)", zsto,zout)
500        endif
501
502#endif
503         CALL histdef(nid_tra, solsym(it), solsym(it), "VMR",
504     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
505     .                "ave(X)", zsto,zout)
506#else
507         CALL histdef(nid_tra, "tr"//str2, "Tracer No."//str2, "U/kga",
508     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
509     .                "ave(X)", zsto,zout)
510         IF (lessivage) THEN
511         CALL histdef(nid_tra, "fl"//str2, "Tracer Flux No."//str2,
512     .                "U/m2/s",iim,jjm+1,nhori, klev,1,klev,nvert, 32,
513     .                "ave(X)", zsto,zout)
514         ENDIF
515#endif
516         ELSE
517         PRINT*, "Trop de traceurs"
518         CALL abort
519         ENDIF
520         ENDDO
521
522#ifdef INCA_CH4
523         CALL histdef(nid_tra, "O3_column", "O3_column",
524     .           "DU", iim,jjm+1,nhori, 1,1,1, -99, 32,
525     .           "ave(X)", zsto,zout)
526         CALL histdef(nid_tra, "CO_column", "CO_column",
527     .           "10^18 CM-2", iim,jjm+1,nhori, 1,1,1, -99, 32,
528     .           "ave(X)", zsto,zout)
529         CALL histdef(nid_tra, "CH4_column", "CH4_column",
530     .           "10^18 CM-2", iim,jjm+1,nhori, 1,1,1, -99, 32,
531     .           "ave(X)", zsto,zout)
532         CALL histdef(nid_tra, "NO2_column", "NO2_column",
533     .           "10^15 CM-2", iim,jjm+1,nhori, 1,1,1, -99, 32,
534     .           "ave(X)", zsto,zout)
535         CALL histdef(nid_tra, "O3_ste", "O3_ste",
536     .           "CM-2 S-1", iim,jjm+1,nhori, 1,1,1, -99, 32,
537     .           "ave(X)", zsto,zout)
538         CALL histdef(nid_tra, "O3_prod", "O3_prod", "CM-3 S-1",
539     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
540     .                "ave(X)", zsto,zout)
541         CALL histdef(nid_tra, "O3_loss", "O3_loss", "CM-3 S-1",
542     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
543     .                "ave(X)", zsto,zout)
544
545!        Special variables for daytime averaging
546!        CALL histdef(nid_tra, "day_cnt", "day_cnt", "-",
547!    .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
548!    .                "t_sum(X)", zsto,zout)
549!        CALL histdef(nid_tra, "NO_day", "NO_day", "VMR",
550!    .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
551!    .                "t_sum(X)", zsto,zout)
552
553#endif
554
555         CALL histend(nid_tra)
556         ndex(1) = 0
557c
558         i = NINT(zout/zsto)
559         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
560         CALL histwrite(nid_tra,"phis",i,zx_tmp_2d,iim*(jjm+1),ndex)
561C
562         i = NINT(zout/zsto)
563         CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
564         CALL histwrite(nid_tra,"aire",i,zx_tmp_2d,iim*(jjm+1),ndex)
565
566c======================================================================
567c   Initialisation de certaines variables pour le Rn et le Pb
568c======================================================================
569
570c Initialisation du traceur dans le sol (couche limite radonique)
571c
572c        print*,'valeur de debut dans phytrac :',debutphy
573         do it=1,nqmax
574            do i=1,klon
575               trs(i,it) = 0.
576            enddo
577         END DO
578
579         open (99,file='starttrac',status='old',
580     .         err=999,form='formatted')
581         read(99,*) (trs(i,1),i=1,klon)
582999      close(99)
583!DH      print*, 'apres starttrac'
584
585c Initialisation de la fraction d'aerosols lessivee
586c
587         DO it=1,nqmax
588           DO k = 1, nlev
589              DO i = 1, klon
590                 d_tr_lessi_impa(i,k,it) =  0.0
591                 d_tr_lessi_nucl(i,k,it) =  0.0
592             ENDDO
593           ENDDO
594         ENDDO
595c
596c Initialisation de la nature des traceurs
597c
598         DO it = 1, nqmax
599            aerosol(it) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
600            radio(it) = .FALSE.    ! Par defaut pas de passage par radiornpb
601            clsol(it) = .FALSE.  ! Par defaut couche limite avec flux prescrit
602         ENDDO
603c
604      ENDIF  ! fin debutphy
605c Initialisation du traceur dans le sol (couche limite radonique)
606      if(inirnpb) THEN
607c
608         radio(1)= .true.
609         radio(2)= .true.
610         clsol(1)= .true.
611         clsol(2)= .true.
612         aerosol(2) = .TRUE. ! le Pb est un aerosol
613c
614         call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
615     .                   ,vdeptr,scavtr)
616         inirnpb=.false.
617      endif
618
619#ifdef INCA
620!======================================================================
621!     Chimie
622!======================================================================
623
624        calday = FLOAT(julien) + gmtime
625        ncsec  = NINT (86400.*gmtime)
626
627        DO k = 1, nlev
628        pdel(:,k) = paprs(:,k) - paprs (:,k+1)
629        END DO
630
631#ifdef INCAINFO
632        PRINT *, 'CHEMMAIN @ ', calday, ' ... '
633        DO it = 1, nbtr
634        PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
635     $                       MAXVAL(tr_seri(:,:,it))
636      END DO
637#endif
638
639
640#ifdef INCA_AER
641        CALL aerosolmain (tr_seri,
642     $                 pdtphys,
643     $                 pplay,
644     $                 prfl,
645     $                 pmflxr,
646     $                 psfl,
647     $                 pmflxs,
648     $                 pmfu,
649     $                 itop_con,
650     $                 ibas_con,
651     $                 pphi,
652     $                 paire,
653     $                 nstep)
654#endif
655
656        CALL chemmain (tr_seri,
657     $                 nas,
658     $                 nstep,
659     $                 calday,
660     $                 julien,
661     $                 ncsec,
662     $                 1,
663     $                 pdtphys,
664     $                 paprs(1,1),
665     $                 pplay,
666     $                 pdel,
667     $                 pctsrf(1,3),
668     $                 ftsol,
669     $                 albsol,
670     $                 pphi,
671     $                 pphis,
672     $                 cldfra,
673     $                 rneb,
674     $                 diafra,
675     $                 itop_con,
676     $                 cldliq,
677     $                 prfl,
678     $                 pmflxr,
679     $                 psfl,
680     $                 pmflxs,
681     $                 pmfu,
682     $                 flxmass_w,
683     $                 t_seri,
684     $                 sh,
685     $                 rh,
686     $                 .false.,
687     $                 hbuf,
688     $                 obuf,
689     $                 iip1,
690     $                 jjp1)
691#ifdef INCAINFO
692      PRINT *, 'OK.'
693      DO it = 1, nbtr
694      PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
695     $                     MAXVAL(tr_seri(:,:,it))
696      END DO
697#endif
698#endif
699C
700
701c======================================================================
702c   Calcul de l'effet de la convection
703c======================================================================
704!DH   print*,'Avant convection'
705      do it=1,nqmax
706         WRITE(itn,'(i1)') it
707c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
708      enddo
709
710      if (convection) then
711
712!DH   print*,'Pas de temps dans phytrac : ',pdtphys
713      DO it=1, nqmax
714#ifdef INCA
715      IF ( conv_flg(it) == 0 ) CYCLE
716#endif
717      CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
718     .            pplay, paprs, tr_seri(1,1,it), d_tr_cv)
719      DO k = 1, nlev
720      DO i = 1, klon
721         tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k)
722      ENDDO
723      ENDDO
724c     WRITE(itn,'(i2)') it
725#ifdef INCA
726      CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//solsym(it))
727#else
728      CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
729#endif
730      ENDDO
731c     print*,'apres nflxtr'
732
733      endif ! convection
734c        print*,'Apres convection'
735c      do it=1,nqmax
736c         WRITE(itn,'(i1)') it
737c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
738c      enddo
739
740c======================================================================
741c   Calcul de l'effet de la couche limite
742c======================================================================
743c       print *,'Avant couchelimite'
744c      do it=1,nqmax
745c         WRITE(itn,'(i1)') it
746c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL  '//itn)
747c      enddo
748
749      if (couchelimite) then
750
751      DO k = 1, nlev
752      DO i = 1, klon
753         delp(i,k) = paprs(i,k)-paprs(i,k+1)
754      ENDDO
755      ENDDO
756
757C maf modif pour tenir compte du cas rnpb + traceur
758      DO it=1, nqmax
759#ifdef INCA
760      IF ( pbl_flg(it) == 0 ) CYCLE
761#endif
762C     print *,'it',it,clsol(it)
763      if (clsol(it)) then  ! couche limite avec quantite dans le sol calculee
764          CALL cltracrn(it, pdtphys, yu1, yv1,
765     e                    coefh,t_seri,ftsol,pctsrf,
766     e                    tr_seri(1,1,it),trs(1,it),
767     e                    paprs, pplay, delp,
768     e                    masktr(1,it),fshtr(1,it),hsoltr(it),
769     e                    tautr(it),vdeptr(it),
770     e                    xlat,
771     s                    d_tr_cl,d_trs)
772          DO k = 1, nlev
773            DO i = 1, klon
774              tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k)
775            ENDDO
776          ENDDO
777c
778c Traceur ds sol
779c
780          DO i = 1, klon
781            trs(i,it) = trs(i,it) + d_trs(i)
782          END DO
783C
784C maf provisoire suppression des prints
785C         WRITE(itn,'(i1)') it
786C         CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
787      else ! couche limite avec flux prescrit
788
789#ifdef INCA
790        DO k =  1, klon
791          source(k) = eflux(k,it)-dflux(k,it)
792        END DO
793#else
794Cmaf provisoire source / traceur a creer
795        DO i=1, klon
796          source(i) = 0.0 ! pas de source, pour l'instant
797        ENDDO
798#endif
799          CALL cltrac(pdtphys, coefh,t_seri,
800     s               tr_seri(1,1,it), source,
801     e               paprs, pplay, delp,
802     s               d_tr )
803            DO k = 1, nlev
804               DO i = 1, klon
805                  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr(i,k)
806               ENDDO
807            ENDDO
808Cmaf provisoire suppression des prints
809Cmaf          WRITE(itn,'(i1)') it
810cmaf          CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
811      endif
812      ENDDO
813c
814      endif ! couche limite
815
816c      print*,'Apres couchelimite'
817c      do it=1,nqmax
818c         WRITE(itn,'(i1)') it
819c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL  '//itn)
820c      enddo
821
822c======================================================================
823c   Calcul de l'effet du puits radioactif
824c======================================================================
825
826C MAF il faudrait faire une modification pour passer dans radiornpb
827c si radio=true mais pour l'instant radiornpb propre au cas rnpb
828      if(rnpb) then
829        print *, 'decroissance radiactive activee'
830        call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
831C
832        DO it=1,nqmax
833            if(radio(it)) then
834            DO k = 1, nlev
835               DO i = 1, klon
836                  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
837               ENDDO
838            ENDDO
839            WRITE(itn,'(i1)') it
840            CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
841            endif
842        ENDDO
843c
844      endif ! rnpb decroissance  radioactive
845
846c======================================================================
847c   Calcul de l'effet de la precipitation
848c======================================================================
849
850!DH   print*,'LESSIVAGE =',lessivage
851      IF (lessivage) THEN
852
853c     print*,'avant lessivage'
854
855       DO it = 1, nqmax
856           DO k = 1, nlev
857              DO i = 1, klon
858               d_tr_lessi_nucl(i,k,it) = 0.0
859               d_tr_lessi_impa(i,k,it) = 0.0
860               flestottr(i,k,it) = 0.0
861              ENDDO
862           ENDDO
863       ENDDO
864c
865c tendance des aerosols nuclees et impactes
866c
867       DO it = 1, nqmax
868         IF (aerosol(it)) THEN
869           DO k = 1, nlev
870              DO i = 1, klon
871               d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
872     s                  ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
873               d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
874     s                  ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
875              ENDDO
876           ENDDO
877         ENDIF
878       ENDDO
879c
880c Mises a jour des traceurs + calcul des flux de lessivage
881c Mise a jour due a l'impaction et a la nucleation
882c
883c      call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
884c      call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
885c      call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
886       DO it = 1, nqmax
887c         print*,'IT=',it,aerosol(it)
888         IF (aerosol(it)) THEN
889c           print*,'IT=',it,' On lessive'
890           DO k = 1, nlev
891              DO i = 1, klon
892               tr_seri(i,k,it)=tr_seri(i,k,it)
893     s         *frac_impa(i,k)*frac_nucl(i,k)
894              ENDDO
895           ENDDO
896         ENDIF
897       ENDDO
898c      call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
899c
900c Flux lessivage total
901c
902      DO it = 1, nqmax
903           DO k = 1, nlev
904            DO i = 1, klon
905               flestottr(i,k,it) = flestottr(i,k,it) -
906     s                   ( d_tr_lessi_nucl(i,k,it)   +
907     s                     d_tr_lessi_impa(i,k,it) ) *
908     s                   ( paprs(i,k)-paprs(i,k+1) ) /
909     s                   (RG * pdtphys)
910            ENDDO
911           ENDDO
912c
913Cmaf suppression provisoire
914Cmaf        WRITE(itn,'(i1)') it
915Cmaf    CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
916      ENDDO
917c
918c     print*,'apres lessivage'
919      ENDIF
920Cc
921      DO k = 1, nlev
922         DO i = 1, klon
923            fluxrn(i,k) = flestottr(i,k,1)
924            fluxpb(i,k) = flestottr(i,k,2)
925            rn(i,k) = tr_seri(i,k,1)
926            pb(i,k) = tr_seri(i,k,2)
927            pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
928            pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
929         ENDDO
930      ENDDO
931
932
933
934      itra=itra+1
935      ndex(1) = 0
936
937      i = NINT(zout/zsto)
938      CALL gr_fi_ecrit(1, klon,iim,jjm+1, paire,zx_tmp_2d)
939      CALL histwrite(nid_tra,"aire",itra,zx_tmp_2d,
940     .     iim*(jjm+1),ndex)
941
942#ifdef INCA
943      CALL gr_fi_ecrit(1, klon,iim,jjm+1, ps,zx_tmp_2d)
944      CALL histwrite(nid_tra,"ps",itra,zx_tmp_2d,
945     .     iim*(jjm+1),ndex)
946
947      CALL gr_fi_ecrit(1, klon,iim,jjm+1, ptrop,zx_tmp_2d)
948      CALL histwrite(nid_tra,"ptrop",itra,zx_tmp_2d,
949     .     iim*(jjm+1),ndex)
950
951      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,t_seri, zx_tmp_3d)
952      CALL histwrite(nid_tra,"temp",itra,zx_tmp_3d,
953     .                                   iim*(jjm+1)*klev,ndex)
954
955      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,u, zx_tmp_3d)
956      CALL histwrite(nid_tra,"u",itra,zx_tmp_3d,
957     .                                   iim*(jjm+1)*klev,ndex)
958
959      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,v, zx_tmp_3d)
960      CALL histwrite(nid_tra,"v",itra,zx_tmp_3d,
961     .                                   iim*(jjm+1)*klev,ndex)
962
963      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,sh, zx_tmp_3d)
964      CALL histwrite(nid_tra,"h2o",itra,zx_tmp_3d,
965     .                                   iim*(jjm+1)*klev,ndex)
966
967      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pdel, zx_tmp_3d)
968      CALL histwrite(nid_tra,"pdel",itra,zx_tmp_3d,
969     .                                   iim*(jjm+1)*klev,ndex)
970
971      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pplay, zx_tmp_3d)
972      CALL histwrite(nid_tra,"pmid",itra,zx_tmp_3d,
973     .                                   iim*(jjm+1)*klev,ndex)
974
975#ifdef INCA_CH4
976#ifdef INCAINFO
977      DO it=1, phtcnt
978      WRITE(str2,'(i2.2)') it
979      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,jrates(1,1,it),
980     .     zx_tmp_3d)
981      CALL histwrite(nid_tra,"j"//str2,itra,zx_tmp_3d,
982     .                                   iim*(jjm+1)*klev,ndex)
983      ENDDO
984
985      DO it=1, hetcnt
986      WRITE(str2,'(i2.2)') it
987      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,hrates(1,1,it),
988     .     zx_tmp_3d)
989      CALL histwrite(nid_tra,"w"//str2,itra,zx_tmp_3d,
990     .                                   iim*(jjm+1)*klev,ndex)
991      ENDDO
992
993 
994      DO it=1, extcnt
995      WRITE(str2,'(i2.2)') it
996
997
998      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,extflx(1,1,it),
999     .     zx_tmp_3d)
1000      CALL histwrite(nid_tra,"ext"//str2,itra,zx_tmp_3d,
1001     .                                   iim*(jjm+1)*klev,ndex)
1002      ENDDO
1003
1004      DO it=1, nfs
1005      WRITE(str2,'(i2.2)') it
1006      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,invariants(1,1,it),
1007     .     zx_tmp_3d)
1008      CALL histwrite(nid_tra,"INV"//str2,itra,zx_tmp_3d,
1009     .                                   iim*(jjm+1)*klev,ndex)
1010      ENDDO
1011
1012
1013#else
1014      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,jrates(1,1,2),
1015     .     zx_tmp_3d)
1016      CALL histwrite(nid_tra,"jO3",itra,zx_tmp_3d,
1017     .                                   iim*(jjm+1)*klev,ndex)
1018
1019      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,jrates(1,1,4),
1020     .     zx_tmp_3d)
1021      CALL histwrite(nid_tra,"jNO2",itra,zx_tmp_3d,
1022     .                                   iim*(jjm+1)*klev,ndex)
1023
1024      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,jrates(1,1,13),
1025     .     zx_tmp_3d)
1026      CALL histwrite(nid_tra,"jH2O2",itra,zx_tmp_3d,
1027     .                                   iim*(jjm+1)*klev,ndex)
1028
1029      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,hrates(1,1,1),
1030     .     zx_tmp_3d)
1031      CALL histwrite(nid_tra,"wHNO3",itra,zx_tmp_3d,
1032     .                                   iim*(jjm+1)*klev,ndex)
1033
1034      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,krates(1,1,1),
1035     .     zx_tmp_3d)
1036      CALL histwrite(nid_tra,"kN2O5",itra,zx_tmp_3d,
1037     .                                   iim*(jjm+1)*klev,ndex)
1038
1039      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,extflx(1,1,1),
1040     .     zx_tmp_3d)
1041      CALL histwrite(nid_tra,"LghtNO",itra,zx_tmp_3d,
1042     .                                   iim*(jjm+1)*klev,ndex)
1043#endif
1044      DO it=1, grpcnt
1045      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,nas(1,1,it),zx_tmp_3d)
1046      zx_tmp_3d = zx_tmp_3d * dry_mass / nadv_mass(it)
1047      CALL histwrite(nid_tra,grpsym(it),itra,zx_tmp_3d,
1048     .                                   iim*(jjm+1)*klev,ndex)
1049      ENDDO
1050#endif
1051
1052#endif
1053#ifdef INCA_AER
1054C   3d FIELDS
1055
1056      it = id_CIDUSTM
1057       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,scavcoef_st(1,1,it),
1058     .                  zx_tmp_3d)
1059       CALL histwrite(nid_tra,"scavcoef_st",itra,zx_tmp_3d,
1060     .                  iim*(jjm+1)*klev,ndex)
1061
1062       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,scavcoef_cv(1,1,it),
1063     .                  zx_tmp_3d)
1064       CALL histwrite(nid_tra,"scavcoef_cv",itra,zx_tmp_3d,
1065     .                  iim*(jjm+1)*klev,ndex)
1066
1067#endif
1068
1069      DO it=1,nqmax
1070      IF (it.LE.99) THEN
1071       WRITE(str2,'(i2.2)') it
1072
1073#ifdef INCA
1074      IF ( prt_flag_ts(it) == 0 ) CYCLE
1075C champs 2D
1076      CALL gr_fi_ecrit(1, klon,iim,jjm+1, eflux(1,it),zx_tmp_2d)
1077      CALL histwrite(nid_tra,"Emi_"//solsym(it),itra,zx_tmp_2d,
1078     .     iim*(jjm+1),ndex)
1079
1080      CALL gr_fi_ecrit(1, klon,iim,jjm+1, dvel(1,it),zx_tmp_2d)
1081      CALL histwrite(nid_tra,"Dep_"//solsym(it),itra,zx_tmp_2d,
1082     .     iim*(jjm+1),ndex)
1083#ifdef INCA_AER
1084        IF (solsym(it) .eq. 'CIDUSTM'.or.solsym(it) .eq. 'CIN'
1085     .  .or.solsym(it) .eq. 'CSSSM'  .or.solsym(it) .eq. 'CSN'
1086     .  .or.solsym(it) .eq. 'ASSSM'  .or.solsym(it) .eq. 'ASN' ) THEN
1087      CALL gr_fi_ecrit(1, klon,iim,jjm+1,sflux(1,it),zx_tmp_2d)
1088      CALL histwrite(nid_tra,"Sed_"//solsym(it),itra,zx_tmp_2d,
1089     .     iim*(jjm+1),ndex)
1090        endif
1091        IF (solsym(it) .eq. 'CIDUSTM') THEN
1092      CALL gr_fi_ecrit(1, klon,iim,jjm+1,tausum(1),zx_tmp_2d)
1093      CALL histwrite(nid_tra,"OD_"//solsym(it),itra,zx_tmp_2d,
1094     .     iim*(jjm+1),ndex)
1095        endif
1096
1097#endif
1098C champs 3D
1099       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
1100
1101       !Prefer vmr to mmr for transported species
1102       if( adv_mass(it) /= 0. ) then
1103       zx_tmp_3d = zx_tmp_3d * dry_mass / adv_mass(it)
1104       else
1105#ifdef INCA_CH4
1106       if ( solsym(it) == 'OX' ) then
1107       zx_tmp_3d = zx_tmp_3d * dry_mass / nadv_mass(id_o3)
1108       end if
1109#endif
1110       end if
1111
1112       CALL histwrite(nid_tra,solsym(it),itra,zx_tmp_3d,
1113     .                                   iim*(jjm+1)*klev,ndex)
1114#else
1115
1116C champs 3D
1117       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
1118       CALL histwrite(nid_tra,"tr"//str2,itra,zx_tmp_3d,
1119     .                                   iim*(jjm+1)*klev,ndex)
1120       IF (lessivage) THEN
1121       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,flestottr(1,1,it),zx_tmp_3d)
1122       CALL histwrite(nid_tra,"fl"//str2,itra,zx_tmp_3d,
1123     .                                   iim*(jjm+1)*klev,ndex)
1124       ENDIF
1125#endif
1126      ELSE
1127         PRINT*, "Trop de traceurs"
1128         CALL abort
1129      ENDIF
1130      ENDDO
1131
1132#ifdef INCA_CH4
1133      CALL gr_fi_ecrit(1, klon,iim,jjm+1, o3_tr_col(1), zx_tmp_2d)
1134      CALL histwrite(nid_tra,"O3_column",itra,zx_tmp_2d,
1135     .     iim*(jjm+1),ndex)
1136
1137      CALL gr_fi_ecrit(1, klon,iim,jjm+1, co_tr_col(1), zx_tmp_2d)
1138      CALL histwrite(nid_tra,"CO_column",itra,zx_tmp_2d,
1139     .     iim*(jjm+1),ndex)
1140
1141      CALL gr_fi_ecrit(1, klon,iim,jjm+1, ch4_tr_col(1), zx_tmp_2d)
1142      CALL histwrite(nid_tra,"CH4_column",itra,zx_tmp_2d,
1143     .     iim*(jjm+1),ndex)
1144
1145      CALL gr_fi_ecrit(1, klon,iim,jjm+1, no2_tr_col(1), zx_tmp_2d)
1146      CALL histwrite(nid_tra,"NO2_column",itra,zx_tmp_2d,
1147     .     iim*(jjm+1),ndex)
1148
1149      CALL gr_fi_ecrit(1, klon,iim,jjm+1, o3_st_flx(1), zx_tmp_2d)
1150      CALL histwrite(nid_tra,"O3_ste",itra,zx_tmp_2d,
1151     .     iim*(jjm+1),ndex)
1152
1153      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,o3_prod(1,1),
1154     .     zx_tmp_3d)
1155      CALL histwrite(nid_tra,"O3_prod",itra,zx_tmp_3d,
1156     .                                   iim*(jjm+1)*klev,ndex)
1157
1158      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,o3_loss(1,1),
1159     .     zx_tmp_3d)
1160      CALL histwrite(nid_tra,"O3_loss",itra,zx_tmp_3d,
1161     .                                   iim*(jjm+1)*klev,ndex)
1162
1163!     ... Special section for daytime averaging
1164!       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,day_cnt(1,1),
1165!    .       zx_tmp_3d)
1166!       CALL histwrite(nid_tra,"day_cnt",itra,zx_tmp_3d,
1167!    .                                  iim*(jjm+1)*klev,ndex)
1168!       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,no_daytime(1,1),
1169!    .       zx_tmp_3d)
1170!       CALL histwrite(nid_tra,"NO_day",itra,zx_tmp_3d,
1171!    .                                  iim*(jjm+1)*klev,ndex)
1172
1173#endif
1174
1175      if (ok_sync) call histsync(nid_tra)
1176
1177      if (lafin) then
1178!DH      print*, 'c est la fin de la physique'
1179         open (99,file='restarttrac',  form='formatted')
1180         do i=1,klon
1181             write(99,*) trs(i,1)
1182         enddo
1183         PRINT*, 'Ecriture du fichier restarttrac'
1184         close(99)
1185      else
1186!DH      print*, 'physique pas fini'
1187      endif
1188
1189
1190      RETURN
1191      END
Note: See TracBrowser for help on using the repository browser.