Context Navigation

phytrac.F @ 345

Last change on this file since 345 was 344, checked in by lmdz, 23 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

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
	212	INTEGER nid_tra
	213	SAVE nid_tra
[204]	214	c REAL x(klon,klev,nbtr+2) ! traceurs
	215	INTEGER ndex(1)
[2]	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
	220	C
	221	C Variables liees a l'ecriture de la bande histoire : phytrac.nc
	222	c
	223	INTEGER ecrit_tra
	224	SAVE ecrit_tra
[53]	225	logical ok_sync
	226	parameter (ok_sync = .true.)
[2]	227	C
	228	C nature du traceur
	229	c
	230	logical aerosol(nbtr) ! Nature du traceur
	231	c ! aerosol(it) = true => aerosol
	232	c ! aerosol(it) = false => gaz
	233	c ! 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
	237	C
	238	c======================================================================
	239	c
	240	c Declaration des procedures appelees
	241	c
	242	c--modif convection tiedtke
[204]	243	INTEGER i, k, it,itap
	244	save itap
[2]	245	REAL delp(klon,klev)
	246	c--end modif
	247	c
	248	c Variables liees a l'ecriture de la bande histoire physique
	249	c
	250	c Variables locales pour effectuer les appels en serie
	251	c----------------------------------------------------
	252	c
	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
	257	c ! radioactive du rn - > pb
	258	REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
	259	c ! par impaction
	260	REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
	261	c ! 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
	269	c ! dans chaque couche
	270
	271	C
	272	character*20 modname
	273	character*80 abort_message
	274	c
	275	c Controles
	276	c-------------
	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./
[344]	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
[2]	306	c
	307	c======================================================================
	308	c ecrit_tra = NINT(86400./pdtphys *1.0) ! tous les jours
	309	modname='phytrac'
[344]	310	ecrit_tra = NINT(86400./pdtphys *ecritphy)
	311	!DH print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
	312	ps(:) = paprs(:,1)
[2]	313
[53]	314	c print*,'DANS PHYTRAC debutphy=',debutphy
[2]	315
	316	if (debutphy) then
	317
[204]	318	print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
	319	ecrit_tra = NINT(86400./pdtphys/2.) ! tous les 12H
	320	c ecrit_tra = NINT(86400./pdtphys) ! tous les 24H
	321
[2]	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
[344]	330	!DH PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra
	331
[2]	332	C
[53]	333	CALL ymds2ju(anne_ini, 1, 1, 0.0, zjulian)
[344]	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
[2]	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
[344]	345	DO ll=1,klev
	346	znivsig(ll)=float(ll)
	347	ENDDO
[2]	348	CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlat,zx_lat)
	349	CALL histbeg("histrac", iim,zx_lon, jjm+1,zx_lat,
[344]	350	. 1,iim,1,jjm+1, itra, zjulian, pdtphys,
[2]	351	. nhori, nid_tra)
[344]	352
	353	call histvert(nid_tra, "presnivs", "presnivs", "mb",
[204]	354	. klev, presnivs, nvert)
[344]	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
[204]	362	zsto = pdtphys
	363	zout = pdtphys * FLOAT(ecrit_tra)
[2]	364	c
	365	CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
	366	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	367	. "once", zsto,zout)
	368	c
	369	CALL histdef(nid_tra, "aire", "Grid area", "-",
	370	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	371	. "once", zsto,zout)
[204]	372
[344]	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)
[204]	377
[344]	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",
[204]	383	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
[344]	384	. "ave(X)", zsto,zout)
	385
	386	CALL histdef(nid_tra, "u", "zonal wind component", "m/s",
[204]	387	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
[344]	388	. "ave(X)", zsto,zout)
	389
	390	CALL histdef(nid_tra, "v", "zonal wind component", "m/s",
[204]	391	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	392	. "ave(X)", zsto,zout)
[344]	393
	394	CALL histdef(nid_tra, "h2o", "Specific Humidity", "MMR",
[204]	395	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	396	. "ave(X)", zsto,zout)
[344]	397
	398	CALL histdef(nid_tra, "pmid", "Pressure", "Pa",
[204]	399	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	400	. "ave(X)", zsto,zout)
[344]	401
	402	CALL histdef(nid_tra, "pdel", "Delta Pressure", "Pa",
[204]	403	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	404	. "ave(X)", zsto,zout)
[344]	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",
[204]	410	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	411	. "ave(X)", zsto,zout)
[344]	412	ENDDO
	413	DO it=1, hetcnt
	414	WRITE(str2,'(i2.2)') it
	415	CALL histdef(nid_tra, "w"//str2,"w"//str2, "S-1",
[204]	416	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	417	. "ave(X)", zsto,zout)
[344]	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",
[204]	423	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	424	. "ave(X)", zsto,zout)
[344]	425	ENDDO
[204]	426
[344]	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
	464	C 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
[2]	474	DO it=1,nqmax
	475	IF (it.LE.99) THEN
	476	WRITE(str2,'(i2.2)') it
[344]	477	C 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
[2]	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
[344]	515	#endif
[2]	516	ELSE
	517	PRINT*, "Trop de traceurs"
	518	CALL abort
	519	ENDIF
	520	ENDDO
[344]	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
[2]	555	CALL histend(nid_tra)
[204]	556	ndex(1) = 0
	557	c
	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)
	561	C
	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)
[53]	565
[2]	566	c======================================================================
	567	c Initialisation de certaines variables pour le Rn et le Pb
	568	c======================================================================
	569
	570	c Initialisation du traceur dans le sol (couche limite radonique)
	571	c
[53]	572	c print*,'valeur de debut dans phytrac :',debutphy
[2]	573	do it=1,nqmax
	574	do i=1,klon
	575	trs(i,it) = 0.
	576	enddo
	577	END DO
[199]	578
	579	open (99,file='starttrac',status='old',
	580	. err=999,form='formatted')
	581	read(99,*) (trs(i,1),i=1,klon)
	582	999 close(99)
[344]	583	!DH print*, 'apres starttrac'
[199]	584
[2]	585	c Initialisation de la fraction d'aerosols lessivee
	586	c
	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
	595	c
	596	c Initialisation de la nature des traceurs
	597	c
	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
	603	c
	604	ENDIF ! fin debutphy
	605	c Initialisation du traceur dans le sol (couche limite radonique)
	606	if(inirnpb) THEN
	607	c
	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
	613	c
	614	call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
	615	. ,vdeptr,scavtr)
	616	inirnpb=.false.
	617	endif
[204]	618
[344]	619	#ifdef INCA
	620	!======================================================================
	621	! Chimie
	622	!======================================================================
[204]	623
[344]	624	calday = FLOAT(julien) + gmtime
	625	ncsec = NINT (86400.*gmtime)
[204]	626
[344]	627	DO k = 1, nlev
	628	pdel(:,k) = paprs(:,k) - paprs (:,k+1)
	629	END DO
[204]	630
[344]	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
[204]	638
	639
[344]	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
[204]	655
[344]	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
	699	C
[204]	700
[2]	701	c======================================================================
	702	c Calcul de l'effet de la convection
	703	c======================================================================
[344]	704	!DH print*,'Avant convection'
[204]	705	do it=1,nqmax
	706	WRITE(itn,'(i1)') it
	707	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
	708	enddo
[2]	709
	710	if (convection) then
	711
[344]	712	!DH print*,'Pas de temps dans phytrac : ',pdtphys
[2]	713	DO it=1, nqmax
[344]	714	#ifdef INCA
	715	IF ( conv_flg(it) == 0 ) CYCLE
	716	#endif
[2]	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
[344]	724	c 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
[2]	730	ENDDO
[344]	731	c print*,'apres nflxtr'
[2]	732
	733	endif ! convection
[204]	734	c print*,'Apres convection'
	735	c do it=1,nqmax
	736	c WRITE(itn,'(i1)') it
	737	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
	738	c enddo
[2]	739
	740	c======================================================================
	741	c Calcul de l'effet de la couche limite
	742	c======================================================================
[204]	743	c print *,'Avant couchelimite'
	744	c do it=1,nqmax
	745	c WRITE(itn,'(i1)') it
	746	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
	747	c enddo
[2]	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
	757	C maf modif pour tenir compte du cas rnpb + traceur
	758	DO it=1, nqmax
[344]	759	#ifdef INCA
	760	IF ( pbl_flg(it) == 0 ) CYCLE
	761	#endif
	762	C print *,'it',it,clsol(it)
[2]	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
	777	c
	778	c Traceur ds sol
	779	c
	780	DO i = 1, klon
	781	trs(i,it) = trs(i,it) + d_trs(i)
	782	END DO
	783	C
	784	C maf provisoire suppression des prints
	785	C WRITE(itn,'(i1)') it
	786	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
	787	else ! couche limite avec flux prescrit
[344]	788
	789	#ifdef INCA
	790	DO k = 1, klon
	791	source(k) = eflux(k,it)-dflux(k,it)
	792	END DO
	793	#else
[2]	794	Cmaf provisoire source / traceur a creer
	795	DO i=1, klon
	796	source(i) = 0.0 ! pas de source, pour l'instant
	797	ENDDO
[344]	798	#endif
[2]	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
	808	Cmaf provisoire suppression des prints
	809	Cmaf WRITE(itn,'(i1)') it
	810	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
	811	endif
	812	ENDDO
	813	c
	814	endif ! couche limite
	815
[204]	816	c print*,'Apres couchelimite'
	817	c do it=1,nqmax
	818	c WRITE(itn,'(i1)') it
	819	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
	820	c enddo
[2]	821
	822	c======================================================================
	823	c Calcul de l'effet du puits radioactif
	824	c======================================================================
	825
	826	C MAF il faudrait faire une modification pour passer dans radiornpb
	827	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
	828	if(rnpb) then
[344]	829	print *, 'decroissance radiactive activee'
[2]	830	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
	831	C
	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
	843	c
	844	endif ! rnpb decroissance radioactive
[344]	845
[2]	846	c======================================================================
	847	c Calcul de l'effet de la precipitation
	848	c======================================================================
	849
[344]	850	!DH print*,'LESSIVAGE =',lessivage
[2]	851	IF (lessivage) THEN
	852
[53]	853	c print*,'avant lessivage'
[2]	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
	864	c
	865	c tendance des aerosols nuclees et impactes
	866	c
	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
	879	c
	880	c Mises a jour des traceurs + calcul des flux de lessivage
	881	c Mise a jour due a l'impaction et a la nucleation
	882	c
[204]	883	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
	884	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
	885	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
[2]	886	DO it = 1, nqmax
[204]	887	c print*,'IT=',it,aerosol(it)
[2]	888	IF (aerosol(it)) THEN
[204]	889	c print*,'IT=',it,' On lessive'
[2]	890	DO k = 1, nlev
	891	DO i = 1, klon
[204]	892	tr_seri(i,k,it)=tr_seri(i,k,it)
	893	s frac_impa(i,k)frac_nucl(i,k)
[2]	894	ENDDO
	895	ENDDO
	896	ENDIF
	897	ENDDO
[204]	898	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
[2]	899	c
	900	c Flux lessivage total
	901	c
	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
	912	c
	913	Cmaf suppression provisoire
	914	Cmaf WRITE(itn,'(i1)') it
	915	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
	916	ENDDO
	917	c
[53]	918	c print*,'apres lessivage'
[2]	919	ENDIF
	920	Cc
	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
[344]	931
	932
	933
[2]	934	itra=itra+1
[204]	935	ndex(1) = 0
[344]	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
	1054	C 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
[2]	1069	DO it=1,nqmax
	1070	IF (it.LE.99) THEN
	1071	WRITE(str2,'(i2.2)') it
[344]	1072
	1073	#ifdef INCA
	1074	IF ( prt_flag_ts(it) == 0 ) CYCLE
	1075	C 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
[2]	1098	C champs 3D
	1099	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
[344]	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
	1116	C champs 3D
	1117	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
[2]	1118	CALL histwrite(nid_tra,"tr"//str2,itra,zx_tmp_3d,
[204]	1119	. iim(jjm+1)klev,ndex)
[344]	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
[2]	1126	ELSE
	1127	PRINT*, "Trop de traceurs"
	1128	CALL abort
	1129	ENDIF
	1130	ENDDO
	1131
[344]	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)
[204]	1136
[344]	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)
[204]	1140
[344]	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)
[204]	1144
[344]	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)
[204]	1148
[344]	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
[199]	1177	if (lafin) then
[344]	1178	!DH print*, 'c est la fin de la physique'
[199]	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
[344]	1186	!DH print*, 'physique pas fini'
[199]	1187	endif
	1188
	1189
[2]	1190	RETURN
	1191	END

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

Download in other formats:

Original Format