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

Last change on this file since 2086 was 348, checked in by lmdz, 23 years ago

Regle le probleme de decalage de 1 jour au debut de chaque simulation
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.5 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 idayref
213      INTEGER nid_tra
214      SAVE nid_tra
215c     REAL x(klon,klev,nbtr+2) ! traceurs
216      INTEGER ndex(1)
217      REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
218      REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
219      INTEGER itra
220      SAVE itra                   ! compteur pour la physique
221C
222C Variables liees a l'ecriture de la bande histoire : phytrac.nc
223c
224      INTEGER ecrit_tra
225      SAVE ecrit_tra   
226      logical ok_sync
227      parameter (ok_sync = .true.)
228C
229C nature du traceur
230c
231      logical aerosol(nbtr)  ! Nature du traceur
232c                            ! aerosol(it) = true  => aerosol
233c                            ! aerosol(it) = false => gaz
234c                            ! nat_trac(it) = 1. aerosolc
235      logical clsol(nbtr)    ! clsol(it) = true => CL sol calculee
236      logical radio(nbtr)    ! radio(it)=true => decroisssance radioactive
237      save aerosol,clsol,radio
238C
239c======================================================================
240c
241c Declaration des procedures appelees
242c
243c--modif convection tiedtke
244      INTEGER i, k, it,itap
245        save itap
246      REAL delp(klon,klev)
247c--end modif
248c
249c Variables liees a l'ecriture de la bande histoire physique
250c
251c Variables locales pour effectuer les appels en serie
252c----------------------------------------------------
253c
254      REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
255      REAL d_tr_cl(klon,klev) ! tendance de traceurs  couche limite
256      REAL d_tr_cv(klon,klev) ! tendance de traceurs  convection
257      REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
258c                                   ! radioactive du rn - > pb
259      REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
260c                                          ! par impaction
261      REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
262c                                          ! par nucleation
263      REAL fluxrn(klon,klev)
264      REAL fluxpb(klon,klev)
265      REAL pbimpa(klon,klev)
266      REAL pbnucl(klon,klev)
267      REAL rn(klon,klev)
268      REAL pb(klon,klev)
269      REAL flestottr(klon,klev,nbtr) ! flux de lessivage
270c                                    ! dans chaque couche
271
272C
273      character*20 modname
274      character*80 abort_message
275c
276c   Controles
277c-------------
278      logical first,couchelimite,convection,lessivage,sorties,
279     s        rnpb,inirnpb
280      save first,couchelimite,convection,lessivage,sorties,
281     s     inirnpb
282      data first,couchelimite,convection,lessivage,sorties
283     s     /.true.,.true.,.true.,.true.,.true./
284
285#ifdef INCA
286      INTEGER           :: ncsec
287      INTEGER           :: prt_flag_ts(nbtr)=(/1,1,1
288#ifdef INCA_CH4
289     .                                              ,0,0,1,1,1,1,1,
290     .                                         0,1,0,0,0,0,0,1,0,0,
291     .                                         0,1,1,1,1,0,1,1,1,0,
292     .                                         1,1,1,1,1,1,1,1,1,1,
293     .                                         1,0,0
294#endif
295#ifdef INCA_AER
296     .                                              ,1,1,1,1,1,1
297#endif
298     .                                         /)
299      REAL, PARAMETER   :: dry_mass = 28.966
300      REAL, POINTER     :: hbuf(:)
301      REAL, ALLOCATABLE :: obuf(:)
302      REAL              :: calday
303      REAL              :: pdel(klon,klev)
304      REAL              :: dummy(klon,klev) = 0.
305#endif
306
307c
308c======================================================================
309c         ecrit_tra = NINT(86400./pdtphys *1.0)  ! tous les jours
310         modname='phytrac'
311         ecrit_tra = NINT(86400./pdtphys *ecritphy)   
312!DH      print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
313         ps(:) = paprs(:,1)
314
315c        print*,'DANS PHYTRAC debutphy=',debutphy
316
317         if (debutphy) then
318
319          print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
320          ecrit_tra = NINT(86400./pdtphys/2.) ! tous les 12H
321c         ecrit_tra = NINT(86400./pdtphys) ! tous les 24H
322
323         if(nbtr.lt.nqmax) then
324           print*,'NQMAX=',nqmax
325           print*,'NBTR=',nbtr
326           abort_message='See above'
327           call abort_gcm(modname,abort_message,1)
328         endif
329
330         inirnpb=rnpb
331!DH      PRINT*, 'La frequence de sortie traceurs est  ', ecrit_tra
332
333C         
334         idayref = 1
335         CALL ymds2ju(anne_ini, 1, idayref, 0.0, zjulian)
336         itra = NINT(FLOAT(day_ini)*86400./pdtphys)
337         itap = NINT(FLOAT(day_ini)*86400./pdtphys)
338!        itra=0
339!        itap=0
340!        zjulian = zjulian + day_ini
341 
342         CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlon,zx_lon)
343         DO i = 1, iim
344            zx_lon(i,1) = xlon(i+1)
345            zx_lon(i,jjm+1) = xlon(i+1)
346         ENDDO
347         DO ll=1,klev
348            znivsig(ll)=float(ll)
349         ENDDO
350         CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlat,zx_lat)
351         CALL histbeg("histrac", iim,zx_lon, jjm+1,zx_lat,
352     .                 1,iim,1,jjm+1, itra, zjulian, pdtphys,
353     .                 nhori, nid_tra)
354
355         call histvert(nid_tra, "presnivs", "presnivs", "mb",
356     .                 klev, presnivs, nvert)
357#ifdef INCA
358!        call histvert(nid_tra, "ap", "Hybrid A parameter", "-",
359!    .                 klev+1, ap, nverta)
360!        call histvert(nid_tra, "bp", "Hybrid B parameter", "-",
361!    .                 klev+1, bp, nvertb)
362#endif
363
364         zsto = pdtphys
365         zout = pdtphys * FLOAT(ecrit_tra)
366c
367         CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
368     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
369     .                "once",  zsto,zout)
370c
371         CALL histdef(nid_tra, "aire", "Grid area", "-",
372     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
373     .                "once",  zsto,zout)
374
375#ifdef INCA
376         CALL histdef(nid_tra, "ps", "Surface pressure", "Pa",
377     .                iim,jjm+1,nhori, 1,1,1,-99, 32,
378     .                "ave(X)", zsto,zout)
379
380         CALL histdef(nid_tra, "ptrop", "Tropopause pressure", "Pa",
381     .                iim,jjm+1,nhori, 1,1,1,-99, 32,
382     .                "ave(X)", zsto,zout)
383
384         CALL histdef(nid_tra, "temp", "Air temperature", "K",
385     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
386     .                "ave(X)", zsto,zout)
387
388         CALL histdef(nid_tra, "u", "zonal wind component", "m/s",
389     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
390     .                "ave(X)", zsto,zout)
391
392         CALL histdef(nid_tra, "v", "zonal wind component", "m/s",
393     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
394     .                "ave(X)", zsto,zout)
395
396         CALL histdef(nid_tra, "h2o", "Specific Humidity", "MMR",
397     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
398     .                "ave(X)", zsto,zout)
399
400         CALL histdef(nid_tra, "pmid", "Pressure", "Pa",
401     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
402     .                "ave(X)", zsto,zout)
403
404         CALL histdef(nid_tra, "pdel", "Delta Pressure", "Pa",
405     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
406     .                "ave(X)", zsto,zout)
407#ifdef INCA_CH4
408#ifdef INCAINFO
409         DO it=1, phtcnt
410         WRITE(str2,'(i2.2)') it
411         CALL histdef(nid_tra, "j"//str2,"j"//str2, "CM-3 S-1",
412     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
413     .                "ave(X)", zsto,zout)
414         ENDDO
415         DO it=1, hetcnt
416         WRITE(str2,'(i2.2)') it
417         CALL histdef(nid_tra, "w"//str2,"w"//str2, "S-1",
418     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
419     .                "ave(X)", zsto,zout)
420         ENDDO
421
422         DO it=1, extcnt
423         WRITE(str2,'(i2.2)') it
424        CALL histdef(nid_tra, "ext"//str2,"ext"//str2, "CM-3 S-1",
425     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
426     .                "ave(X)", zsto,zout)
427         ENDDO
428
429         DO it=1, nfs
430         WRITE(str2,'(i2.2)') it
431         CALL histdef(nid_tra, "INV"//str2, "INV"//str2, "CM-3",
432     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
433     .                "ave(X)", zsto,zout)
434         ENDDO
435
436#else
437         CALL histdef(nid_tra, "jO3","jO3", "CM-3 S-1",
438     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
439     .                "ave(X)", zsto,zout)
440         CALL histdef(nid_tra, "jNO2","jNO2", "CM-3 S-1",
441     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
442     .                "ave(X)", zsto,zout)
443         CALL histdef(nid_tra, "jH2O2","jH2O2", "CM-3 S-1",
444     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
445     .                "ave(X)", zsto,zout)
446         CALL histdef(nid_tra, "wHNO3","wHNO3", "S-1",
447     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
448     .                "ave(X)", zsto,zout)
449         CALL histdef(nid_tra, "kN2O5", "kN2O5","CM-3 S-1",
450     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
451     .                "ave(X)", zsto,zout)
452         CALL histdef(nid_tra, "LghtNO","LghtNO", "CM-3 S-1",
453     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
454     .                "ave(X)", zsto,zout)
455#endif
456
457         DO it=1, grpcnt
458         CALL histdef(nid_tra, grpsym(it), grpsym(it), "VMR",
459     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
460     .                "ave(X)", zsto,zout)
461         ENDDO
462#endif
463#endif
464
465#ifdef INCA_AER
466C   3d FIELDS
467
468        CALL histdef(nid_tra, "scavcoef_st","scavcoef_st", "S-1",
469     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
470     .                "ave(X)", zsto,zout)
471        CALL histdef(nid_tra, "scavcoef_cv","scavcoef_cv", "S-1",
472     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
473     .                "ave(X)", zsto,zout)
474#endif
475 
476         DO it=1,nqmax
477         IF (it.LE.99) THEN
478         WRITE(str2,'(i2.2)') it
479C champ 2D
480
481#ifdef INCA
482         IF ( prt_flag_ts(it) == 0 ) CYCLE
483
484         CALL histdef(nid_tra, "Emi_"//solsym(it), "Emi_"//solsym(it),
485     .           "kg/m2/s", iim,jjm+1,nhori, 1,1,1, -99, 32,
486     .           "ave(X)", zsto,zout)
487         CALL histdef(nid_tra, "Dep_"//solsym(it), "Dep_"//solsym(it),
488     .           "cm/s", iim,jjm+1,nhori, 1,1,1, -99, 32,
489     .           "ave(X)", zsto,zout)
490#ifdef INCA_AER
491        IF (solsym(it) .eq. 'CIDUSTM'.or.solsym(it) .eq. 'CIN'
492     .  .or.solsym(it) .eq. 'CSSSM'  .or.solsym(it) .eq. 'CSN'
493     .  .or.solsym(it) .eq. 'ASSSM'  .or.solsym(it) .eq. 'ASN' ) THEN
494         CALL histdef(nid_tra, "Sed_"//solsym(it), "Sed_"//solsym(it),
495     .           "kg/m2", iim,jjm+1,nhori, 1,1,1, -99, 32,
496     .           "ave(X)", zsto,zout)
497        endif
498        IF (solsym(it) .eq. 'CIDUSTM') THEN
499         CALL histdef(nid_tra, "OD_"//solsym(it), "OD_"//solsym(it),
500     .           "opt. depth", iim,jjm+1,nhori, 1,1,1, -99, 32,
501     .           "ave(X)", zsto,zout)
502        endif
503
504#endif
505         CALL histdef(nid_tra, solsym(it), solsym(it), "VMR",
506     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
507     .                "ave(X)", zsto,zout)
508#else
509         CALL histdef(nid_tra, "tr"//str2, "Tracer No."//str2, "U/kga",
510     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
511     .                "ave(X)", zsto,zout)
512         IF (lessivage) THEN
513         CALL histdef(nid_tra, "fl"//str2, "Tracer Flux No."//str2,
514     .                "U/m2/s",iim,jjm+1,nhori, klev,1,klev,nvert, 32,
515     .                "ave(X)", zsto,zout)
516         ENDIF
517#endif
518         ELSE
519         PRINT*, "Trop de traceurs"
520         CALL abort
521         ENDIF
522         ENDDO
523
524#ifdef INCA_CH4
525         CALL histdef(nid_tra, "O3_column", "O3_column",
526     .           "DU", iim,jjm+1,nhori, 1,1,1, -99, 32,
527     .           "ave(X)", zsto,zout)
528         CALL histdef(nid_tra, "CO_column", "CO_column",
529     .           "10^18 CM-2", iim,jjm+1,nhori, 1,1,1, -99, 32,
530     .           "ave(X)", zsto,zout)
531         CALL histdef(nid_tra, "CH4_column", "CH4_column",
532     .           "10^18 CM-2", iim,jjm+1,nhori, 1,1,1, -99, 32,
533     .           "ave(X)", zsto,zout)
534         CALL histdef(nid_tra, "NO2_column", "NO2_column",
535     .           "10^15 CM-2", iim,jjm+1,nhori, 1,1,1, -99, 32,
536     .           "ave(X)", zsto,zout)
537         CALL histdef(nid_tra, "O3_ste", "O3_ste",
538     .           "CM-2 S-1", iim,jjm+1,nhori, 1,1,1, -99, 32,
539     .           "ave(X)", zsto,zout)
540         CALL histdef(nid_tra, "O3_prod", "O3_prod", "CM-3 S-1",
541     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
542     .                "ave(X)", zsto,zout)
543         CALL histdef(nid_tra, "O3_loss", "O3_loss", "CM-3 S-1",
544     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
545     .                "ave(X)", zsto,zout)
546
547!        Special variables for daytime averaging
548!        CALL histdef(nid_tra, "day_cnt", "day_cnt", "-",
549!    .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
550!    .                "t_sum(X)", zsto,zout)
551!        CALL histdef(nid_tra, "NO_day", "NO_day", "VMR",
552!    .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
553!    .                "t_sum(X)", zsto,zout)
554
555#endif
556
557         CALL histend(nid_tra)
558         ndex(1) = 0
559c
560         i = NINT(zout/zsto)
561         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
562         CALL histwrite(nid_tra,"phis",i,zx_tmp_2d,iim*(jjm+1),ndex)
563C
564         i = NINT(zout/zsto)
565         CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
566         CALL histwrite(nid_tra,"aire",i,zx_tmp_2d,iim*(jjm+1),ndex)
567
568c======================================================================
569c   Initialisation de certaines variables pour le Rn et le Pb
570c======================================================================
571
572c Initialisation du traceur dans le sol (couche limite radonique)
573c
574c        print*,'valeur de debut dans phytrac :',debutphy
575         do it=1,nqmax
576            do i=1,klon
577               trs(i,it) = 0.
578            enddo
579         END DO
580
581         open (99,file='starttrac',status='old',
582     .         err=999,form='formatted')
583         read(99,*) (trs(i,1),i=1,klon)
584999      close(99)
585!DH      print*, 'apres starttrac'
586
587c Initialisation de la fraction d'aerosols lessivee
588c
589         DO it=1,nqmax
590           DO k = 1, nlev
591              DO i = 1, klon
592                 d_tr_lessi_impa(i,k,it) =  0.0
593                 d_tr_lessi_nucl(i,k,it) =  0.0
594             ENDDO
595           ENDDO
596         ENDDO
597c
598c Initialisation de la nature des traceurs
599c
600         DO it = 1, nqmax
601            aerosol(it) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
602            radio(it) = .FALSE.    ! Par defaut pas de passage par radiornpb
603            clsol(it) = .FALSE.  ! Par defaut couche limite avec flux prescrit
604         ENDDO
605c
606      ENDIF  ! fin debutphy
607c Initialisation du traceur dans le sol (couche limite radonique)
608      if(inirnpb) THEN
609c
610         radio(1)= .true.
611         radio(2)= .true.
612         clsol(1)= .true.
613         clsol(2)= .true.
614         aerosol(2) = .TRUE. ! le Pb est un aerosol
615c
616         call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
617     .                   ,vdeptr,scavtr)
618         inirnpb=.false.
619      endif
620
621#ifdef INCA
622!======================================================================
623!     Chimie
624!======================================================================
625
626        calday = FLOAT(julien) + gmtime
627        ncsec  = NINT (86400.*gmtime)
628
629        DO k = 1, nlev
630        pdel(:,k) = paprs(:,k) - paprs (:,k+1)
631        END DO
632
633#ifdef INCAINFO
634        PRINT *, 'CHEMMAIN @ ', calday, ' ... '
635        DO it = 1, nbtr
636        PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
637     $                       MAXVAL(tr_seri(:,:,it))
638      END DO
639#endif
640
641
642#ifdef INCA_AER
643        CALL aerosolmain (tr_seri,
644     $                 pdtphys,
645     $                 pplay,
646     $                 prfl,
647     $                 pmflxr,
648     $                 psfl,
649     $                 pmflxs,
650     $                 pmfu,
651     $                 itop_con,
652     $                 ibas_con,
653     $                 pphi,
654     $                 paire,
655     $                 nstep)
656#endif
657
658        CALL chemmain (tr_seri,
659     $                 nas,
660     $                 nstep,
661     $                 calday,
662     $                 julien,
663     $                 ncsec,
664     $                 1,
665     $                 pdtphys,
666     $                 paprs(1,1),
667     $                 pplay,
668     $                 pdel,
669     $                 pctsrf(1,3),
670     $                 ftsol,
671     $                 albsol,
672     $                 pphi,
673     $                 pphis,
674     $                 cldfra,
675     $                 rneb,
676     $                 diafra,
677     $                 itop_con,
678     $                 cldliq,
679     $                 prfl,
680     $                 pmflxr,
681     $                 psfl,
682     $                 pmflxs,
683     $                 pmfu,
684     $                 flxmass_w,
685     $                 t_seri,
686     $                 sh,
687     $                 rh,
688     $                 .false.,
689     $                 hbuf,
690     $                 obuf,
691     $                 iip1,
692     $                 jjp1)
693#ifdef INCAINFO
694      PRINT *, 'OK.'
695      DO it = 1, nbtr
696      PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
697     $                     MAXVAL(tr_seri(:,:,it))
698      END DO
699#endif
700#endif
701C
702
703c======================================================================
704c   Calcul de l'effet de la convection
705c======================================================================
706!DH   print*,'Avant convection'
707      do it=1,nqmax
708         WRITE(itn,'(i1)') it
709c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
710      enddo
711
712      if (convection) then
713
714!DH   print*,'Pas de temps dans phytrac : ',pdtphys
715      DO it=1, nqmax
716#ifdef INCA
717      IF ( conv_flg(it) == 0 ) CYCLE
718#endif
719      CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
720     .            pplay, paprs, tr_seri(1,1,it), d_tr_cv)
721      DO k = 1, nlev
722      DO i = 1, klon
723         tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k)
724      ENDDO
725      ENDDO
726c     WRITE(itn,'(i2)') it
727#ifdef INCA
728      CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//solsym(it))
729#else
730      CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
731#endif
732      ENDDO
733c     print*,'apres nflxtr'
734
735      endif ! convection
736c        print*,'Apres convection'
737c      do it=1,nqmax
738c         WRITE(itn,'(i1)') it
739c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
740c      enddo
741
742c======================================================================
743c   Calcul de l'effet de la couche limite
744c======================================================================
745c       print *,'Avant couchelimite'
746c      do it=1,nqmax
747c         WRITE(itn,'(i1)') it
748c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL  '//itn)
749c      enddo
750
751      if (couchelimite) then
752
753      DO k = 1, nlev
754      DO i = 1, klon
755         delp(i,k) = paprs(i,k)-paprs(i,k+1)
756      ENDDO
757      ENDDO
758
759C maf modif pour tenir compte du cas rnpb + traceur
760      DO it=1, nqmax
761#ifdef INCA
762      IF ( pbl_flg(it) == 0 ) CYCLE
763#endif
764C     print *,'it',it,clsol(it)
765      if (clsol(it)) then  ! couche limite avec quantite dans le sol calculee
766          CALL cltracrn(it, pdtphys, yu1, yv1,
767     e                    coefh,t_seri,ftsol,pctsrf,
768     e                    tr_seri(1,1,it),trs(1,it),
769     e                    paprs, pplay, delp,
770     e                    masktr(1,it),fshtr(1,it),hsoltr(it),
771     e                    tautr(it),vdeptr(it),
772     e                    xlat,
773     s                    d_tr_cl,d_trs)
774          DO k = 1, nlev
775            DO i = 1, klon
776              tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k)
777            ENDDO
778          ENDDO
779c
780c Traceur ds sol
781c
782          DO i = 1, klon
783            trs(i,it) = trs(i,it) + d_trs(i)
784          END DO
785C
786C maf provisoire suppression des prints
787C         WRITE(itn,'(i1)') it
788C         CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
789      else ! couche limite avec flux prescrit
790
791#ifdef INCA
792        DO k =  1, klon
793          source(k) = eflux(k,it)-dflux(k,it)
794        END DO
795#else
796Cmaf provisoire source / traceur a creer
797        DO i=1, klon
798          source(i) = 0.0 ! pas de source, pour l'instant
799        ENDDO
800#endif
801          CALL cltrac(pdtphys, coefh,t_seri,
802     s               tr_seri(1,1,it), source,
803     e               paprs, pplay, delp,
804     s               d_tr )
805            DO k = 1, nlev
806               DO i = 1, klon
807                  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr(i,k)
808               ENDDO
809            ENDDO
810Cmaf provisoire suppression des prints
811Cmaf          WRITE(itn,'(i1)') it
812cmaf          CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
813      endif
814      ENDDO
815c
816      endif ! couche limite
817
818c      print*,'Apres couchelimite'
819c      do it=1,nqmax
820c         WRITE(itn,'(i1)') it
821c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL  '//itn)
822c      enddo
823
824c======================================================================
825c   Calcul de l'effet du puits radioactif
826c======================================================================
827
828C MAF il faudrait faire une modification pour passer dans radiornpb
829c si radio=true mais pour l'instant radiornpb propre au cas rnpb
830      if(rnpb) then
831        print *, 'decroissance radiactive activee'
832        call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
833C
834        DO it=1,nqmax
835            if(radio(it)) then
836            DO k = 1, nlev
837               DO i = 1, klon
838                  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
839               ENDDO
840            ENDDO
841            WRITE(itn,'(i1)') it
842            CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
843            endif
844        ENDDO
845c
846      endif ! rnpb decroissance  radioactive
847
848c======================================================================
849c   Calcul de l'effet de la precipitation
850c======================================================================
851
852!DH   print*,'LESSIVAGE =',lessivage
853      IF (lessivage) THEN
854
855c     print*,'avant lessivage'
856
857       DO it = 1, nqmax
858           DO k = 1, nlev
859              DO i = 1, klon
860               d_tr_lessi_nucl(i,k,it) = 0.0
861               d_tr_lessi_impa(i,k,it) = 0.0
862               flestottr(i,k,it) = 0.0
863              ENDDO
864           ENDDO
865       ENDDO
866c
867c tendance des aerosols nuclees et impactes
868c
869       DO it = 1, nqmax
870         IF (aerosol(it)) THEN
871           DO k = 1, nlev
872              DO i = 1, klon
873               d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
874     s                  ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
875               d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
876     s                  ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
877              ENDDO
878           ENDDO
879         ENDIF
880       ENDDO
881c
882c Mises a jour des traceurs + calcul des flux de lessivage
883c Mise a jour due a l'impaction et a la nucleation
884c
885c      call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
886c      call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
887c      call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
888       DO it = 1, nqmax
889c         print*,'IT=',it,aerosol(it)
890         IF (aerosol(it)) THEN
891c           print*,'IT=',it,' On lessive'
892           DO k = 1, nlev
893              DO i = 1, klon
894               tr_seri(i,k,it)=tr_seri(i,k,it)
895     s         *frac_impa(i,k)*frac_nucl(i,k)
896              ENDDO
897           ENDDO
898         ENDIF
899       ENDDO
900c      call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
901c
902c Flux lessivage total
903c
904      DO it = 1, nqmax
905           DO k = 1, nlev
906            DO i = 1, klon
907               flestottr(i,k,it) = flestottr(i,k,it) -
908     s                   ( d_tr_lessi_nucl(i,k,it)   +
909     s                     d_tr_lessi_impa(i,k,it) ) *
910     s                   ( paprs(i,k)-paprs(i,k+1) ) /
911     s                   (RG * pdtphys)
912            ENDDO
913           ENDDO
914c
915Cmaf suppression provisoire
916Cmaf        WRITE(itn,'(i1)') it
917Cmaf    CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
918      ENDDO
919c
920c     print*,'apres lessivage'
921      ENDIF
922Cc
923      DO k = 1, nlev
924         DO i = 1, klon
925            fluxrn(i,k) = flestottr(i,k,1)
926            fluxpb(i,k) = flestottr(i,k,2)
927            rn(i,k) = tr_seri(i,k,1)
928            pb(i,k) = tr_seri(i,k,2)
929            pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
930            pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
931         ENDDO
932      ENDDO
933
934
935
936      itra=itra+1
937      ndex(1) = 0
938
939      i = NINT(zout/zsto)
940      CALL gr_fi_ecrit(1, klon,iim,jjm+1, paire,zx_tmp_2d)
941      CALL histwrite(nid_tra,"aire",itra,zx_tmp_2d,
942     .     iim*(jjm+1),ndex)
943
944#ifdef INCA
945      CALL gr_fi_ecrit(1, klon,iim,jjm+1, ps,zx_tmp_2d)
946      CALL histwrite(nid_tra,"ps",itra,zx_tmp_2d,
947     .     iim*(jjm+1),ndex)
948
949      CALL gr_fi_ecrit(1, klon,iim,jjm+1, ptrop,zx_tmp_2d)
950      CALL histwrite(nid_tra,"ptrop",itra,zx_tmp_2d,
951     .     iim*(jjm+1),ndex)
952
953      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,t_seri, zx_tmp_3d)
954      CALL histwrite(nid_tra,"temp",itra,zx_tmp_3d,
955     .                                   iim*(jjm+1)*klev,ndex)
956
957      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,u, zx_tmp_3d)
958      CALL histwrite(nid_tra,"u",itra,zx_tmp_3d,
959     .                                   iim*(jjm+1)*klev,ndex)
960
961      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,v, zx_tmp_3d)
962      CALL histwrite(nid_tra,"v",itra,zx_tmp_3d,
963     .                                   iim*(jjm+1)*klev,ndex)
964
965      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,sh, zx_tmp_3d)
966      CALL histwrite(nid_tra,"h2o",itra,zx_tmp_3d,
967     .                                   iim*(jjm+1)*klev,ndex)
968
969      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pdel, zx_tmp_3d)
970      CALL histwrite(nid_tra,"pdel",itra,zx_tmp_3d,
971     .                                   iim*(jjm+1)*klev,ndex)
972
973      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pplay, zx_tmp_3d)
974      CALL histwrite(nid_tra,"pmid",itra,zx_tmp_3d,
975     .                                   iim*(jjm+1)*klev,ndex)
976
977#ifdef INCA_CH4
978#ifdef INCAINFO
979      DO it=1, phtcnt
980      WRITE(str2,'(i2.2)') it
981      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,jrates(1,1,it),
982     .     zx_tmp_3d)
983      CALL histwrite(nid_tra,"j"//str2,itra,zx_tmp_3d,
984     .                                   iim*(jjm+1)*klev,ndex)
985      ENDDO
986
987      DO it=1, hetcnt
988      WRITE(str2,'(i2.2)') it
989      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,hrates(1,1,it),
990     .     zx_tmp_3d)
991      CALL histwrite(nid_tra,"w"//str2,itra,zx_tmp_3d,
992     .                                   iim*(jjm+1)*klev,ndex)
993      ENDDO
994
995 
996      DO it=1, extcnt
997      WRITE(str2,'(i2.2)') it
998
999
1000      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,extflx(1,1,it),
1001     .     zx_tmp_3d)
1002      CALL histwrite(nid_tra,"ext"//str2,itra,zx_tmp_3d,
1003     .                                   iim*(jjm+1)*klev,ndex)
1004      ENDDO
1005
1006      DO it=1, nfs
1007      WRITE(str2,'(i2.2)') it
1008      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,invariants(1,1,it),
1009     .     zx_tmp_3d)
1010      CALL histwrite(nid_tra,"INV"//str2,itra,zx_tmp_3d,
1011     .                                   iim*(jjm+1)*klev,ndex)
1012      ENDDO
1013
1014
1015#else
1016      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,jrates(1,1,2),
1017     .     zx_tmp_3d)
1018      CALL histwrite(nid_tra,"jO3",itra,zx_tmp_3d,
1019     .                                   iim*(jjm+1)*klev,ndex)
1020
1021      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,jrates(1,1,4),
1022     .     zx_tmp_3d)
1023      CALL histwrite(nid_tra,"jNO2",itra,zx_tmp_3d,
1024     .                                   iim*(jjm+1)*klev,ndex)
1025
1026      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,jrates(1,1,13),
1027     .     zx_tmp_3d)
1028      CALL histwrite(nid_tra,"jH2O2",itra,zx_tmp_3d,
1029     .                                   iim*(jjm+1)*klev,ndex)
1030
1031      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,hrates(1,1,1),
1032     .     zx_tmp_3d)
1033      CALL histwrite(nid_tra,"wHNO3",itra,zx_tmp_3d,
1034     .                                   iim*(jjm+1)*klev,ndex)
1035
1036      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,krates(1,1,1),
1037     .     zx_tmp_3d)
1038      CALL histwrite(nid_tra,"kN2O5",itra,zx_tmp_3d,
1039     .                                   iim*(jjm+1)*klev,ndex)
1040
1041      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,extflx(1,1,1),
1042     .     zx_tmp_3d)
1043      CALL histwrite(nid_tra,"LghtNO",itra,zx_tmp_3d,
1044     .                                   iim*(jjm+1)*klev,ndex)
1045#endif
1046      DO it=1, grpcnt
1047      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,nas(1,1,it),zx_tmp_3d)
1048      zx_tmp_3d = zx_tmp_3d * dry_mass / nadv_mass(it)
1049      CALL histwrite(nid_tra,grpsym(it),itra,zx_tmp_3d,
1050     .                                   iim*(jjm+1)*klev,ndex)
1051      ENDDO
1052#endif
1053
1054#endif
1055#ifdef INCA_AER
1056C   3d FIELDS
1057
1058      it = id_CIDUSTM
1059       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,scavcoef_st(1,1,it),
1060     .                  zx_tmp_3d)
1061       CALL histwrite(nid_tra,"scavcoef_st",itra,zx_tmp_3d,
1062     .                  iim*(jjm+1)*klev,ndex)
1063
1064       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,scavcoef_cv(1,1,it),
1065     .                  zx_tmp_3d)
1066       CALL histwrite(nid_tra,"scavcoef_cv",itra,zx_tmp_3d,
1067     .                  iim*(jjm+1)*klev,ndex)
1068
1069#endif
1070
1071      DO it=1,nqmax
1072      IF (it.LE.99) THEN
1073       WRITE(str2,'(i2.2)') it
1074
1075#ifdef INCA
1076      IF ( prt_flag_ts(it) == 0 ) CYCLE
1077C champs 2D
1078      CALL gr_fi_ecrit(1, klon,iim,jjm+1, eflux(1,it),zx_tmp_2d)
1079      CALL histwrite(nid_tra,"Emi_"//solsym(it),itra,zx_tmp_2d,
1080     .     iim*(jjm+1),ndex)
1081
1082      CALL gr_fi_ecrit(1, klon,iim,jjm+1, dvel(1,it),zx_tmp_2d)
1083      CALL histwrite(nid_tra,"Dep_"//solsym(it),itra,zx_tmp_2d,
1084     .     iim*(jjm+1),ndex)
1085#ifdef INCA_AER
1086        IF (solsym(it) .eq. 'CIDUSTM'.or.solsym(it) .eq. 'CIN'
1087     .  .or.solsym(it) .eq. 'CSSSM'  .or.solsym(it) .eq. 'CSN'
1088     .  .or.solsym(it) .eq. 'ASSSM'  .or.solsym(it) .eq. 'ASN' ) THEN
1089      CALL gr_fi_ecrit(1, klon,iim,jjm+1,sflux(1,it),zx_tmp_2d)
1090      CALL histwrite(nid_tra,"Sed_"//solsym(it),itra,zx_tmp_2d,
1091     .     iim*(jjm+1),ndex)
1092        endif
1093        IF (solsym(it) .eq. 'CIDUSTM') THEN
1094      CALL gr_fi_ecrit(1, klon,iim,jjm+1,tausum(1),zx_tmp_2d)
1095      CALL histwrite(nid_tra,"OD_"//solsym(it),itra,zx_tmp_2d,
1096     .     iim*(jjm+1),ndex)
1097        endif
1098
1099#endif
1100C champs 3D
1101       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
1102
1103       !Prefer vmr to mmr for transported species
1104       if( adv_mass(it) /= 0. ) then
1105       zx_tmp_3d = zx_tmp_3d * dry_mass / adv_mass(it)
1106       else
1107#ifdef INCA_CH4
1108       if ( solsym(it) == 'OX' ) then
1109       zx_tmp_3d = zx_tmp_3d * dry_mass / nadv_mass(id_o3)
1110       end if
1111#endif
1112       end if
1113
1114       CALL histwrite(nid_tra,solsym(it),itra,zx_tmp_3d,
1115     .                                   iim*(jjm+1)*klev,ndex)
1116#else
1117
1118C champs 3D
1119       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
1120       CALL histwrite(nid_tra,"tr"//str2,itra,zx_tmp_3d,
1121     .                                   iim*(jjm+1)*klev,ndex)
1122       IF (lessivage) THEN
1123       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,flestottr(1,1,it),zx_tmp_3d)
1124       CALL histwrite(nid_tra,"fl"//str2,itra,zx_tmp_3d,
1125     .                                   iim*(jjm+1)*klev,ndex)
1126       ENDIF
1127#endif
1128      ELSE
1129         PRINT*, "Trop de traceurs"
1130         CALL abort
1131      ENDIF
1132      ENDDO
1133
1134#ifdef INCA_CH4
1135      CALL gr_fi_ecrit(1, klon,iim,jjm+1, o3_tr_col(1), zx_tmp_2d)
1136      CALL histwrite(nid_tra,"O3_column",itra,zx_tmp_2d,
1137     .     iim*(jjm+1),ndex)
1138
1139      CALL gr_fi_ecrit(1, klon,iim,jjm+1, co_tr_col(1), zx_tmp_2d)
1140      CALL histwrite(nid_tra,"CO_column",itra,zx_tmp_2d,
1141     .     iim*(jjm+1),ndex)
1142
1143      CALL gr_fi_ecrit(1, klon,iim,jjm+1, ch4_tr_col(1), zx_tmp_2d)
1144      CALL histwrite(nid_tra,"CH4_column",itra,zx_tmp_2d,
1145     .     iim*(jjm+1),ndex)
1146
1147      CALL gr_fi_ecrit(1, klon,iim,jjm+1, no2_tr_col(1), zx_tmp_2d)
1148      CALL histwrite(nid_tra,"NO2_column",itra,zx_tmp_2d,
1149     .     iim*(jjm+1),ndex)
1150
1151      CALL gr_fi_ecrit(1, klon,iim,jjm+1, o3_st_flx(1), zx_tmp_2d)
1152      CALL histwrite(nid_tra,"O3_ste",itra,zx_tmp_2d,
1153     .     iim*(jjm+1),ndex)
1154
1155      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,o3_prod(1,1),
1156     .     zx_tmp_3d)
1157      CALL histwrite(nid_tra,"O3_prod",itra,zx_tmp_3d,
1158     .                                   iim*(jjm+1)*klev,ndex)
1159
1160      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,o3_loss(1,1),
1161     .     zx_tmp_3d)
1162      CALL histwrite(nid_tra,"O3_loss",itra,zx_tmp_3d,
1163     .                                   iim*(jjm+1)*klev,ndex)
1164
1165!     ... Special section for daytime averaging
1166!       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,day_cnt(1,1),
1167!    .       zx_tmp_3d)
1168!       CALL histwrite(nid_tra,"day_cnt",itra,zx_tmp_3d,
1169!    .                                  iim*(jjm+1)*klev,ndex)
1170!       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,no_daytime(1,1),
1171!    .       zx_tmp_3d)
1172!       CALL histwrite(nid_tra,"NO_day",itra,zx_tmp_3d,
1173!    .                                  iim*(jjm+1)*klev,ndex)
1174
1175#endif
1176
1177      if (ok_sync) call histsync(nid_tra)
1178
1179      if (lafin) then
1180!DH      print*, 'c est la fin de la physique'
1181         open (99,file='restarttrac',  form='formatted')
1182         do i=1,klon
1183             write(99,*) trs(i,1)
1184         enddo
1185         PRINT*, 'Ecriture du fichier restarttrac'
1186         close(99)
1187      else
1188!DH      print*, 'physique pas fini'
1189      endif
1190
1191
1192      RETURN
1193      END
Note: See TracBrowser for help on using the repository browser.