Context Navigation

phytrac.F @ 5415

Last change on this file since 5415 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

Rev	Line
[199]	1	c
[207]	2	c $Header$
[199]	3	c
[2]	4	SUBROUTINE phytrac (rnpb,
[344]	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)
[2]	52	USE ioipsl
[344]	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
[2]	68	IMPLICIT none
	69	c======================================================================
	70	c Auteur(s) FH
	71	c Objet: Moniteur general des tendances traceurs
	72	c
	73	cAA Remarques en vrac:
	74	cAA--------------------
	75	cAA 1/ le call phytrac se fait avec nqmax-2 donc nous avons bien
	76	cAA les vrais traceurs (nbtr) dans phytrac (pas la vapeur ni eau liquide)
	77	cAA 2/ Le choix du radon et du pb se fait juste avec un data
	78	cAA (peu propre). Peut-etre pourrait-on prevoir dans l'avenir
	79	cAA une variable "type de traceur"
	80	c======================================================================
	81	#include "YOMCST.h"
	82	#include "dimensions.h"
[344]	83	#include "paramet.h"
[2]	84	#include "dimphy.h"
	85	#include "indicesol.h"
[344]	86	#include "comvert.h"
[209]	87	#include "temps.h"
	88	#include "control.h"
[2]	89	c======================================================================
	90
	91	c Arguments:
	92	c
	93	c EN ENTREE:
	94	c ==========
	95	c
	96	c divers:
	97	c -------
	98	c
[344]	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
[2]	107	real pdtphys ! pas d'integration pour la physique (seconde)
	108	real t_seri(nlon,nlev) ! temperature
[344]	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
[2]	118	real tr_seri(nlon,nlev,nbtr) ! traceur
	119	real paprs(nlon,nlev+1) ! pression pour chaque inter-couche (en Pa)
[344]	120	real ps(nlon) ! pression surface
[2]	121	real pplay(nlon,nlev) ! pression pour le mileu de chaque couche (en Pa)
[344]	122	real pphi(nlon,klev) ! geopotentiel
	123	real znivsig(klev) ! niveaux sigma
[2]	124	real paire(klon)
	125	real pphis(klon)
[344]	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)
[2]	129	logical debutphy ! le flag de l'initialisation de la physique
[199]	130	logical lafin ! le flag de la fin de la physique
	131
[2]	132	integer ll
	133	c
	134	cAA Rem : nbtr : nombre de vrais traceurs est defini dans dimphy.h
	135	c
	136	c convection:
	137	c -----------
	138	c
	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
	145	c
	146	c Couche limite:
	147	c --------------
	148	c
	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
	155	c
	156	c Lessivage:
	157	c ----------
	158	c
	159	c pour le ON-LINE
	160	c
	161	REAL frac_impa(nlon,nlev) ! fraction d'aerosols impactes
	162	REAL frac_nucl(nlon,nlev) ! fraction d'aerosols nuclees
	163	c
	164	cAA
	165	cAA Arguments necessaires pour les sources et puits de traceur:
	166	cAA ----------------
	167	cAA
	168	real ftsol(nlon,nbsrf) ! Temperature du sol (surf)(Kelvin)
	169	real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol)
[204]	170	c abder
	171	real pftsol1(nlon),pftsol2(nlon),pftsol3(nlon),pftsol4(nlon)
	172	real ppsrf1(nlon),ppsrf2(nlon),ppsrf3(nlon),ppsrf4(nlon)
	173	c fin
[2]	174	cAA ----------------------------
	175	cAA VARIABLES LOCALES TRACEURS
	176	cAA ----------------------------
	177	cAA
	178	cAA Sources et puits des traceurs:
	179	cAA ------------------------------
	180	cAA
	181	cAA 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
	185	cAA
	186	cAA Pour la source de radon et son reservoir de sol
	187	cAA ................................................
	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
	193	c Masque de l'echange avec la surface
	194	c (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
	206	cAA
[344]	207	CHARACTER*2 itn
[2]	208	C maf ioipsl
	209	CHARACTER*2 str2
[344]	210	INTEGER nhori, nvert, nverta, nvertb, nverts
[2]	211	REAL zsto, zout, zjulian
[348]	212	integer idayref
[2]	213	INTEGER nid_tra
	214	SAVE nid_tra
[204]	215	c REAL x(klon,klev,nbtr+2) ! traceurs
	216	INTEGER ndex(1)
[2]	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
	221	C
	222	C Variables liees a l'ecriture de la bande histoire : phytrac.nc
	223	c
	224	INTEGER ecrit_tra
	225	SAVE ecrit_tra
[53]	226	logical ok_sync
	227	parameter (ok_sync = .true.)
[2]	228	C
	229	C nature du traceur
	230	c
	231	logical aerosol(nbtr) ! Nature du traceur
	232	c ! aerosol(it) = true => aerosol
	233	c ! aerosol(it) = false => gaz
	234	c ! 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
	238	C
	239	c======================================================================
	240	c
	241	c Declaration des procedures appelees
	242	c
	243	c--modif convection tiedtke
[204]	244	INTEGER i, k, it,itap
	245	save itap
[2]	246	REAL delp(klon,klev)
	247	c--end modif
	248	c
	249	c Variables liees a l'ecriture de la bande histoire physique
	250	c
	251	c Variables locales pour effectuer les appels en serie
	252	c----------------------------------------------------
	253	c
	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
	258	c ! radioactive du rn - > pb
	259	REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
	260	c ! par impaction
	261	REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
	262	c ! 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
	270	c ! dans chaque couche
	271
	272	C
	273	character*20 modname
	274	character*80 abort_message
	275	c
	276	c Controles
	277	c-------------
	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./
[344]	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
[2]	307	c
	308	c======================================================================
	309	c ecrit_tra = NINT(86400./pdtphys *1.0) ! tous les jours
	310	modname='phytrac'
[344]	311	ecrit_tra = NINT(86400./pdtphys *ecritphy)
	312	!DH print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
	313	ps(:) = paprs(:,1)
[2]	314
[53]	315	c print*,'DANS PHYTRAC debutphy=',debutphy
[2]	316
	317	if (debutphy) then
	318
[204]	319	print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
	320	ecrit_tra = NINT(86400./pdtphys/2.) ! tous les 12H
	321	c ecrit_tra = NINT(86400./pdtphys) ! tous les 24H
	322
[2]	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
[344]	331	!DH PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra
	332
[2]	333	C
[348]	334	idayref = 1
	335	CALL ymds2ju(anne_ini, 1, idayref, 0.0, zjulian)
[344]	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
[2]	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
[344]	347	DO ll=1,klev
	348	znivsig(ll)=float(ll)
	349	ENDDO
[2]	350	CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlat,zx_lat)
	351	CALL histbeg("histrac", iim,zx_lon, jjm+1,zx_lat,
[344]	352	. 1,iim,1,jjm+1, itra, zjulian, pdtphys,
[2]	353	. nhori, nid_tra)
[344]	354
	355	call histvert(nid_tra, "presnivs", "presnivs", "mb",
[204]	356	. klev, presnivs, nvert)
[344]	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
[204]	364	zsto = pdtphys
	365	zout = pdtphys * FLOAT(ecrit_tra)
[2]	366	c
	367	CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
	368	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	369	. "once", zsto,zout)
	370	c
	371	CALL histdef(nid_tra, "aire", "Grid area", "-",
	372	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	373	. "once", zsto,zout)
[204]	374
[344]	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)
[204]	379
[344]	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",
[204]	385	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
[344]	386	. "ave(X)", zsto,zout)
	387
	388	CALL histdef(nid_tra, "u", "zonal wind component", "m/s",
[204]	389	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
[344]	390	. "ave(X)", zsto,zout)
	391
	392	CALL histdef(nid_tra, "v", "zonal wind component", "m/s",
[204]	393	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	394	. "ave(X)", zsto,zout)
[344]	395
	396	CALL histdef(nid_tra, "h2o", "Specific Humidity", "MMR",
[204]	397	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	398	. "ave(X)", zsto,zout)
[344]	399
	400	CALL histdef(nid_tra, "pmid", "Pressure", "Pa",
[204]	401	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	402	. "ave(X)", zsto,zout)
[344]	403
	404	CALL histdef(nid_tra, "pdel", "Delta Pressure", "Pa",
[204]	405	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	406	. "ave(X)", zsto,zout)
[344]	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",
[204]	412	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	413	. "ave(X)", zsto,zout)
[344]	414	ENDDO
	415	DO it=1, hetcnt
	416	WRITE(str2,'(i2.2)') it
	417	CALL histdef(nid_tra, "w"//str2,"w"//str2, "S-1",
[204]	418	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	419	. "ave(X)", zsto,zout)
[344]	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",
[204]	425	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	426	. "ave(X)", zsto,zout)
[344]	427	ENDDO
[204]	428
[344]	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
	466	C 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
[2]	476	DO it=1,nqmax
	477	IF (it.LE.99) THEN
	478	WRITE(str2,'(i2.2)') it
[344]	479	C 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
[2]	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
[344]	517	#endif
[2]	518	ELSE
	519	PRINT*, "Trop de traceurs"
	520	CALL abort
	521	ENDIF
	522	ENDDO
[344]	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
[2]	557	CALL histend(nid_tra)
[204]	558	ndex(1) = 0
	559	c
	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)
	563	C
	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)
[53]	567
[2]	568	c======================================================================
	569	c Initialisation de certaines variables pour le Rn et le Pb
	570	c======================================================================
	571
	572	c Initialisation du traceur dans le sol (couche limite radonique)
	573	c
[53]	574	c print*,'valeur de debut dans phytrac :',debutphy
[2]	575	do it=1,nqmax
	576	do i=1,klon
	577	trs(i,it) = 0.
	578	enddo
	579	END DO
[199]	580
	581	open (99,file='starttrac',status='old',
	582	. err=999,form='formatted')
	583	read(99,*) (trs(i,1),i=1,klon)
	584	999 close(99)
[344]	585	!DH print*, 'apres starttrac'
[199]	586
[2]	587	c Initialisation de la fraction d'aerosols lessivee
	588	c
	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
	597	c
	598	c Initialisation de la nature des traceurs
	599	c
	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
	605	c
	606	ENDIF ! fin debutphy
	607	c Initialisation du traceur dans le sol (couche limite radonique)
	608	if(inirnpb) THEN
	609	c
	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
	615	c
	616	call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
	617	. ,vdeptr,scavtr)
	618	inirnpb=.false.
	619	endif
[204]	620
[344]	621	#ifdef INCA
	622	!======================================================================
	623	! Chimie
	624	!======================================================================
[204]	625
[344]	626	calday = FLOAT(julien) + gmtime
	627	ncsec = NINT (86400.*gmtime)
[204]	628
[344]	629	DO k = 1, nlev
	630	pdel(:,k) = paprs(:,k) - paprs (:,k+1)
	631	END DO
[204]	632
[344]	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
[204]	640
	641
[344]	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
[204]	657
[344]	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
	701	C
[204]	702
[2]	703	c======================================================================
	704	c Calcul de l'effet de la convection
	705	c======================================================================
[344]	706	!DH print*,'Avant convection'
[204]	707	do it=1,nqmax
	708	WRITE(itn,'(i1)') it
	709	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
	710	enddo
[2]	711
	712	if (convection) then
	713
[344]	714	!DH print*,'Pas de temps dans phytrac : ',pdtphys
[2]	715	DO it=1, nqmax
[344]	716	#ifdef INCA
	717	IF ( conv_flg(it) == 0 ) CYCLE
	718	#endif
[2]	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
[344]	726	c 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
[2]	732	ENDDO
[344]	733	c print*,'apres nflxtr'
[2]	734
	735	endif ! convection
[204]	736	c print*,'Apres convection'
	737	c do it=1,nqmax
	738	c WRITE(itn,'(i1)') it
	739	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
	740	c enddo
[2]	741
	742	c======================================================================
	743	c Calcul de l'effet de la couche limite
	744	c======================================================================
[204]	745	c print *,'Avant couchelimite'
	746	c do it=1,nqmax
	747	c WRITE(itn,'(i1)') it
	748	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
	749	c enddo
[2]	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
	759	C maf modif pour tenir compte du cas rnpb + traceur
	760	DO it=1, nqmax
[344]	761	#ifdef INCA
	762	IF ( pbl_flg(it) == 0 ) CYCLE
	763	#endif
	764	C print *,'it',it,clsol(it)
[2]	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
	779	c
	780	c Traceur ds sol
	781	c
	782	DO i = 1, klon
	783	trs(i,it) = trs(i,it) + d_trs(i)
	784	END DO
	785	C
	786	C maf provisoire suppression des prints
	787	C WRITE(itn,'(i1)') it
	788	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
	789	else ! couche limite avec flux prescrit
[344]	790
	791	#ifdef INCA
	792	DO k = 1, klon
	793	source(k) = eflux(k,it)-dflux(k,it)
	794	END DO
	795	#else
[2]	796	Cmaf provisoire source / traceur a creer
	797	DO i=1, klon
	798	source(i) = 0.0 ! pas de source, pour l'instant
	799	ENDDO
[344]	800	#endif
[2]	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
	810	Cmaf provisoire suppression des prints
	811	Cmaf WRITE(itn,'(i1)') it
	812	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
	813	endif
	814	ENDDO
	815	c
	816	endif ! couche limite
	817
[204]	818	c print*,'Apres couchelimite'
	819	c do it=1,nqmax
	820	c WRITE(itn,'(i1)') it
	821	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
	822	c enddo
[2]	823
	824	c======================================================================
	825	c Calcul de l'effet du puits radioactif
	826	c======================================================================
	827
	828	C MAF il faudrait faire une modification pour passer dans radiornpb
	829	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
	830	if(rnpb) then
[344]	831	print *, 'decroissance radiactive activee'
[2]	832	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
	833	C
	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
	845	c
	846	endif ! rnpb decroissance radioactive
[344]	847
[2]	848	c======================================================================
	849	c Calcul de l'effet de la precipitation
	850	c======================================================================
	851
[344]	852	!DH print*,'LESSIVAGE =',lessivage
[2]	853	IF (lessivage) THEN
	854
[53]	855	c print*,'avant lessivage'
[2]	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
	866	c
	867	c tendance des aerosols nuclees et impactes
	868	c
	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
	881	c
	882	c Mises a jour des traceurs + calcul des flux de lessivage
	883	c Mise a jour due a l'impaction et a la nucleation
	884	c
[204]	885	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
	886	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
	887	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
[2]	888	DO it = 1, nqmax
[204]	889	c print*,'IT=',it,aerosol(it)
[2]	890	IF (aerosol(it)) THEN
[204]	891	c print*,'IT=',it,' On lessive'
[2]	892	DO k = 1, nlev
	893	DO i = 1, klon
[204]	894	tr_seri(i,k,it)=tr_seri(i,k,it)
	895	s frac_impa(i,k)frac_nucl(i,k)
[2]	896	ENDDO
	897	ENDDO
	898	ENDIF
	899	ENDDO
[204]	900	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
[2]	901	c
	902	c Flux lessivage total
	903	c
	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
	914	c
	915	Cmaf suppression provisoire
	916	Cmaf WRITE(itn,'(i1)') it
	917	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
	918	ENDDO
	919	c
[53]	920	c print*,'apres lessivage'
[2]	921	ENDIF
	922	Cc
	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
[344]	933
	934
	935
[2]	936	itra=itra+1
[204]	937	ndex(1) = 0
[344]	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
	1056	C 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
[2]	1071	DO it=1,nqmax
	1072	IF (it.LE.99) THEN
	1073	WRITE(str2,'(i2.2)') it
[344]	1074
	1075	#ifdef INCA
	1076	IF ( prt_flag_ts(it) == 0 ) CYCLE
	1077	C 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
[2]	1100	C champs 3D
	1101	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
[344]	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
	1118	C champs 3D
	1119	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
[2]	1120	CALL histwrite(nid_tra,"tr"//str2,itra,zx_tmp_3d,
[204]	1121	. iim(jjm+1)klev,ndex)
[344]	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
[2]	1128	ELSE
	1129	PRINT*, "Trop de traceurs"
	1130	CALL abort
	1131	ENDIF
	1132	ENDDO
	1133
[344]	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)
[204]	1138
[344]	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)
[204]	1142
[344]	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)
[204]	1146
[344]	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)
[204]	1150
[344]	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
[199]	1179	if (lafin) then
[344]	1180	!DH print*, 'c est la fin de la physique'
[199]	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
[344]	1188	!DH print*, 'physique pas fini'
[199]	1189	endif
	1190
	1191
[2]	1192	RETURN
	1193	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format