Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

physiq.F_bugT @ 3558

Last change on this file since 3558 was 3175, checked in by emillour, 19 months ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
File size: 83.0 KB

Rev	Line
[3175]	1	subroutine physiq(ngrid,nlayer,nq,
	2	* firstcall,lastcall,
	3	* pday,ptime,ptimestep,
	4	* pplev,pplay,pphi,
	5	* pu,pv,pt,pq,
	6	* pw,
	7	* pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
	8
	9	use radinc_h, only : naerkind
	10	use planet_h
	11	USE ioipsl_getincom
	12	use comgeomfi_h
	13	use comsoil_h
	14	use aerosol_mod
	15	implicit none
	16
	17	!==================================================================
	18	! Purpose
	19	! -------
	20	! Central subroutine for all the physics parameterisations in the
	21	! universal model. Originally adapted from the Mars LMDZ model.
	22	!
	23	! The model can be run without or with tracer transport
	24	! depending on the value of "tracer" in file "callphys.def".
	25	!
	26	! It includes:
	27	!
	28	! 1. Initialization:
	29	! 1.1 Firstcall initializations
	30	! 1.2 Initialization for every call to physiq
	31	! 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
	32	! 2. Compute radiative transfer tendencies
	33	! (longwave and shortwave).
	34	! 4. Vertical diffusion (turbulent mixing):
	35	! 5. Convective adjustment
	36	! 6. Condensation and sublimation of gases (currently just CO2).
	37	! 7. TRACERS :
	38	! 7c. other schemes for tracer transport (lifting, sedimentation)
	39	! 7d. updates (pressure variations, surface budget)
	40	! 9. Surface and sub-surface temperature calculations
	41	! 10. Write outputs :
	42	! - "startfi", "histfi" if it's time
	43	! - Saving statistics if "callstats = .true."
	44	! - Output any needed variables in "diagfi"
	45	! 10. Diagnostics: mass conservation of tracers, radiative energy balance etc.
	46	!
	47	! arguments
	48	! ---------
	49	!
	50	! input
	51	! -----
	52	! ecri period (in dynamical timestep) to write output
	53	! ngrid Size of the horizontal grid.
	54	! All internal loops are performed on that grid.
	55	! nlayer Number of vertical layers.
	56	! nq Number of advected fields
	57	! firstcall True at the first call
	58	! lastcall True at the last call
	59	! pday Number of days counted from the North. Spring
	60	! equinoxe.
	61	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
	62	! ptimestep timestep (s)
	63	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
	64	! pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
	65	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
	66	! pu(ngrid,nlayer) u component of the wind (ms-1)
	67	! pv(ngrid,nlayer) v component of the wind (ms-1)
	68	! pt(ngrid,nlayer) Temperature (K)
	69	! pq(ngrid,nlayer,nq) Advected fields
	70	! pudyn(ngrid,nlayer) \
	71	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
	72	! ptdyn(ngrid,nlayer) / corresponding variables
	73	! pqdyn(ngrid,nlayer,nq) /
	74	! pw(ngrid,?) vertical velocity
	75	!
	76	! output
	77	! ------
	78	!
	79	! pdu(ngrid,nlayermx) \
	80	! pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding
	81	! pdt(ngrid,nlayermx) / variables due to physical processes.
	82	! pdq(ngrid,nlayermx) /
	83	! pdpsrf(ngrid) /
	84	! tracerdyn call tracer in dynamical part of GCM ?
	85	!
	86	! Authors
	87	! -------
	88	! Frederic Hourdin 15/10/93
	89	! Francois Forget 1994
	90	! Christophe Hourdin 02/1997
	91	! Subroutine completly rewritten by F.Forget (01/2000)
	92	! Water ice clouds: Franck Montmessin (update 06/2003)
	93	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
	94	! New correlated-k radiative scheme: R. Wordsworth (2009)
	95	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
	96	! Pluto version improved :T. Bertrand (2015,2020)
	97	!==================================================================
	98
	99	c 0. Declarations :
	100	c ------------------
	101	#include "dimensions.h"
	102	#include "dimphys.h"
	103	#include "surfdat.h"
	104	#include "comdiurn.h"
	105	#include "callkeys.h"
	106	#include "comcstfi.h"
	107	#include "comsaison.h"
	108	#include "control.h"
	109	#include "comg1d.h"
	110	#include "tracer.h"
	111	#include "netcdf.inc"
	112
	113	c Arguments :
	114	c -----------
	115
	116	c inputs:
	117	c -------
	118	INTEGER ngrid,nlayer,nq
	119	REAL ptimestep
	120	REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
	121	REAL pphi(ngridmx,nlayer)
	122	REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
	123	REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
	124	REAL pw(ngridmx,nlayer) ! pvervel transmitted by dyn3d
	125	REAL zh(ngridmx,nlayermx) ! potential temperature (K)
	126	LOGICAL firstcall,lastcall
	127	REAL pday
	128	REAL ptime
	129	logical tracerdyn
	130
	131	c outputs:
	132	c --------
	133	c physical tendencies
	134	REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
	135	REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
	136	REAL pdpsrf(ngridmx) ! surface pressure tendency
	137	REAL picen2(ngrid)
	138
	139	c Local saved variables:
	140	c ----------------------
	141	INTEGER day_ini ! Initial date of the run (sol since Ls=0)
	142	INTEGER icount ! counter of calls to physiq during the run.
	143	REAL tsurf(ngridmx) ! Surface temperature (K)
	144	REAL tsoil(ngridmx,nsoilmx) ! sub-surface temperatures (K)
	145	REAL tidat(ngridmx,nsoilmx) ! thermal inertia soil
	146	REAL tidat_out(ngridmx,nlayermx) ! thermal inertia soil but output on vertical grid
	147	REAL tsub(ngridmx) ! temperature of the deepest layer
	148	REAL albedo(ngridmx) ! Surface albedo
	149
	150	REAL emis(ngridmx) ! Thermal IR surface emissivity
	151	REAL dtrad(ngridmx,nlayermx) ! Net atm. radiative heating rate (K.s-1)
	152	REAL fluxrad_sky(ngridmx) ! rad. flux from sky absorbed by surface (W.m-2)
	153	REAL fluxrad(ngridmx) ! Net radiative surface flux (W.m-2)
	154	REAL capcal(ngridmx) ! surface heat capacity (J m-2 K-1)
	155	REAL dplanck(ngridmx) ! for output (W.s/m2/K)
	156	REAL fluxgrd(ngridmx) ! surface conduction flux (W.m-2)
	157	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
	158	REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy
	159	REAL qsurf1(ngridmx,nqmx) ! saving qsurf to calculate flux over long timescales kg.m-2
	160	REAL flusurf(ngridmx,nqmx) ! flux cond/sub kg.m-2.s-1
	161	REAL flusurfold(ngridmx,nqmx) ! old flux cond/sub kg.m-2.s-1
	162	REAL totmass ! total mass of n2 (atm + surf)
	163	REAL globave ! globalaverage 2D ps
	164	REAL globaveice(nqmx) ! globalaverage 2D ice
	165	REAL globavenewice(nqmx) ! globalaverage 2D ice
	166
	167	REAL ptime0 ! store the first time
	168	SAVE ptime0
	169
	170	REAL dstep
	171	REAL glastep ! step en pluto day pour etaler le glacier
	172	SAVE glastep
	173	DATA glastep/20/
	174
	175	SAVE day_ini,icount
	176	SAVE tsurf,tsoil,tsub
	177	SAVE albedo,emis,q2
	178	SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
	179	SAVE flusurfold
	180
	181	REAL stephan
	182	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
	183	SAVE stephan
	184
	185	REAL acond,bcond
	186	SAVE acond,bcond
	187	REAL tcond1p4Pa
	188	DATA tcond1p4Pa/38/
	189
	190	c Local variables :
	191	c -----------------
	192	! Tendencies for the paleoclimate mode
	193	REAL,SAVE :: qsurfyear(ngridmx,nqmx) ! kg.m-2 averaged mass of ice lost/gained in the last Pluto year of the run
	194	REAL,SAVE :: phisfinew(ngridmx) ! geopotential of the bedrock (= phisfi-qsurf/1000*g)
	195	REAL qsurfpal(ngridmx,nqmx) ! qsurf after a paleoclimate step : for physdem1 and restartfi
	196	REAL phisfipal(ngridmx) ! geopotential after a paleoclimate step : for physdem1 and restartfi
	197	REAL oblipal ! change of obliquity
	198	REAL peri_daypal ! new periday
	199	REAL eccpal ! change of eccentricity
	200	REAL tpalnew ! change of time
	201	REAL adjustnew ! change in N2 ice albedo
	202	REAL pdaypal ! new pday = day_ini + step
	203	REAL zdt_tot ! time range corresponding to the flux of qsurfyear
	204	REAL massacc(nqmx) ! accumulated mass
	205	REAL masslost(nqmx) ! accumulated mass
	206	REAL realarea(ngridmx) ! real area (correction poles)
	207
	208	! aerosol (ice or haze) extinction optical depth at reference wavelength
	209	! for the "naerkind" optically active aerosols:
	210	REAL aerosol(ngridmx,nlayermx,naerkind)
	211
	212	CHARACTER*80 fichier
	213	INTEGER l,ig,ierr,igout,iq,i, tapphys
	214	INTEGER lecttsoil ! lecture of tsoil from proftsoil
	215	! INTEGER iqmin ! Used if iceparty engaged
	216
	217	REAL fluxsurf_lw(ngridmx) ! incident LW (IR) surface flux (W.m-2)
	218	REAL fluxsurf_sw(ngridmx) ! incident SW (solar) surface flux (W.m-2)
	219	REAL fluxtop_lw(ngridmx) ! Outgoing LW (IR) flux to space (W.m-2)
	220	REAL fluxtop_sw(ngridmx) ! Outgoing SW (solar) flux to space (W.m-2)
	221
	222	! included by RW for equilibration test
	223	real fluxtop_dn(ngridmx)
	224	real fluxabs_sw(ngridmx) ! absorbed shortwave radiation
	225	real fluxdyn(ngridmx) ! horizontal heat transport by dynamics
	226
	227	REAL zls ! solar longitude (rad)
	228	REAL zfluxuv ! Lyman alpha flux at 1AU
	229	! ph/cm2/s
	230	REAL zday ! date (time since Ls=0, in martian days)
	231	REAL saveday ! saved date
	232	SAVE saveday
	233	REAL savedeclin ! saved declin
	234	SAVE savedeclin
	235	!REAL adjust ! adjustment for surface albedo
	236	REAL zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
	237	REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
	238
	239	c Aerosol effective radius used for radiative transfer (units=meters)
	240	REAL reffrad(ngridmx,nlayermx,naerkind)
	241
	242	c Tendencies due to various processes:
	243	REAL dqsurf(ngridmx,nqmx)
	244	REAL zdtlw(ngridmx,nlayermx) ! (K/s)
	245	REAL zdtsw(ngridmx,nlayermx) ! (K/s)
	246	REAL cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear areas
	247	REAL cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear areas
	248	REAL zdtsurf(ngridmx) ! (K/s)
	249	REAL zdtlscale(ngridmx,nlayermx)
	250	REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2)
	251	REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s)
	252	REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2)
	253	REAL zdhadj(ngridmx,nlayermx) ! (K/s)
	254	REAL zdtgw(ngridmx,nlayermx) ! (K/s)
	255	REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2)
	256	REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
	257	REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
	258	REAL tconds(ngridmx,nlayermx)
	259
	260	REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
	261	REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
	262	REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
	263	REAL zdqadj(ngridmx,nlayermx,nqmx)
	264	REAL zdqc(ngridmx,nlayermx,nqmx)
	265	REAL zdqlscale(ngridmx,nlayermx,nqmx)
	266	REAL zdqslscale(ngridmx,nqmx), zdqsc(ngridmx,nqmx)
	267	REAL zdqchim(ngridmx,nlayermx,nqmx)
	268	REAL zdqschim(ngridmx,nqmx)
	269	REAL zdqch4cloud(ngridmx,nlayermx,nqmx)
	270	REAL zdqsch4cloud(ngridmx,nqmx)
	271	REAL zdtch4cloud(ngridmx,nlayermx)
	272	REAL zdqcocloud(ngridmx,nlayermx,nqmx)
	273	REAL zdqscocloud(ngridmx,nqmx)
	274	REAL zdtcocloud(ngridmx,nlayermx)
	275	REAL rice_ch4(ngridmx,nlayermx) ! Methane ice geometric mean radius (m)
	276	REAL rice_co(ngridmx,nlayermx) ! CO ice geometric mean radius (m)
	277
	278	REAL zdqsch4fast(ngridmx) ! used only for fast model nogcm
	279	REAL zdqch4fast(ngridmx) ! used only for fast model nogcm
	280	REAL zdqscofast(ngridmx) ! used only for fast model nogcm
	281	REAL zdqcofast(ngridmx) ! used only for fast model nogcm
	282	REAL zdqflow(ngridmx,nqmx)
	283
	284	REAL zdteuv(ngridmx,nlayermx) ! (K/s)
	285	REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
	286	REAL zdumolvis(ngridmx,nlayermx)
	287	REAL zdvmolvis(ngridmx,nlayermx)
	288	real zdqmoldiff(ngridmx,nlayermx,nqmx)
	289
	290	! Haze relatated tendancies
	291	REAL zdqhaze(ngridmx,nlayermx,nqmx)
	292	REAL zdqprodhaze(ngridmx,nqmx)
	293	REAL zdqsprodhaze(ngridmx)
	294	REAL zdqhaze_col(ngridmx)
	295	REAL zdqphot_prec(ngrid,nlayer)
	296	REAL zdqphot_ch4(ngrid,nlayer)
	297	REAL zdqconv_prec(ngrid,nlayer)
	298	REAL zdq_source(ngridmx,nlayermx,nqmx)
	299	! Fast Haze relatated tendancies
	300	REAL fluxbot(ngridmx)
	301	REAL gradflux(ngridmx)
	302	REAL fluxlym_sol_bot(ngridmx) ! Solar flux Lyman alpha ph.m-2.s-1 reaching the surface
	303	REAL fluxlym_ipm_bot(ngridmx) ! IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 reaching the surface
	304	REAL flym_sol(ngridmx) ! Incident Solar flux Lyman alpha ph.m-2.s-1
	305	REAL flym_ipm(ngridmx) ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1
	306
	307	real nconsMAX
	308	real vdifcncons(ngridmx)
	309	real cadjncons(ngridmx)
	310	real ch4csMAX
	311	real ch4cncons(ngridmx)
	312	real ch4sedncons(ngridmx)
	313	real ch4sedsMAX
	314	real condncons(ngridmx)
	315	real dWtot, dWtots
	316	real dWtotn2, dWtotsn2
	317	real nconsMAXn2
	318	real vdifcnconsn2(ngridmx)
	319	real cadjnconsn2(ngridmx)
	320	real condnconsn2(ngridmx)
	321
	322	c Local variable for local intermediate calcul:
	323	REAL zflubid(ngridmx)
	324	REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
	325	REAL zdum1(ngridmx,nlayermx)
	326	REAL zdum2(ngridmx,nlayermx)
	327	REAL ztim1,ztim2,ztim3, z1,z2
	328	REAL ztime_fin
	329	REAL zdh(ngridmx,nlayermx)
	330	INTEGER length
	331	PARAMETER (length=100)
	332
	333	c local variables only used for diagnostic (output in file "diagfi" or "stats")
	334	c -----------------------------------------------------------------------------
	335	REAL ps(ngridmx), zt(ngridmx,nlayermx)
	336	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
	337	REAL zq(ngridmx,nlayermx,nqmx)
	338	REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
	339	character*2 str2
	340	character*5 str5
	341	real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
	342	real zdtdyn(ngridmx,nlayermx),ztprevious(ngridmx,nlayermx)
	343	save ztprevious
	344	real qtot1,qtot2 ! total aerosol mass
	345	integer igmin, lmin
	346	logical tdiag
	347
	348	REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
	349	REAL cd
	350	real vmr(ngridmx,nlayermx) ! volume mixing ratio
	351
	352	REAL time_phys
	353
	354	REAL tau_col(ngrid)
	355	real flux1,flux2,flux3,ts1,ts2,ts3
	356
	357	real qcol(ngridmx,nqmx)
	358	! Pluto adding variables
	359	real vmr_ch4(ngridmx) ! vmr ch4
	360	real vmr_co(ngridmx) ! vmr co
	361	real rho(ngridmx,nlayermx) ! density
	362	real zrho_ch4(ngridmx,nlayermx) ! density methane kg.m-3
	363	real zrho_co(ngridmx,nlayermx) ! density CO kg.m-3
	364	real zrho_haze(ngridmx,nlayermx) ! density haze kg.m-3
	365	real zdqrho_photprec(ngridmx,nlayermx) !photolysis rate kg.m-3.s-1
	366	real zq1temp_ch4(ngridmx) !
	367	real qsat_ch4(ngridmx) !
	368	real qsat_ch4_l1(ngridmx) !
	369
	370	! CHARACTER(LEN=20) :: txt ! to temporarly store text for eddy tracers
	371	real profmmr(ngridmx,nlayermx) ! fixed profile of haze if haze_proffix
	372	real sensiblehf1(ngridmx) ! sensible heat flux
	373	real sensiblehf2(ngridmx) ! sensible heat flux
	374
	375	! included by RW to test energy conservation
	376	REAL dEtot, dEtots, masse, vabs, dvabs, ang0
	377
	378	! included by RW to allow variations in cp with location
	379	REAL cpp3D(ngridmx,nlayermx) ! specific heat capacity at const. pressure
	380	REAL rcp3D(ngridmx,nlayermx) ! R / specific heat capacity
	381	real cppNI, rcpNI, nnu ! last one just for Seb version
	382	REAL zpopskNI(ngridmx,nlayermx)
	383
	384	! included by RW to make 1D saves not every timestep
	385	integer countG1D
	386	! integer countG3D,saveG3D
	387	save countG1D
	388	! save countG3D,saveG3D
	389	CHARACTER(len=10) :: tname
	390
	391	c=======================================================================
	392
	393	c 1. Initialisation:
	394	c -----------------
	395
	396	c 1.1 Initialisation only at first call
	397	c ---------------------------------------
	398	IF (firstcall) THEN
	399	if(ngrid.eq.1)then
	400	saveG1D=1
	401	countG1D=1
	402	endif
	403
	404	c variables set to 0
	405	c ~~~~~~~~~~~~~~~~~~
	406	call zerophys(ngrid*nlayer,dtrad)
	407	call zerophys(ngrid,fluxrad)
	408	call zerophys(ngridnlayernq,zdqsed)
	409	call zerophys(ngrid*nq, zdqssed)
	410	call zerophys(ngridnlayernq,zdqadj)
	411
	412	c initialize tracer names, indexes and properties
	413	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	414	tracerdyn=tracer ! variable tracer in dynamics
	415	IF (tracer) THEN
	416	CALL initracer()
	417	ENDIF ! end tracer
	418
	419	! Initialize aerosol indexes: done in initracer
	420	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
	421	! call iniaerosol()
	422
	423	c read startfi (initial parameters)
	424	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	425	CALL phyetat0 ("startfi.nc",0,0,
	426	& nsoilmx,nq,
	427	& day_ini,time_phys,
	428	& tsurf,tsoil,emis,q2,qsurf)
	429
	430	if (pday.ne.day_ini) then
	431	write(,) "PHYSIQ: ERROR: bad synchronization between ",
	432	& "physics and dynamics"
	433	write(,) "dynamics day: ",pday
	434	write(,) "physics day: ",day_ini
	435	stop
	436	endif
	437
	438	write (,) 'In physiq day_ini =', day_ini
	439	ptime0=ptime
	440	write (,) 'In physiq ptime0 =', ptime
	441
	442	c Initialize albedo and orbital calculation
	443	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	444	CALL surfini(ngrid,albedo)
	445	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
	446
	447	saveday=pday
	448	savedeclin=0.
	449	!adjust=0. ! albedo adjustment for convergeps
	450
	451	c Initialize surface inventory and geopotential for paleo mode
	452	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	453	if (paleo) then
	454	write(,) 'Paleo time tpal = ',tpal
	455	DO ig=1,ngrid
	456	phisfinew(ig)=phisfi(ig)-qsurf(ig,1)*g/1000. ! topo of bedrock below the ice
	457	ENDDO
	458	endif
	459
	460	if (conservn2) then
	461	write(,) 'conservn2 initial'
	462	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
	463	endif
	464
	465	c initialize soil
	466	c ~~~~~~~~~~~~~~~
	467	IF (callsoil) THEN
	468	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
	469	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	470	ELSE
	471	PRINT*,'WARNING! Thermal conduction in the soil turned off'
	472	DO ig=1,ngrid
	473	capcal(ig)=1.e5
	474	fluxgrd(ig)=0.
	475	ENDDO
	476	ENDIF
	477	icount=1
	478
	479	ENDIF ! (end of "if firstcall")
	480
	481	c ---------------------------------------------------
	482	c 1.2 Initializations done at every physical timestep:
	483	c ---------------------------------------------------
	484	c
	485	IF (ngrid.NE.ngridmx) THEN
	486	PRINT*,'STOP in PHYSIQ'
	487	PRINT*,'Probleme de dimensions :'
	488	PRINT*,'ngrid = ',ngrid
	489	PRINT*,'ngridmx = ',ngridmx
	490	STOP
	491	ENDIF
	492
	493	c Initialize various variables
	494	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	495	call zerophys(ngrid*nlayer, pdv)
	496	call zerophys(ngrid*nlayer, pdu)
	497	call zerophys(ngrid*nlayer,pdt)
	498	call zerophys(ngridnlayernq, pdq)
	499	call zerophys(ngrid, pdpsrf)
	500	call zerophys(ngrid, zflubid)
	501	call zerophys(ngrid, zdtsurf)
	502	call zerophys(ngrid*nq, dqsurf)
	503
	504	if (conservn2) then
	505	write(,) 'conservn2 iniloop'
	506	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
	507	endif
	508
	509	igout=ngrid/2+1
	510
	511	zday=pday+ptime ! compute time, in sols (and fraction thereof)
	512
	513	c Compute Solar Longitude (Ls) :
	514	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	515	if (season) then
	516	call solarlong(zday,zls)
	517	else
	518	call solarlong(float(day_ini),zls)
	519	end if
	520
	521	c Get Lyman alpha flux at specific Ls
	522	if (haze) then
	523	call lymalpha(zls,zfluxuv)
	524	print*, 'Haze lyman-alpha zls,zfluxuv=',zls,zfluxuv
	525	end if
	526
	527	c Compute geopotential between layers
	528	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	529	DO l=1,nlayer
	530	DO ig=1,ngrid
	531	zzlay(ig,l)=pphi(ig,l)/g
	532	ENDDO
	533	ENDDO
	534	DO ig=1,ngrid
	535	zzlev(ig,1)=0.
	536	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
	537	ENDDO
	538	DO l=2,nlayer
	539	DO ig=1,ngrid
	540	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
	541	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
	542	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	543	ENDDO
	544	ENDDO
	545
	546	! Potential temperature calculation not the same in physiq and dynamic...
	547	c Compute potential temperature zh
	548	DO l=1,nlayer
	549	DO ig=1,ngrid
	550	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
	551	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
	552	ENDDO
	553	ENDDO
	554
	555	c --------------------------------------------------------
	556	c 1.3 thermosphere
	557	c --------------------------------------------------------
	558
	559	c ajout de la conduction depuis la thermosphere
	560	IF (callconduct) THEN
	561
	562	call conduction (ngrid,nlayer,ptimestep,
	563	$ pplay,pt,zzlay,zzlev,zdtconduc,tsurf)
	564	DO l=1,nlayer
	565	DO ig=1,ngrid
	566	pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
	567	ENDDO
	568	ENDDO
	569
	570	ENDIF
	571
	572	! ajout de la viscosite moleculaire
	573	IF (callmolvis) THEN
	574	call molvis(ngrid,nlayer,ptimestep,
	575	$ pplay,pt,zzlay,zzlev,
	576	$ zdtconduc,pu,tsurf,zdumolvis)
	577	call molvis(ngrid,nlayer,ptimestep,
	578	$ pplay,pt,zzlay,zzlev,
	579	$ zdtconduc,pv,tsurf,zdvmolvis)
	580
	581	DO l=1,nlayer
	582	DO ig=1,ngrid
	583	! pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
	584	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
	585	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
	586	ENDDO
	587	ENDDO
	588	ENDIF
	589
	590	IF (callmoldiff) THEN
	591	call moldiff_red(ngrid,nlayer,nq,
	592	& pplay,pplev,pt,pdt,pq,pdq,ptimestep,
	593	& zzlay,zdtconduc,zdqmoldiff)
	594
	595	DO l=1,nlayer
	596	DO ig=1,ngrid
	597	DO iq=1, nq
	598	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
	599	ENDDO
	600	ENDDO
	601	ENDDO
	602	ENDIF
	603
	604	if (conservn2) then
	605	write(,) 'conservn2 thermo'
	606	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
	607	endif
	608
	609	c-----------------------------------------------------------------------
	610	c 2. Compute radiative tendencies :
	611	c--------- --------------------------------------------------------------
	612	c Saving qsurf to compute paleo flux condensation/sublimation
	613	DO iq=1, nq
	614	DO ig=1,ngrid
	615	IF (qsurf(ig,iq).lt.0.) then
	616	qsurf(ig,iq)=0.
	617	ENDIF
	618	qsurf1(ig,iq)=qsurf(ig,iq)
	619	ENDDO
	620	ENDDO
	621
	622	IF (callrad) THEN
	623	IF( MOD(icount-1,iradia).EQ.0) THEN
	624
	625	c Local stellar zenith angle
	626	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	627	IF (triton) then
	628	CALL orbitetriton(zls,zday,dist_sol,declin)
	629	ELSE
	630	CALL orbite(zls,dist_sol,declin)
	631	ENDIF
	632
	633	IF (diurnal) THEN
	634	ztim1=SIN(declin)
	635	ztim2=COS(declin)COS(2.pi*(zday-.5))
	636	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
	637
	638	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
	639	s ztim1,ztim2,ztim3,mu0,fract)
	640
	641	ELSE
	642	CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
	643	! WARNING: this function appears not to work in 1D
	644	ENDIF
	645
	646	c Albedo /IT changes depending of surface ices (only in 3D)
	647	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	648	if (ngrid.ne.1) then
	649
	650	!! Specific Triton to change albedo of N2 so that Psurf
	651	!! converges toward 1.4 Pa in "1989" seasons
	652	if (convergeps.and.triton) then
	653	! 1989 declination
	654	if (declin180./pi.gt.-46..and.declin180./pi.lt.-45.
	655	& .and.zday.gt.saveday+1000..and.declin.lt.savedeclin) then
	656	call globalaverage2d(ngrid,pplev(:,1),globave)
	657	if (globave.gt.1.5) then
	658	adjust=adjust+0.005
	659	else if (globave.lt.1.3) then
	660	adjust=adjust-0.005
	661	endif
	662	saveday=zday
	663	endif
	664	call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,
	665	& capcal,adjust,dist_sol,albedo,emis,flusurfold,ptimestep,
	666	& zls)
	667	savedeclin=declin
	668
	669	else
	670	! Surface properties
	671	call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,
	672	& capcal,adjust,dist_sol,albedo,emis,flusurfold,ptimestep,
	673	& zls)
	674	end if
	675	!TB22
	676	else ! here ngrid=1
	677	call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,
	678	& capcal,adjust,dist_sol,albedo,emis,flusurfold,ptimestep,
	679	& zls)
	680
	681	end if ! if ngrid ne 1
	682
	683	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	684	c Fixed zenith angle in 1D
	685	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	686	if(ngrid.eq.1.and.globmean1d) then
	687	!global mean 1D radiative equiilium
	688	ig=1
	689	mu0(ig) = 2**(-0.5)
	690	fract(ig)= 1/(4*mu0(ig))
	691	!print*,'WARNING GLOBAL MEAN CALCULATION'
	692	endif
	693
	694	c Call main radiative transfer scheme
	695	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	696
	697	c Radiative transfer
	698	c ------------------
	699	if (corrk) then
	700	c print*,'TB22 start corrk',tsurf ! Bug testphys1d nodebug
	701	CALL callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,
	702	& albedo,emis,mu0,pplev,pplay,pt,
	703	& tsurf,fract,dist_sol,igout,aerosol,cpp3D,
	704	& zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,
	705	& fluxtop_sw,fluxtop_dn,reffrad,tau_col,ptime,pday,
	706	& firstcall,lastcall,zzlay)
	707
	708	c Radiative flux from the sky absorbed by the surface (W.m-2)
	709	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	710	DO ig=1,ngrid
	711	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
	712	& +fluxsurf_sw(ig)*(1.-albedo(ig))
	713	ENDDO
	714
	715	c Net atmospheric radiative heating/cooling rate (K.s-1):
	716	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	717	DO l=1,nlayer
	718	DO ig=1,ngrid
	719	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
	720	ENDDO
	721	ENDDO
	722
	723	else ! corrk = F
	724
	725	c Atmosphere has no radiative effect
	726	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	727	DO ig=1,ngrid
	728	fluxtop_dn(ig) = fract(ig)mu0(ig)Fat1AU/dist_sol**2
	729	fluxrad_sky(ig) = fluxtop_dn(ig)*(1.-albedo(ig))
	730	fluxtop_sw(ig) = fluxtop_dn(ig)*albedo(ig)
	731	fluxtop_lw(ig) = emis(ig)stephantsurf(ig)**4
	732	ENDDO ! radiation skips the atmosphere entirely
	733	DO l=1,nlayer
	734	DO ig=1,ngrid
	735	dtrad(ig,l)=0.0
	736	ENDDO
	737	ENDDO ! hence no atmospheric radiative heating
	738
	739	endif ! if corrk
	740
	741	ENDIF ! of if(mod(icount-1,iradia).eq.0) ie radiative timestep
	742
	743	! Transformation of the radiative tendencies:
	744	! -------------------------------------------
	745
	746	! Net radiative surface flux (W.m-2)
	747	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	748
	749	DO ig=1,ngrid
	750	zplanck(ig)=tsurf(ig)*tsurf(ig)
	751	zplanck(ig)=emis(ig)*
	752	& stephanzplanck(ig)zplanck(ig)
	753	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
	754	ENDDO
	755
	756	DO l=1,nlayer
	757	DO ig=1,ngrid
	758	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
	759	ENDDO
	760	ENDDO
	761
	762	ENDIF ! of IF (callrad)
	763
	764	if (conservn2) then
	765	write(,) 'conservn2 radiat'
	766	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
	767	endif
	768
	769	!-----------------------------------------------------------------------
	770	! 4. Vertical diffusion (turbulent mixing):
	771	!-----------------------------------------------------------------------
	772
	773	IF (calldifv) THEN
	774
	775	DO ig=1,ngrid
	776	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
	777	ENDDO
	778
	779	CALL zerophys(ngrid*nlayer,zdum1)
	780	CALL zerophys(ngrid*nlayer,zdum2)
	781	do l=1,nlayer
	782	do ig=1,ngrid
	783	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	784	enddo
	785	enddo
	786
	787	c Calling vdif (Martian version WITH N2 condensation)
	788	CALL vdifc(ngrid,nlayer,nq,zpopsk,
	789	& ptimestep,capcal,lwrite,
	790	& pplay,pplev,zzlay,zzlev,z0,
	791	& pu,pv,zh,pq,pt,tsurf,emis,qsurf,
	792	& zdum1,zdum2,zdh,pdq,pdt,zflubid,
	793	& zdudif,zdvdif,zdhdif,zdtsdif,q2,
	794	& zdqdif,zdqsdif,qsat_ch4,qsat_ch4_l1) !,zq1temp_ch4,qsat_ch4)
	795
	796	DO l=1,nlayer
	797	DO ig=1,ngrid
	798	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
	799	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
	800	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
	801
	802	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
	803	ENDDO
	804	ENDDO
	805
	806	DO ig=1,ngrid
	807	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
	808	ENDDO
	809
	810	bcond=1./tcond1p4Pa
	811	acond=r/lw_n2
	812
	813	if (tracer) then
	814	DO iq=1, nq
	815	DO l=1,nlayer
	816	DO ig=1,ngrid
	817	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
	818	ENDDO
	819	ENDDO
	820	ENDDO
	821
	822	DO iq=1, nq
	823	DO ig=1,ngrid
	824	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
	825	ENDDO
	826	ENDDO
	827
	828	end if ! of if (tracer)
	829
	830	!------------------------------------------------------------------
	831	! test methane conservation
	832	c if(methane)then
	833	c call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
	834	c & ptimestep,pplev,zdqdif,zdqsdif,'CH4',' vdifc ')
	835	c endif ! methane
	836	!------------------------------------------------------------------
	837	! test CO conservation
	838	c if(carbox)then
	839	c call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
	840	c & ptimestep,pplev,zdqdif,zdqsdif,'CO ',' vdifc ')
	841	c endif ! carbox
	842	!------------------------------------------------------------------
	843
	844	ELSE ! case with calldifv=F
	845	ztim1=4.stephanptimestep
	846	DO ig=1,ngrid
	847	ztim2=ztim1emis(ig)tsurf(ig)**3
	848	z1=capcal(ig)*tsurf(ig)+
	849	s ztim2tsurf(ig)+ (fluxrad(ig)+fluxgrd(ig))ptimestep
	850	z2= capcal(ig)+ztim2
	851	zdtsurf(ig)=(z1/z2 - tsurf(ig))/ptimestep
	852	if (ig.eq.475) print*, 'TB22 z1,z2=',z1,z2,z1/z2
	853	if (ig.eq.475) print*, 'TB22 cc,ee=',capcal(ig),emis(ig)
	854	if (ig.eq.475) print*, 'TB22 fr,fg=',fluxrad(ig),fluxgrd(ig)
	855
	856	! for output:
	857	!dplanck(ig)=4.stephanptimestepemis(ig)tsurf(ig)**3
	858	ENDDO
	859
	860	c ------------------------------------------------------------------
	861	c Methane surface sublimation and condensation in fast model (nogcm)
	862	c ------------------------------------------------------------------
	863	print*, 'TB22 tsurf ini=',tsurf(475)
	864	print, 'TB22 dt ini=',zdtsurf(475)ptimestep
	865	print*, 'TB22 qsurfn2=',qsurf(475,1)
	866	print*, 'TB22 qsurfch4=',qsurf(475,igcm_ch4_ice)
	867	if ((methane).and.(fast).and.condmetsurf) THEN
	868
	869	call ch4surf(ngrid,nlayer,nq,ptimestep,
	870	& tsurf,zdtsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf,
	871	& zdqch4fast,zdqsch4fast)
	872
	873	DO ig=1,ngrid
	874	dqsurf(ig,igcm_ch4_ice)= dqsurf(ig,igcm_ch4_ice) +
	875	& zdqsch4fast(ig)
	876	pdq(ig,1,igcm_ch4_gas)= pdq(ig,1,igcm_ch4_gas) +
	877	& zdqch4fast(ig)
	878	zdtsurf(ig)=zdtsurf(ig)+lw_ch4*zdqsch4fast(ig)/capcal(ig)
	879	ENDDO
	880	end if
	881	print, lati(475)180./pi,long(475)*180./pi,capcal(475)
	882	print, 'dtch4=',lw_ch4zdqsch4fast(475)/capcal(475)*ptimestep
	883	print, 'dqch4s=',zdqsch4fast(475)ptimestep
	884	print*, 'newqsurfch4=',qsurf(475,igcm_ch4_ice)+
	885	& dqsurf(475,igcm_ch4_ice)*ptimestep
	886	print, 'TB22 tsurf new=',tsurf(475)+zdtsurf(475)ptimestep
	887	c ------------------------------------------------------------------
	888	c CO surface sublimation and condensation in fast model (nogcm)
	889	c ------------------------------------------------------------------
	890	if ((carbox).and.(fast).and.condcosurf) THEN
	891
	892	call cosurf(ngrid,nlayer,nq,ptimestep,
	893	& tsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf,
	894	& zdqcofast,zdqscofast)
	895
	896	DO ig=1,ngrid
	897	dqsurf(ig,igcm_co_ice)= dqsurf(ig,igcm_co_ice) +
	898	& zdqscofast(ig)
	899	pdq(ig,1,igcm_co_gas)= pdq(ig,1,igcm_co_gas) +
	900	& zdqcofast(ig)
	901	zdtsurf(ig)=zdtsurf(ig)+lw_co*zdqscofast(ig)/capcal(ig)
	902	ENDDO
	903	end if
	904
	905	ENDIF ! of IF (calldifv)
	906
	907	if (conservn2) then
	908	write(,) 'conservn2 diff'
	909	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)+
	910	& dqsurf(:,1)*ptimestep)
	911	endif
	912
	913	!------------------------------------------------------------------
	914	! test CO conservation
	915	! if(carbox)then
	916	! call testconservfast(ngrid,nlayer,nq,
	917	! & ptimestep,pplev,zdqcofast,zdqscofast,'CO ',' vdifc ')
	918	! endif ! carbox
	919	!------------------------------------------------------------------
	920
	921	!-----------------------------------------------------------------------
	922	! 5. Dry convective adjustment:
	923	! -----------------------------
	924
	925	IF(calladj) THEN
	926
	927	DO l=1,nlayer
	928	DO ig=1,ngrid
	929	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	930	ENDDO
	931	ENDDO
	932	CALL zerophys(ngrid*nlayer,zduadj)
	933	CALL zerophys(ngrid*nlayer,zdvadj)
	934	CALL zerophys(ngrid*nlayer,zdhadj)
	935
	936	CALL convadj(ngrid,nlayer,nq,ptimestep,
	937	& pplay,pplev,zpopsk,
	938	& pu,pv,zh,pq,
	939	& pdu,pdv,zdh,pdq,
	940	& zduadj,zdvadj,zdhadj,
	941	& zdqadj)
	942
	943
	944	DO l=1,nlayer
	945	DO ig=1,ngrid
	946	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
	947	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
	948	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
	949	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
	950	ENDDO
	951	ENDDO
	952
	953	if(tracer) then
	954	DO iq=1, nq
	955	DO l=1,nlayer
	956	DO ig=1,ngrid
	957	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
	958	ENDDO
	959	ENDDO
	960	ENDDO
	961	end if
	962
	963	ENDIF ! of IF(calladj)
	964
	965	!------------------------------------------------------------------
	966	! test methane conservation
	967	c if(methane)then
	968	c call testchange(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
	969	c & ptimestep,pplev,zdqadj,'CH4','calladj')
	970	c endif ! methane
	971	!------------------------------------------------------------------
	972	! test CO conservation
	973	c if(carbox)then
	974	c call testchange(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
	975	c & ptimestep,pplev,zdqadj,'CO ','calladj')
	976	c endif ! carbox
	977	!------------------------------------------------------------------
	978
	979	!-----------------------------------------------------------------------
	980	! 6. Nitrogen condensation-sublimation:
	981	! -------------------------------------------
	982
	983	IF (N2cond) THEN
	984	if (.not.tracer.or.igcm_n2.eq.0) then
	985	print*,'We need a n2 ice tracer to condense n2'
	986	call abort
	987	endif
	988
	989	c write(,) 'before N2 condens :'
	990	c call massn2(ngrid,nlayer,pplev,qsurf(:,1),dqsurf(:,igcm_n2),
	991	c & pdpsrf,ptimestep)
	992
	993	if (tracer) then
	994	zdqc(:,:,:)=0.
	995	zdqsc(:,:)=0.
	996	call condens_n2(ngrid,nlayer,nq,ptimestep,
	997	& capcal,pplay,pplev,tsurf,pt,
	998	& pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
	999	& qsurf(1,igcm_n2),albedo,emis,
	1000	& zdtc,zdtsurfc,pdpsrf,zduc,zdvc,
	1001	& zdqc,zdqsc(1,igcm_n2))
	1002	endif
	1003
	1004	!!! TB temporaire
	1005	!zdtc=0.
	1006	!zdvc=0.
	1007	!zduc=0.
	1008	!zdqc=0.
	1009	!zdqsc=0.
	1010	!zdtsurfc=0.
	1011
	1012	DO l=1,nlayer
	1013	DO ig=1,ngrid
	1014	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
	1015	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
	1016	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
	1017	ENDDO
	1018	ENDDO
	1019	DO ig=1,ngrid
	1020	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
	1021	ENDDO
	1022
	1023	DO iq=1, nq
	1024	c TB: option eddy not condensed with N2
	1025	c txt=noms(iq)
	1026	c IF (txt(1:4).ne."eddy") THEN
	1027	DO l=1,nlayer
	1028	DO ig=1,ngrid
	1029	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
	1030	ENDDO
	1031	ENDDO
	1032	c ENDIF
	1033	ENDDO
	1034	DO ig=1,ngrid
	1035	dqsurf(ig,igcm_n2)= dqsurf(ig,igcm_n2) + zdqsc(ig,igcm_n2)
	1036	ENDDO
	1037
	1038	ENDIF ! of IF (N2cond)
	1039
	1040	if (conservn2) then
	1041	write(,) 'conservn2 n2cond'
	1042	call testconservmass(ngrid,nlayer,pplev(:,1)+
	1043	& pdpsrf(:)ptimestep,qsurf(:,1)+dqsurf(:,1)ptimestep)
	1044	endif
	1045
	1046	!------------------------------------------------------------------
	1047	! test n2 conservation for one tendency / pplevis not updated here w pdpsrf
	1048	! if(tracer)then
	1049	! call testconserv(ngrid,nlayer,nq,igcm_n2,igcm_n2,
	1050	! & ptimestep,pplev,zdqc,zdqsc,'N2 ',' n2cond')
	1051	! call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
	1052	! endif ! n2
	1053	!------------------------------------------------------------------
	1054	! test methane conservation
	1055	! if(methane)then
	1056	! call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
	1057	! & ptimestep,pplev,zdqc,zdqsc,'CH4',' n2cond')
	1058	! endif ! methane
	1059	!------------------------------------------------------------------
	1060	! test CO conservation
	1061	c if(carbox)then
	1062	c call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
	1063	c & ptimestep,pplev,zdqc,zdqsc,'CO ',' n2cond')
	1064	c endif ! carbox
	1065	!------------------------------------------------------------------
	1066
	1067	c-----------------------------------------------------------------------
	1068	c 7. Specific parameterizations for tracers
	1069	c -----------------------------------------
	1070
	1071	if (tracer) then
	1072
	1073	c 7a. Methane, CO, and ice
	1074	c ---------------------------------------
	1075	c Methane ice condensation in the atmosphere
	1076	c ----------------------------------------
	1077	zdqch4cloud(:,:,:)=0.
	1078	if ((methane).and.(metcloud).and.(.not.fast)) THEN
	1079	call ch4cloud(ngrid,nlayer,naerkind,ptimestep,
	1080	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
	1081	& pq,pdq,zdqch4cloud,zdqsch4cloud,zdtch4cloud,
	1082	& nq,rice_ch4)
	1083
	1084	DO l=1,nlayer
	1085	DO ig=1,ngrid
	1086	pdq(ig,l,igcm_ch4_gas)=pdq(ig,l,igcm_ch4_gas)+
	1087	& zdqch4cloud(ig,l,igcm_ch4_gas)
	1088	pdq(ig,l,igcm_ch4_ice)=pdq(ig,l,igcm_ch4_ice)+
	1089	& zdqch4cloud(ig,l,igcm_ch4_ice)
	1090	ENDDO
	1091	ENDDO
	1092
	1093	! Increment methane ice surface tracer tendency
	1094	DO ig=1,ngrid
	1095	dqsurf(ig,igcm_ch4_ice)=dqsurf(ig,igcm_ch4_ice)+
	1096	& zdqsch4cloud(ig,igcm_ch4_ice)
	1097	ENDDO
	1098
	1099	! update temperature tendancy
	1100	DO ig=1,ngrid
	1101	DO l=1,nlayer
	1102	pdt(ig,l)=pdt(ig,l)+zdtch4cloud(ig,l)
	1103	ENDDO
	1104	ENDDO
	1105	end if
	1106
	1107	c ---------------------------------------
	1108	c CO ice condensation in the atmosphere
	1109	c ----------------------------------------
	1110	zdqcocloud(:,:,:)=0.
	1111	IF ((carbox).and.(monoxcloud).and.(.not.fast)) THEN
	1112	call cocloud(ngrid,nlayer,naerkind,ptimestep,
	1113	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
	1114	& pq,pdq,zdqcocloud,zdqscocloud,zdtcocloud,
	1115	& nq,rice_co,qsurf(1,igcm_n2),dqsurf(1,igcm_n2))
	1116
	1117	DO l=1,nlayer
	1118	DO ig=1,ngrid
	1119	pdq(ig,l,igcm_co_gas)=pdq(ig,l,igcm_co_gas)+
	1120	& zdqcocloud(ig,l,igcm_co_gas)
	1121	pdq(ig,l,igcm_co_ice)=pdq(ig,l,igcm_co_ice)+
	1122	& zdqcocloud(ig,l,igcm_co_ice)
	1123	ENDDO
	1124	ENDDO
	1125
	1126	! Increment CO ice surface tracer tendency
	1127	DO ig=1,ngrid
	1128	dqsurf(ig,igcm_co_ice)=dqsurf(ig,igcm_co_ice)+
	1129	& zdqscocloud(ig,igcm_co_ice)
	1130	ENDDO
	1131
	1132	! update temperature tendancy
	1133	DO ig=1,ngrid
	1134	DO l=1,nlayer
	1135	pdt(ig,l)=pdt(ig,l)+zdtcocloud(ig,l)
	1136	ENDDO
	1137	ENDDO
	1138
	1139	END IF ! of IF (carbox)
	1140
	1141	!------------------------------------------------------------------
	1142	! test methane conservation
	1143	c if(methane)then
	1144	c call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
	1145	c & ptimestep,pplev,zdqch4cloud,zdqsch4cloud,'CH4','ch4clou')
	1146	c endif ! methane
	1147	!------------------------------------------------------------------
	1148	! test CO conservation
	1149	c if(carbox)then
	1150	c call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
	1151	c & ptimestep,pplev,zdqcocloud,zdqscocloud,'CO ','cocloud')
	1152	c endif ! carbox
	1153	!------------------------------------------------------------------
	1154
	1155	c 7b. Haze particle production
	1156	c -------------------
	1157	IF (haze) THEN
	1158
	1159	zdqphot_prec(:,:)=0.
	1160	zdqphot_ch4(:,:)=0.
	1161	! Forcing to a fixed haze profile if haze_proffix
	1162	if (haze_proffix.and.i_haze.gt.0.) then
	1163	call haze_prof(ngrid,nlayer,zzlay,pplay,pt,
	1164	& reffrad,profmmr)
	1165	zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze))
	1166	& /ptimestep
	1167	else
	1168	call hazecloud(ngrid,nlayer,nq,ptimestep,
	1169	& pplay,pplev,pq,pdq,dist_sol,mu0,zfluxuv,zdqhaze,
	1170	& zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
	1171	endif
	1172
	1173	DO iq=1, nq ! should be updated
	1174	DO l=1,nlayer
	1175	DO ig=1,ngrid
	1176	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqhaze(ig,l,iq)
	1177	ENDDO
	1178	ENDDO
	1179	ENDDO
	1180
	1181	ENDIF
	1182
	1183	IF (fast.and.fasthaze) THEN
	1184	call prodhaze(ngrid,nlayer,nq,ptimestep,pplev,pq,pdq,dist_sol,
	1185	& mu0,declin,zdqprodhaze,zdqsprodhaze,gradflux,fluxbot,
	1186	& fluxlym_sol_bot,fluxlym_ipm_bot,flym_sol,flym_ipm)
	1187
	1188	DO ig=1,ngrid
	1189	pdq(ig,1,igcm_ch4_gas)=pdq(ig,1,igcm_ch4_gas)+
	1190	& zdqprodhaze(ig,igcm_ch4_gas)
	1191	pdq(ig,1,igcm_prec_haze)=pdq(ig,1,igcm_prec_haze)+
	1192	& zdqprodhaze(ig,igcm_prec_haze)
	1193	pdq(ig,1,igcm_haze)=abs(pdq(ig,1,igcm_haze)+
	1194	& zdqprodhaze(ig,igcm_haze))
	1195	qsurf(ig,igcm_haze)= qsurf(ig,igcm_haze)+
	1196	& zdqsprodhaze(ig)*ptimestep
	1197	ENDDO
	1198
	1199	ENDIF
	1200
	1201	c 7c. Aerosol particles
	1202	c -------------------
	1203
	1204	c -------------
	1205	c Sedimentation
	1206	c -------------
	1207	IF (sedimentation) THEN
	1208	call zerophys(ngridnlayernq, zdqsed)
	1209	call zerophys(ngrid*nq, zdqssed)
	1210
	1211	call callsedim(ngrid,nlayer, ptimestep,
	1212	& pplev,zzlev, pt,rice_ch4,rice_co,
	1213	& pq, pdq, zdqsed, zdqssed,nq,pphi)
	1214
	1215	DO iq=1, nq
	1216	DO l=1,nlayer
	1217	DO ig=1,ngrid
	1218	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
	1219	ENDDO
	1220	ENDDO
	1221	ENDDO
	1222	DO iq=2, nq
	1223	DO ig=1,ngrid
	1224	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
	1225	ENDDO
	1226	ENDDO
	1227
	1228	END IF ! of IF (sedimentation)
	1229
	1230	!------------------------------------------------------------------
	1231	! test methane conservation
	1232	c if(methane)then
	1233	c call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
	1234	c & ptimestep,pplev,zdqsed,zdqssed,'CH4',' sedim ')
	1235	c endif ! methane
	1236	!------------------------------------------------------------------
	1237	! test CO conservation
	1238	c if(carbox)then
	1239	c call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
	1240	c & ptimestep,pplev,zdqsed,zdqssed,'CO ',' sedim ')
	1241	c endif ! carbox
	1242	!------------------------------------------------------------------
	1243
	1244	c 7d. Updates
	1245	c ---------
	1246	! write(,) 'before update qsurf:'
	1247	! call massn2(ngrid,nlayer,pplev,qsurf(:,1),dqsurf(:,igcm_n2),
	1248	! & pdpsrf,ptimestep)
	1249
	1250	c ---------------------------------
	1251	c Updating tracer budget on surface (before spread of glacier)
	1252	c ---------------------------------
	1253	DO iq=1, nq
	1254	DO ig=1,ngrid
	1255	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
	1256	ENDDO
	1257	ENDDO
	1258
	1259	endif ! of if (tracer)
	1260
	1261	if (conservn2) then
	1262	write(,) 'conservn2 tracer'
	1263	call testconservmass(ngrid,nlayer,pplev(:,1)+
	1264	& pdpsrf(:)*ptimestep,qsurf(:,1))
	1265	endif
	1266
	1267	DO ig=1,ngrid
	1268	flusurf(ig,igcm_n2)=(qsurf(ig,igcm_n2)-
	1269	& qsurf1(ig,igcm_n2))/ptimestep
	1270	flusurfold(ig,igcm_n2)=flusurf(ig,igcm_n2)
	1271	if (methane) then
	1272	flusurf(ig,igcm_ch4_ice)=(qsurf(ig,igcm_ch4_ice)-
	1273	& qsurf1(ig,igcm_ch4_ice))/ptimestep
	1274	flusurfold(ig,igcm_ch4_ice)=flusurf(ig,igcm_ch4_ice)
	1275	endif
	1276	if (carbox) then
	1277	flusurf(ig,igcm_co_ice)=(qsurf(ig,igcm_co_ice)-
	1278	& qsurf1(ig,igcm_co_ice))/ptimestep
	1279	!flusurfold(ig,igcm_co_ice)=flusurf(ig,igcm_co_ice)
	1280	endif
	1281	ENDDO
	1282
	1283	!! Special source of haze particle !
	1284	! todo: should be placed in haze and use tendency of n2 instead of flusurf
	1285	IF (source_haze) THEN
	1286	call hazesource(ngrid,nlayer,nq,ptimestep,
	1287	& pplev,flusurf,mu0,zdq_source)
	1288
	1289	DO iq=1, nq
	1290	DO l=1,nlayer
	1291	DO ig=1,ngrid
	1292	pdq(ig,l,iq)=pdq(ig,l,iq)+zdq_source(ig,l,iq)
	1293	ENDDO
	1294	ENDDO
	1295	ENDDO
	1296	ENDIF
	1297
	1298	!-----------------------------------------------------------------------
	1299	! 9. Surface and sub-surface soil temperature
	1300	!-----------------------------------------------------------------------
	1301	! For diagnostic
	1302	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1303	DO ig=1,ngrid
	1304	sensiblehf1(ig)=rho(ig,1)cpp(0.4/log(zzlay(ig,1)/z0))*2
	1305	& (pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1))*0.5
	1306	& (tsurf(ig)-pt(ig,1))
	1307	sensiblehf2(ig)=zflubid(ig)-capcal(ig)*zdtsdif(ig)
	1308
	1309	ENDDO
	1310
	1311	! 9.1 Increment Surface temperature:
	1312	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1313	DO ig=1,ngrid
	1314	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
	1315	if (tsurf(ig).lt.0.) then
	1316	print*,'TB22 tsurf=',tsurf(ig),ig
	1317	stop
	1318	endif
	1319	ENDDO
	1320	!
	1321	! 9.2 Compute soil temperatures and subsurface heat flux:
	1322	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1323	IF (callsoil) THEN
	1324	CALL soil(ngrid,nsoilmx,.false.,tidat,
	1325	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
	1326	ENDIF
	1327
	1328	! For output :
	1329	tidat_out(:,:)=0.
	1330	DO l=1,min(nlayermx,nsoilmx)
	1331	tidat_out(:,l)=tidat(:,l)
	1332	ENDDO
	1333
	1334	! 9.3 multiply tendencies of cond/subli for paleo loop only in the
	1335	! last Pluto year of the simulation
	1336	! Year day must be adapted in the startfi for each object
	1337	! Paleo uses year_day to calculate the annual mean tendancies
	1338	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1339	IF (paleo) then
	1340	if (zday.gt.day_ini+ptime0+nday-year_day) then
	1341	DO iq=1,nq
	1342	DO ig=1,ngrid
	1343	qsurfyear(ig,iq)=qsurfyear(ig,iq)+
	1344	& (qsurf(ig,iq)-qsurf1(ig,iq)) !kg m-2 !ptimestep
	1345	ENDDO
	1346	ENDDO
	1347	endif
	1348	endif
	1349
	1350	! 9.4 Glacial flow at each timestep glastep or at lastcall
	1351	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1352	IF (fast.and.glaflow) THEN
	1353	if ((mod(zday-day_ini-ptime0,glastep)).lt.1.
	1354	& .or.lastcall) then
	1355	IF (lastcall) then
	1356	dstep=mod(zday-day_ini-ptime0,glastep)*daysec
	1357	else
	1358	dstep=glastep*daysec
	1359	endif
	1360	zdqflow(:,:)=qsurf(:,:)
	1361	IF (paleo) then
	1362	!print*, 'calling glacier',dstep
	1363	call spreadglacier_paleo(ngrid,nq,qsurf,
	1364	& phisfinew,dstep,tsurf)
	1365	else
	1366	call spreadglacier_simple(ngrid,nq,qsurf,dstep)
	1367	endif
	1368	zdqflow(:,:)=(zdqflow(:,:)-qsurf(:,:))/dstep
	1369
	1370	if (conservn2) then
	1371	write(,) 'conservn2 glaflow'
	1372	call testconservmass(ngrid,nlayer,pplev(:,1)+
	1373	& pdpsrf(:)*ptimestep,qsurf(:,1))
	1374	endif
	1375
	1376	endif
	1377	ENDIF
	1378
	1379	!-----------------------------------------------------------------------
	1380	! 10. Perform diagnostics and write output files
	1381	!-----------------------------------------------------------------------
	1382	! ---------------------------------------------------------
	1383	! Check the energy balance of the simulation during the run
	1384	! ---------------------------------------------------------
	1385	flux1 = 0
	1386	flux2 = 0
	1387	flux3 = 0
	1388	do ig=1,ngrid
	1389
	1390	flux1 = flux1 +
	1391	& area(ig)*(fluxtop_dn(ig) - fluxtop_sw(ig))
	1392	flux2 = flux2 + area(ig)*fluxtop_lw(ig)
	1393	flux3 = flux3 + area(ig)*fluxtop_dn(ig)
	1394	fluxabs_sw(ig)=fluxtop_dn(ig) - fluxtop_sw(ig)
	1395
	1396	end do
	1397	! print*,'Incident solar flux, absorbed solar flux, OLR (W/m2)'
	1398	! print*, flux3/totarea, ' ', flux1/totarea ,
	1399	! & ' = ', flux2/totarea
	1400
	1401	if(meanOLR)then
	1402	! to record global radiative balance
	1403	open(92,file="rad_bal.out",form='formatted',access='append')
	1404	write(92,*) zday,flux3/totarea,flux1/totarea,flux2/totarea
	1405	close(92)
	1406	endif
	1407
	1408	! Surface temperature information
	1409	ts1 = 0
	1410	ts2 = 999
	1411	ts3 = 0
	1412	do ig=1,ngrid
	1413	ts1 = ts1 + area(ig)*tsurf(ig)
	1414	ts2 = min(ts2,tsurf(ig))
	1415	ts3 = max(ts3,tsurf(ig))
	1416	end do
	1417	! print*,'Mean Tsurf =',ts1/totarea , ' Min Tsurf=',ts2,
	1418	! & ' Max Tsurf =',ts3
	1419	! print*,'Max Tsurf=',ts3,'Min Tsurf=',ts2
	1420
	1421	! -------------------------------
	1422	! Dynamical fields incrementation
	1423	! -------------------------------
	1424	! (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
	1425	! temperature, zonal and meridional wind
	1426	DO l=1,nlayer
	1427	DO ig=1,ngrid
	1428	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
	1429	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
	1430	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
	1431	! diagnostic
	1432	zdtdyn(ig,l)=ztprevious(ig,l)-pt(ig,l)
	1433	ztprevious(ig,l)=zt(ig,l)
	1434	ENDDO
	1435	ENDDO
	1436	! if(firstcall)call zerophys(zdtdyn)
	1437	if(firstcall) zdtdyn(:,:)=0.
	1438	! dynamical heating diagnostic
	1439	DO ig=1,ngrid
	1440	fluxdyn(ig)=0.
	1441	DO l=1,nlayer
	1442	fluxdyn(ig)=fluxdyn(ig) - (zdtdyn(ig,l)/ptimestep)
	1443	& (pplev(ig,l)-pplev(ig,l+1))cpp/g
	1444	ENDDO
	1445	ENDDO
	1446
	1447	! -------------------------------
	1448	! Other diagnostics
	1449	! -------------------------------
	1450	!! diagnostic qsat_lev1 (ps, T1)
	1451	! CALL methanesat(ngridmx,zt(:,1),pplev(1,1),qsat_lev1(:),
	1452	! & qsurf(:,igcm_n2))
	1453
	1454	!! eddy tracers : decay time
	1455	! DO l=1,nlayer
	1456	! DO ig=1,ngrid
	1457	! pdq(ig,l,igcm_eddy1e6)=pdq(ig,l,igcm_eddy1e6)-
	1458	! & pq(ig,l,igcm_eddy1e6)*(1-exp(-ptimestep/1.e6))/ptimestep
	1459	! pdq(ig,l,igcm_eddy1e7)=pdq(ig,l,igcm_eddy1e7)-
	1460	! & pq(ig,l,igcm_eddy1e7)*(1-exp(-ptimestep/1.e7))/ptimestep
	1461	! pdq(ig,l,igcm_eddy5e7)=pdq(ig,l,igcm_eddy5e7)-
	1462	! & pq(ig,l,igcm_eddy5e7)*(1-exp(-ptimestep/5.e7))/ptimestep
	1463	! pdq(ig,l,igcm_eddy1e8)=pdq(ig,l,igcm_eddy1e8)-
	1464	! & pq(ig,l,igcm_eddy1e8)*(1-exp(-ptimestep/1.e8))/ptimestep
	1465	! pdq(ig,l,igcm_eddy5e8)=pdq(ig,l,igcm_eddy5e8)-
	1466	! & pq(ig,l,igcm_eddy5e8)*(1-exp(-ptimestep/5.e8))/ptimestep
	1467	! ENDDO
	1468	! ENDDO
	1469
	1470	! -------------------------------
	1471	! Update tracers
	1472	! -------------------------------
	1473	DO iq=1, nq
	1474	DO l=1,nlayer
	1475	DO ig=1,ngrid
	1476	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
	1477	ENDDO
	1478	ENDDO
	1479	ENDDO
	1480
	1481	DO ig=1,ngrid
	1482	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
	1483	ENDDO
	1484	call globalaverage2d(ngrid,ps,globave)
	1485
	1486	! pressure density
	1487	IF (.not.fast) then !
	1488	DO ig=1,ngrid
	1489	DO l=1,nlayer
	1490	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
	1491	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
	1492	rho(ig,l) = zplay(ig,l)/(r*zt(ig,l))
	1493	ENDDO
	1494	zplev(ig,nlayer+1)=pplev(ig,nlayer+1)/pplev(ig,1)*ps(ig)
	1495	ENDDO
	1496	ENDIF
	1497
	1498	! ---------------------------------------------------------
	1499	! Compute column amounts (kg m-2) if tracers are enabled
	1500	! ---------------------------------------------------------
	1501	call zerophys(ngrid*nq,qcol)
	1502	if(tracer.and..not.fast)then
	1503	do iq=1,nq
	1504	DO ig=1,ngrid
	1505	DO l=1,nlayer
	1506	qcol(ig,iq) = qcol(ig,iq) + zq(ig,l,iq) *
	1507	& (zplev(ig,l) - zplev(ig,l+1)) / g
	1508	enddo
	1509	enddo
	1510	enddo
	1511	endif
	1512
	1513	if (methane) then
	1514	IF (fast) then ! zq is the mixing ratio supposingly mixed in all atmosphere
	1515	DO ig=1,ngrid
	1516	vmr_ch4(ig)=zq(ig,1,igcm_ch4_gas)*
	1517	& mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.
	1518	ENDDO
	1519	ELSE
	1520	DO ig=1,ngrid
	1521	! compute vmr methane
	1522	vmr_ch4(ig)=qcol(ig,igcm_ch4_gas)*
	1523	& g/ps(ig)mmol(igcm_n2)/mmol(igcm_ch4_gas)100.
	1524	! compute density methane
	1525	DO l=1,nlayer
	1526	zrho_ch4(ig,l)=zq(ig,l,igcm_ch4_gas)*rho(ig,l)
	1527	ENDDO
	1528	ENDDO
	1529	ENDIF
	1530	endif
	1531
	1532	if (carbox) then
	1533	IF (fast) then
	1534	DO ig=1,ngrid
	1535	vmr_co(ig)=zq(ig,1,igcm_co_gas)*
	1536	& mmol(igcm_n2)/mmol(igcm_co_gas)*100.
	1537	ENDDO
	1538	ELSE
	1539	DO ig=1,ngrid
	1540	! compute vmr CO
	1541	vmr_co(ig)=qcol(ig,igcm_co_gas)*
	1542	& g/ps(ig)mmol(igcm_n2)/mmol(igcm_co_gas)100.
	1543	! compute density CO
	1544	DO l=1,nlayer
	1545	zrho_co(ig,l)=zq(ig,l,igcm_co_gas)*rho(ig,l)
	1546	ENDDO
	1547	ENDDO
	1548	ENDIF
	1549	endif
	1550
	1551	zrho_haze(:,:)=0.
	1552	zdqrho_photprec(:,:)=0.
	1553	IF (haze.and.aerohaze) then
	1554	DO ig=1,ngrid
	1555	DO l=1,nlayer
	1556	zrho_haze(ig,l)=zq(ig,l,igcm_haze)*rho(ig,l)
	1557	zdqrho_photprec(ig,l)=zdqphot_prec(ig,l)*rho(ig,l)
	1558	ENDDO
	1559	ENDDO
	1560	ENDIF
	1561
	1562	IF (fasthaze) then
	1563	DO ig=1,ngrid
	1564	qcol(ig,igcm_haze)=zq(ig,1,igcm_haze)*pplev(ig,1)/g
	1565	qcol(ig,igcm_prec_haze)=zq(ig,1,igcm_prec_haze)*pplev(ig,1)/g
	1566	ENDDO
	1567	ENDIF
	1568
	1569	c Info about Ls, declin...
	1570	IF (fast) THEN
	1571	write(,),'Ls=',zls180./pi,' dec=',declin180./pi
	1572	write(,),'zday=',zday,' ps=',globave
	1573	IF (lastcall) then
	1574	write(,),'lastcall'
	1575	ENDIF
	1576	ELSE
	1577	write(,),'Ls=',zls180./pi,'decli=',declin180./pi,'zday=',zday
	1578	ENDIF
	1579
	1580	lecttsoil=0 ! default value for lecttsoil
	1581	call getin("lecttsoil",lecttsoil)
	1582	IF (lastcall.and.(ngrid.EQ.1).and.(lecttsoil.eq.1)) THEN
	1583	! save tsoil temperature profile for 1D profile
	1584	OPEN(13,file='proftsoil.out',form='formatted')
	1585	DO i=1,nsoilmx
	1586	write(13,*) tsoil(1,i)
	1587	ENDDO
	1588	CLOSE(13)
	1589	ENDIF
	1590
	1591
	1592	IF (ngrid.NE.1) THEN
	1593	! -------------------------------------------------------------------
	1594	! Writing NetCDF file "RESTARTFI" at the end of the run
	1595	! -------------------------------------------------------------------
	1596	! Note: 'restartfi' is stored just before dynamics are stored
	1597	! in 'restart'. Between now and the writting of 'restart',
	1598	! there will have been the itau=itau+1 instruction and
	1599	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	1600	! thus we store for time=time+dtvr
	1601
	1602	IF(lastcall) THEN
	1603
	1604	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
	1605	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
	1606
	1607	if (paleo) then
	1608	! time range for tendencies of ice flux qsurfyear
	1609	zdt_tot=year_day ! Last year of simulation
	1610
	1611	masslost(:)=0.
	1612	massacc(:)=0.
	1613	DO ig=2,ngrid-1
	1614	realarea(ig)=area(ig)
	1615	ENDDO
	1616	realarea(1)=area(1)*iim
	1617	realarea(ngrid)=area(ngrid)*iim
	1618
	1619
	1620	DO ig=1,ngrid
	1621	! update new reservoir of ice on the surface
	1622	DO iq=1,nq
	1623	! kg/m2 to be sublimed or condensed during paleoyears
	1624	qsurfyear(ig,iq)=qsurfyear(ig,iq)*
	1625	& paleoyears365.25/(zdt_totdaysec/86400.)
	1626
	1627	! special case if we sublime the entire reservoir
	1628	IF (-qsurfyear(ig,iq).gt.qsurf(ig,iq)) THEN
	1629	masslost(iq)=masslost(iq)+realarea(ig)*
	1630	& (-qsurfyear(ig,iq)-qsurf(ig,iq))
	1631	qsurfyear(ig,iq)=-qsurf(ig,iq)
	1632	ENDIF
	1633
	1634	IF (qsurfyear(ig,iq).gt.0.) THEN
	1635	massacc(iq)=massacc(iq)+realarea(ig)*qsurfyear(ig,iq)
	1636	ENDIF
	1637
	1638	ENDDO
	1639	ENDDO
	1640
	1641	DO ig=1,ngrid
	1642	DO iq=1,nq
	1643	qsurfpal(ig,iq)=qsurf(ig,iq)+qsurfyear(ig,iq)
	1644	IF (qsurfyear(ig,iq).gt.0.) THEN
	1645	qsurfpal(ig,iq)=qsurfpal(ig,iq)-
	1646	& qsurfyear(ig,iq)*masslost(iq)/massacc(iq)
	1647	ENDIF
	1648	ENDDO
	1649	ENDDO
	1650	! Finally ensure conservation of qsurf
	1651	DO iq=1,nq
	1652	call globalaverage2d(ngrid,qsurf(:,iq),globaveice(iq))
	1653	call globalaverage2d(ngrid,qsurfpal(:,iq),
	1654	& globavenewice(iq))
	1655	IF (globavenewice(iq).gt.0.) THEN
	1656	qsurfpal(:,iq)=qsurfpal(:,iq)*
	1657	& globaveice(iq)/globavenewice(iq)
	1658	ENDIF
	1659	ENDDO
	1660
	1661	! update new geopotential depending on the ice reservoir
	1662	phisfipal(:)=phisfinew(:)+qsurfpal(:,igcm_n2)*g/1000.
	1663	!phisfipal(ig)=phisfi(ig)
	1664
	1665
	1666	if (kbo.or.triton) then ! case of Triton : we do not change the orbital parameters
	1667
	1668	pdaypal=pday ! no increment of pdaypal to keep same evolution of the subsolar point
	1669	eccpal=1.-periheli/((periheli+aphelie)/2.) !no change of ecc
	1670	peri_daypal=peri_day ! no change
	1671	oblipal=obliquit ! no change
	1672	tpalnew=tpal
	1673	adjustnew=adjust
	1674
	1675	else ! Pluto
	1676	! update new pday and tpal (Myr) to be set in startfi controle
	1677	pdaypal=int(day_ini+paleoyears*365.25/6.3872)
	1678	tpalnew=tpal+paleoyears*1.e-6 ! Myrs
	1679
	1680	! update new N2 ice adjustment (not tested yet on Pluto)
	1681	adjustnew=adjust
	1682
	1683	! update new obliquity
	1684	oblipal=23./2.*
	1685	& sin(2.pi/2.77(tpalnew+0.16))+115.5
	1686
	1687	! update new peri_day
	1688	call calcperiday(tpalnew,peri_daypal)
	1689	!peri_daypal=peri_day
	1690
	1691	! update new eccentricity
	1692	call calcecc(tpalnew,eccpal)
	1693	!eccpal=0.009
	1694
	1695	endif
	1696
	1697	write(,) "Paleo peri=",peri_daypal," tpal=",tpalnew
	1698	write(,) "Paleo eccpal=",eccpal," tpal=",tpalnew
	1699
	1700	! create restartfi
	1701	call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq,
	1702	. ptimestep,pdaypal,
	1703	. ztime_fin,tsurf,tsoil,emis,q2,qsurfpal,
	1704	. area,albedodat,tidat,zmea,zstd,zsig,
	1705	. zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal,
	1706	. peri_daypal)
	1707	else
	1708
	1709	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
	1710	. ptimestep,pday,
	1711	. ztime_fin,tsurf,tsoil,emis,q2,qsurf,
	1712	. area,albedodat,tidat,zmea,zstd,zsig,
	1713	. zgam,zthe)
	1714
	1715	endif
	1716
	1717	ENDIF
	1718
	1719	! -----------------------------------------------------------------
	1720	! Saving statistics :
	1721	! -----------------------------------------------------------------
	1722	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
	1723	! which can later be used to make the statistic files of the run:
	1724	! "stats") only possible in 3D runs !
	1725
	1726
	1727	IF (callstats) THEN
	1728	write (,) "stats have been turned off in the code.
	1729	& You can manage your own output in physiq.F "
	1730	stop
	1731
	1732	! call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
	1733	! call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1734	! call wstats(ngrid,"fluxsurf_lw",
	1735	! . "Thermal IR radiative flux to surface","W.m-2",2,
	1736	! . fluxsurf_lw)
	1737	! call wstats(ngrid,"fluxsurf_sw",
	1738	! . "Solar radiative flux to surface","W.m-2",2,
	1739	! . fluxsurf_sw_tot)
	1740	! call wstats(ngrid,"fluxtop_lw",
	1741	! . "Thermal IR radiative flux to space","W.m-2",2,
	1742	! . fluxtop_lw)
	1743	! call wstats(ngrid,"fluxtop_sw",
	1744	! . "Solar radiative flux to space","W.m-2",2,
	1745	! . fluxtop_sw_tot)
	1746
	1747	! call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
	1748	! call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
	1749	! call wstats(ngrid,"v","Meridional (North-South) wind",
	1750	! . "m.s-1",3,zv)
	1751
	1752	IF(lastcall) THEN
	1753	write (,) "Writing stats..."
	1754	call mkstats(ierr)
	1755	ENDIF
	1756	ENDIF !if callstats
	1757
	1758
	1759	! ---------------------------------------------------------------------
	1760	! 3D OUTPUTS
	1761	! ----------------------------------------------------------------------
	1762	! output in netcdf file "DIAGFI", containing any variable for diagnostic
	1763	! (output with period "ecritphy", set in "run.def")
	1764	! ----------------------------------------------------------------------
	1765
	1766	! if(MOD(countG3D,saveG3D).eq.0)then
	1767	! print*,'countG3D',countG3D
	1768
	1769	!! BASIC
	1770
	1771	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
	1772	& tsurf)
	1773	call WRITEDIAGFI(ngrid,"rice_ch4","ch4 ice mass mean radius",
	1774	& "m",3,rice_ch4)
	1775	call WRITEDIAGFI(ngrid,"area","area","m",2,area)
	1776	call WRITEDIAGFI(ngrid,"mu0","cos sza","m",2,mu0)
	1777	call WRITEDIAGFI(ngrid,"fract","fractional day","",2,fract)
	1778	call WRITEDIAGFI(ngrid,"declin","declin","deg",0,declin*180./pi)
	1779	call WRITEDIAGFI(ngrid,"ls","ls","deg",0,zls*180./pi)
	1780	call WRITEDIAGFI(ngrid,"dist_sol","dist_sol","AU",0,dist_sol)
	1781	call WRITEDIAGFI(ngrid,"phisfi","phisfi"," ",2,phisfi)
	1782	! call WRITEDIAGFI(ngrid,"tsub","subsurface temperature","K",
	1783	! & 1,tsub)
	1784	! call WRITEDIAGFI(ngrid,"phitop","phitop"," ",2,phitop)
	1785
	1786	! ENERGY - Total energy balance diagnostics
	1787
	1788	call WRITEDIAGFI(ngrid,"albedo","albedo","sans dim",2,
	1789	& albedo)
	1790	call WRITEDIAGFI(ngrid,"emissivite","emissivite","sans dim",2,
	1791	& emis)
	1792	call WRITEDIAGFI(ngrid,"fluxtop_dn","fluxtop_dn","sans dim",2,
	1793	& fluxtop_dn)
	1794	call WRITEDIAGFI(ngrid,"ISR","incoming stellar rad.","W m-2",
	1795	& 2,fluxtop_dn)
	1796	call WRITEDIAGFI(ngrid,"ASR","absorbed stellar rad.","W m-2",
	1797	& 2,fluxabs_sw)
	1798	call WRITEDIAGFI(ngrid,"OLR","outgoing longwave rad.","W m-2",
	1799	& 2,fluxtop_lw)
	1800
	1801	IF (fast) then
	1802	! pressure reprocessed in nogcm.F
	1803	call WRITEDIAGFI(ngrid,"globave","surf press","Pa",0,globave)
	1804	call WRITEDIAGFI(ngrid,"ps","surf press","Pa",2,ps)
	1805	call WRITEDIAGFI(ngrid,"fluxrad","fluxrad",
	1806	& "W m-2",2,fluxrad)
	1807	call WRITEDIAGFI(ngrid,"fluxgrd","fluxgrd",
	1808	& "W m-2",2,fluxgrd)
	1809	call WRITEDIAGFI(ngrid,"capcal","capcal",
	1810	& "W.s m-2 K-1",2,capcal)
	1811	! call WRITEDIAGFI(ngrid,"dplanck","dplanck",
	1812	! & "W.s m-2 K-1",2,dplanck)
	1813	call WRITEDIAGFI(ngrid,"tsoil","tsoil","K",3,tsoil)
	1814
	1815	! If we want to see tidat evolution with time:
	1816	call WRITEDIAGFI(ngrid,"tidat","tidat","SI",3,tidat_out)
	1817
	1818
	1819	IF (fasthaze) then
	1820	call WRITEDIAGFI(ngrid,"gradflux","gradflux",
	1821	& "Ph m-2 s-1",2,gradflux)
	1822	call WRITEDIAGFI(ngrid,"fluxbot","flux bottom",
	1823	& "Ph m-2 s-1",2,fluxbot)
	1824	call WRITEDIAGFI(ngrid,"fluxbotsol","flux bottom sol",
	1825	& "Ph m-2 s-1",2,fluxlym_sol_bot)
	1826	call WRITEDIAGFI(ngrid,"fluxbotipm","flux bottom ipm",
	1827	& "Ph m-2 s-1",2,fluxlym_ipm_bot)
	1828	call WRITEDIAGFI(ngrid,"fluxlymipm","flux incid ipm",
	1829	& "Ph m-2 s-1",2,flym_ipm)
	1830	call WRITEDIAGFI(ngrid,"fluxlymsol","flux incid sol",
	1831	& "Ph m-2 s-1",2,flym_sol)
	1832	call WRITEDIAGFI(ngrid,"tend_prodhaze","prod haze",
	1833	& "kg m-2 s-1",2,zdqsprodhaze)
	1834	call WRITEDIAGFI(ngrid,"tend_lostch4","tend ch4",
	1835	& "kg kg-1 s-1",2,zdqprodhaze(1,igcm_ch4_gas))
	1836	call WRITEDIAGFI(ngrid,"tend_prodhaze","tend prod haze",
	1837	& "kg kg-1 s-1",2,zdqprodhaze(1,igcm_haze))
	1838	ENDIF
	1839	IF (paleo) then
	1840	call WRITEDIAGFI(ngrid,"qsurfn2_year","qsurfn2_year",
	1841	& "kg m-2.plutoyear-1",2,qsurfyear(:,1))
	1842	call WRITEDIAGFI(ngrid,"phisfipal","phisfipal",
	1843	& " ",2,phisfipal)
	1844	call WRITEDIAGFI(ngrid,"zdqflow","zdqflow",
	1845	& " ",2,zdqflow(1,igcm_n2))
	1846	ENDIF
	1847	ELSE
	1848	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
	1849	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
	1850	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
	1851	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
	1852	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1853	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
	1854	! call WRITEDIAGFI(ngrid,"pressure","Pression","Pa",3,pplay)
	1855	call WRITEDIAGFI(ngrid,"fluxrad","fluxrad",
	1856	& "W m-2",2,fluxrad)
	1857	call WRITEDIAGFI(ngrid,"fluxgrd","fluxgrd",
	1858	& "W m-2",2,fluxgrd)
	1859	ENDIF
	1860
	1861	! Usually not used :
	1862
	1863	! call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
	1864	call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
	1865	! call WRITEDIAGFI(ngrid,"DYN","dynamical heat input","W m-2",
	1866	! & 2,fluxdyn)
	1867
	1868	! TENDANCIES
	1869
	1870	call WRITEDIAGFI(ngrid,"dps","dps","K",2,pdpsrf)
	1871	call WRITEDIAGFI(ngrid,"zdtsw","SW heating","K s-1",3,zdtsw)
	1872	call WRITEDIAGFI(ngrid,"zdtlw","LW heating","K s-1",3,zdtlw)
	1873	call WRITEDIAGFI(ngrid,"zdtc","tendancy T cond N2","K",3,zdtc)
	1874	call WRITEDIAGFI(ngrid,"zdtch4cloud","tendancy T ch4cloud",
	1875	& "K",3,zdtch4cloud)
	1876	call WRITEDIAGFI(ngrid,"zdtcocloud","tendancy T cocloud",
	1877	& "K",3,zdtcocloud)
	1878	call WRITEDIAGFI(ngrid,"zdtsurfc","zdtsurfc","K",2,zdtsurfc)
	1879	call WRITEDIAGFI(ngrid,"zdtdif","zdtdif","K",3,zdtdif)
	1880	call WRITEDIAGFI(ngrid,"zdtconduc","tendancy T conduc",
	1881	& "K",3,zdtconduc)
	1882	call WRITEDIAGFI(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
	1883	call WRITEDIAGFI(ngrid,"zdtrad","rad heating","T s-1",3,dtrad)
	1884	call WRITEDIAGFI(ngrid,"zdtadj","conv adj","T s-1",3,zdtadj)
	1885	call WRITEDIAGFI(ngrid,"zdqcloud","ch4 cloud","T s-1",
	1886	& 3,zdqch4cloud(1,1,igcm_ch4_ice))
	1887	call WRITEDIAGFI(ngrid,"zdqcloudgas","ch4 cloud","T s-1",
	1888	& 3,zdqch4cloud(1,1,igcm_ch4_gas))
	1889	call WRITEDIAGFI(ngrid,"zdqcn2_c","zdq condn2","",3,zdqc(:,:,3))
	1890	call WRITEDIAGFI(ngrid,"zdqdif_c","zdqdif","",3,zdqdif(:,:,3))
	1891	call WRITEDIAGFI(ngrid,"zdqadj_c","zdqadj","",3,zdqadj(:,:,3))
	1892	!call WRITEDIAGFI(ngrid,"zdqsed_h","zdqsed","",3,zdqsed(:,:,7))
	1893	!call WRITEDIAGFI(ngrid,"zdqssed","zdqssed","",2,zdqssed)
	1894	call WRITEDIAGFI(ngrid,"zq1temp_ch4"," "," ",2,zq1temp_ch4)
	1895	call WRITEDIAGFI(ngrid,"qsat_ch4"," "," ",2,qsat_ch4)
	1896	call WRITEDIAGFI(ngrid,"qsat_ch4_l1"," "," ",2,qsat_ch4_l1)
	1897	call WRITEDIAGFI(ngrid,"senshf1","senshf1"," ",2,sensiblehf1)
	1898	call WRITEDIAGFI(ngrid,"senshf2","senshf2"," ",2,sensiblehf2)
	1899
	1900
	1901	! OTHERS
	1902
	1903	! call WRITEDIAGFI(ngrid,"dWtotdifv","dWtotdifv","kg m-2",1,dWtot)
	1904	! call WRITEDIAGFI(ngrid,"dWtotsdifv","dWtotsdifv","kgm-2",1,dWtots)
	1905	! call WRITEDIAGFI(ngrid,"nconsMX","nconsMX","kgm-2s-1",1,nconsMAX)
	1906
	1907
	1908	! TRACERS
	1909
	1910	if (tracer) then
	1911
	1912	!!! afficher la tendance sur le flux de la glace
	1913	call WRITEDIAGFI(ngridmx,'n2_iceflux',
	1914	& 'n2_iceflux',
	1915	& "kg m^-2 s^-1",2,flusurf(1,igcm_n2) )
	1916
	1917	do iq=1,nq
	1918	call WRITEDIAGFI(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
	1919	call WRITEDIAGFI(ngrid,trim(noms(iq))//'_Tcondn2',
	1920	& trim(noms(iq))//'_Tcondn2',
	1921	& 'kg/kg',3,zdqc(1,1,iq))
	1922
	1923	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_surf',
	1924	& trim(noms(iq))//'_surf',
	1925	& 'kg m^-2',2,qsurf(1,iq) )
	1926
	1927	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_col',
	1928	& trim(noms(iq))//'_col',
	1929	& 'kg m^-2',2,qcol(1,iq) )
	1930
	1931	! call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_aero',
	1932	! & trim(noms(iq))//'_aero',
	1933	! & 'kg/kg',3,aerosol(1,1,iq))
	1934
	1935
	1936	enddo
	1937
	1938	call WRITEDIAGFI(ngridmx,'haze_reff',
	1939	& 'haze_reff',
	1940	& 'm',3,reffrad(1,1,1))
	1941
	1942	if (methane) then
	1943
	1944	call WRITEDIAGFI(ngridmx,'ch4_iceflux',
	1945	& 'ch4_iceflux',
	1946	& "kg m^-2 s^-1",2,flusurf(1,igcm_ch4_ice) )
	1947
	1948	call WRITEDIAGFI(ngrid,"vmr_ch4","vmr_ch4","%",
	1949	& 2,vmr_ch4)
	1950	if (.not.fast) THEN
	1951	call WRITEDIAGFI(ngrid,"zrho_ch4","zrho_ch4","kg.m-3",
	1952	& 3,zrho_ch4(:,:))
	1953	endif
	1954
	1955	! if (fast) THEN
	1956	! call WRITEDIAGFI(ngrid,"zq1_ch4","zq ch4","kg m-2",
	1957	! & 2,zq(1,1,igcm_ch4_gas))
	1958	! endif
	1959
	1960	! Tendancies
	1961	call WRITEDIAGFI(ngrid,"zdqch4cloud","ch4 cloud","T s-1",
	1962	& 3,zdqch4cloud(1,1,igcm_ch4_gas))
	1963	call WRITEDIAGFI(ngrid,"zdqcn2_ch4","zdq condn2 ch4","",
	1964	& 3,zdqc(:,:,igcm_ch4_gas))
	1965	call WRITEDIAGFI(ngrid,"zdqdif_ch4","zdqdif ch4","",
	1966	& 3,zdqdif(:,:,igcm_ch4_gas))
	1967	call WRITEDIAGFI(ngrid,"zdqsdif_ch4_ice","zdqsdif ch4","",
	1968	& 2,zdqsdif(:,igcm_ch4_ice))
	1969	call WRITEDIAGFI(ngrid,"zdqadj_ch4","zdqadj ch4","",
	1970	& 3,zdqadj(:,:,igcm_ch4_gas))
	1971	call WRITEDIAGFI(ngrid,"zdqsed_ch4","zdqsed ch4","",
	1972	& 3,zdqsed(:,:,igcm_ch4_gas))
	1973	call WRITEDIAGFI(ngrid,"zdqssed_ch4","zdqssed ch4","",
	1974	& 2,zdqssed(:,igcm_ch4_gas))
	1975	call WRITEDIAGFI(ngrid,"zdtch4cloud","ch4 cloud","T s-1",
	1976	& 3,zdtch4cloud)
	1977
	1978	endif
	1979
	1980	if (carbox) then
	1981
	1982	call WRITEDIAGFI(ngridmx,'co_iceflux',
	1983	& 'co_iceflux',
	1984	& "kg m^-2 s^-1",2,flusurf(1,igcm_co_ice) )
	1985
	1986	call WRITEDIAGFI(ngrid,"vmr_co","vmr_co","%",
	1987	& 2,vmr_co)
	1988
	1989	if (.not.fast) THEN
	1990	call WRITEDIAGFI(ngrid,"zrho_co","zrho_co","kg.m-3",
	1991	& 3,zrho_co(:,:))
	1992	endif
	1993
	1994	! if (fast) THEN
	1995	! call WRITEDIAGFI(ngrid,"zq1_co","zq co","kg m-2",
	1996	! & 2,zq(1,1,igcm_co_gas))
	1997	! endif
	1998
	1999	! Tendancies
	2000	! call WRITEDIAGFI(ngrid,"zdqcocloud","co cloud","T s-1",
	2001	! & 3,zdqcocloud(1,1,igcm_co_gas))
	2002	! call WRITEDIAGFI(ngrid,"zdqcn2_co","zdq condn2 co","",
	2003	! & 3,zdqc(:,:,igcm_co_gas))
	2004	! call WRITEDIAGFI(ngrid,"zdqdif_co","zdqdif co","",
	2005	! & 3,zdqdif(:,:,igcm_co_gas))
	2006	! call WRITEDIAGFI(ngrid,"zdqsdif_co_ice","zdqsdif co","",
	2007	! & 2,zdqsdif(:,igcm_co_ice))
	2008	! call WRITEDIAGFI(ngrid,"zdqadj_co","zdqadj co","",
	2009	! & 3,zdqadj(:,:,igcm_co_gas))
	2010	! call WRITEDIAGFI(ngrid,"zdqsed_co","zdqsed co","",
	2011	! & 3,zdqsed(:,:,igcm_co_gas))
	2012	! call WRITEDIAGFI(ngrid,"zdqssed_co","zdqssed co","",
	2013	! & 2,zdqssed(:,igcm_co_gas))
	2014	! call WRITEDIAGFI(ngrid,"zdtcocloud","co cloud","T s-1",
	2015	! & 3,zdtcocloud)
	2016
	2017	endif
	2018
	2019	if (haze) then
	2020
	2021	! call WRITEDIAGFI(ngrid,"zrho_haze","zrho_haze","kg.m-3",
	2022	! & 3,zrho_haze(:,:))
	2023	call WRITEDIAGFI(ngrid,"zdqrho_photprec","zdqrho_photprec",
	2024	& "kg.m-3.s-1",3,zdqrho_photprec(:,:))
	2025	call WRITEDIAGFI(ngrid,"zdqphot_prec","zdqphot_prec","",
	2026	& 3,zdqphot_prec(:,:))
	2027	call WRITEDIAGFI(ngrid,"zdqhaze_ch4","zdqhaze_ch4","",
	2028	& 3,zdqhaze(:,:,igcm_ch4_gas))
	2029	call WRITEDIAGFI(ngrid,"zdqhaze_prec","zdqhaze_prec","",
	2030	& 3,zdqhaze(:,:,igcm_prec_haze))
	2031	call WRITEDIAGFI(ngrid,"zdqhaze_haze","zdqhaze_haze","",
	2032	& 3,zdqhaze(:,:,igcm_haze))
	2033	call WRITEDIAGFI(ngrid,"zdqphot_ch4","zdqphot_ch4","",
	2034	& 3,zdqphot_ch4(:,:))
	2035	call WRITEDIAGFI(ngrid,"zdqconv_prec","zdqconv_prec","",
	2036	& 3,zdqconv_prec(:,:))
	2037	call WRITEDIAGFI(ngrid,"zdqssed_haze","zdqssed haze",
	2038	& "kg/m2/s",2,zdqssed(:,igcm_haze))
	2039	! call WRITEDIAGFI(ngrid,"zdqhaze_col","zdqhaze col","kg/m2/s",
	2040	! & 2,zdqhaze_col(:))
	2041	endif
	2042
	2043	if (aerohaze) then
	2044	call WRITEDIAGFI(ngridmx,"tau_col",
	2045	& "Total aerosol optical depth","opacity",2,tau_col)
	2046	endif
	2047
	2048	! call WRITEDIAGFI(ngridmx,"tau_col",
	2049	! & "Total aerosol optical depth","[]",2,tau_col)
	2050	endif
	2051
	2052	! countG3D=1
	2053	! else
	2054	! countG3D=countG3D+1
	2055	! endif ! if time to save
	2056
	2057	! ----------------------------------------------------------------------
	2058	! 1D OUTPUTS
	2059	! ----------------------------------------------------------------------
	2060	ELSE ! if(ngrid.eq.1)
	2061
	2062	if(countG1D.eq.saveG1D)then
	2063
	2064	! BASIC 1D ONLY
	2065
	2066	call WRITEDIAGFI(ngrid,"ISR","incoming stellar rad.","W m-2",
	2067	& 0,fluxtop_dn)
	2068	call WRITEDIAGFI(ngrid,"ASR","absorbed stellar rad.","W m-2",
	2069	& 0,fluxabs_sw)
	2070	call WRITEDIAGFI(ngrid,"OLR","outgoing longwave rad.","W m-2",
	2071	& 0,fluxtop_lw)
	2072	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
	2073	& tsurf)
	2074	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",0,ps)
	2075	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
	2076
	2077	call WRITEDIAGFI(ngrid,"fluxsurf_sw","sw surface flux",
	2078	& "K",0,fluxsurf_sw)
	2079	call WRITEDIAGFI(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
	2080	call WRITEDIAGFI(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
	2081	call WRITEDIAGFI(ngrid,"zdtdif","zdtdif","K",3,zdtdif)
	2082	call WRITEDIAGFI(ngrid,"zdtconduc","tendancy T conduc",
	2083	& "K",3,zdtconduc)
	2084
	2085	call WRITEDIAGFI(ngrid,"fluxsurf_sw","sw surface flux",
	2086	& "K",0,fluxsurf_sw)
	2087	call WRITEDIAGFI(ngrid,"declin","lat subsolaire",
	2088	& "deg",0,declin*180./pi)
	2089	call WRITEDIAGFI(ngrid,"ls","ls","deg",0,zls*180./pi)
	2090	call WRITEDIAGFI(ngrid,"dist_sol","dist_sol","AU",0,dist_sol)
	2091	call WRITEDIAGFI(ngrid,"mu0","solar zenith angle",
	2092	& "deg",0,mu0)
	2093	call WRITEDIAGFI(ngrid,"albedo","albedo",
	2094	& "",0,albedo)
	2095	call WRITEDIAGFI(ngrid,"emis","emis",
	2096	& "",0,emis)
	2097	call WRITEDIAGFI(ngrid,"tsoil","tsoil","K",1,tsoil(1,:))
	2098	call WRITEDIAGFI(ngrid,"tidat","tidat","SI",1,tidat_out(1,:))
	2099
	2100	! OUTPUT OF TENDANCIES
	2101	! call WRITEDIAGFI(ngrid,"zdqcloud","ch4 cloud","T s-1",
	2102	! & 3,zdqcloud(1,1,igcm_ch4_gas))
	2103	! call WRITEDIAGFI(ngrid,"zdqcn2","zdq condn2","",1,zdqc(:,:,3))
	2104	! call WRITEDIAGFI(ngrid,"zdqdif","zdqdif","",1,zdqdif(:,:,3))
	2105	! call WRITEDIAGFI(ngrid,"zdqadj","zdqadj","",1,zdqadj(:,:,3))
	2106	! call WRITEDIAGFI(ngrid,"zdqsed","zdqsed","",1,zdqsed(:,:,3))
	2107	! call WRITEDIAGFI(ngrid,"zdqc","zdqc","",1,zdqc(:,:,3))
	2108	! call WRITEDIAGFI(ngrid,"zdqhaze","zdqhaze","",1,
	2109	! & zdqhaze(:,:,3))
	2110
	2111	! 1D methane
	2112	if (methane) then
	2113
	2114	call WRITEDIAGFI(ngridmx,'ch4_iceflux',
	2115	& 'ch4_iceflux',
	2116	& "kg m^-2 s^-1",0,flusurf(1,igcm_ch4_ice) )
	2117
	2118	call WRITEDIAGFI(ngrid,"vmr_ch4","vmr_ch4","%",
	2119	& 0,vmr_ch4(1))
	2120	call WRITEDIAGFI(ngrid,"zrho_ch4","zrho_ch4","kg.m-3",
	2121	& 1,zrho_ch4(1,:))
	2122
	2123	! Tendancies
	2124	call WRITEDIAGFI(ngrid,"zdqch4cloud","ch4 cloud","T s-1",
	2125	& 1,zdqch4cloud(1,1,igcm_ch4_gas))
	2126	call WRITEDIAGFI(ngrid,"zdqcn2_ch4","zdq condn2 ch4","",
	2127	& 1,zdqc(1,:,igcm_ch4_gas))
	2128	call WRITEDIAGFI(ngrid,"zdqdif_ch4","zdqdif ch4","",
	2129	& 1,zdqdif(1,:,igcm_ch4_gas))
	2130	call WRITEDIAGFI(ngrid,"zdqsdif_ch4_ice","zdqsdif ch4","",
	2131	& 1,zdqsdif(:,igcm_ch4_ice))
	2132	call WRITEDIAGFI(ngrid,"zdqadj_ch4","zdqadj ch4","",
	2133	& 1,zdqadj(1,:,igcm_ch4_gas))
	2134	call WRITEDIAGFI(ngrid,"zdqsed_ch4","zdqsed ch4","",
	2135	& 1,zdqsed(1,:,igcm_ch4_gas))
	2136	call WRITEDIAGFI(ngrid,"zdqssed_ch4","zdqssed ch4","",
	2137	& 0,zdqssed(1,igcm_ch4_gas))
	2138
	2139	endif
	2140
	2141	! 1D Haze
	2142	if (haze) then
	2143
	2144	! call WRITEDIAGFI(ngrid,"zrho_haze","zrho_haze","kg.m-3",
	2145	! & 1,zrho_haze(:,:))
	2146	call WRITEDIAGFI(ngrid,"zdqrho_photprec","zdqrho_photprec",
	2147	& "kg.m-3.s-1",3,zdqrho_photprec)
	2148	call WRITEDIAGFI(ngrid,"zdqphot_prec","zdqphot_prec","",
	2149	& 3,zdqphot_prec)
	2150	call WRITEDIAGFI(ngrid,"zdqhaze_ch4","zdqhaze_ch4","",
	2151	& 1,zdqhaze(1,:,igcm_ch4_gas))
	2152	call WRITEDIAGFI(ngrid,"zdqhaze_prec","zdqhaze_prec","",
	2153	& 1,zdqhaze(1,:,igcm_prec_haze))
	2154	call WRITEDIAGFI(ngrid,"zdqhaze_haze","zdqhaze_haze","",
	2155	& 1,zdqhaze(1,:,igcm_haze))
	2156	call WRITEDIAGFI(ngrid,"zdqphot_ch4","zdqphot_ch4","",
	2157	& 3,zdqphot_ch4)
	2158	! call WRITEDIAGFI(ngrid,"zdqconv_prec","zdqconv_prec","",
	2159	! & 1,zdqconv_prec(1,:))
	2160	call WRITEDIAGFI(ngrid,"zdqssed_haze","zdqssed haze","kg/m2/s",
	2161	& 0,zdqssed(1,igcm_haze))
	2162	! call WRITEDIAGFI(ngrid,"zdqhaze_col","zdqhaze col","kg/m2/s",
	2163	! & 0,zdqhaze_col(:))
	2164	if (haze_radproffix) then
	2165	call WRITEDIAGFI(ngridmx,'haze_reffrad','haze_reffrad','m',
	2166	& 3,reffrad(1,1,1))
	2167	endif
	2168	endif
	2169
	2170	if (aerohaze) then
	2171	call WRITEDIAGFI(ngridmx,"tau_col",
	2172	& "Total aerosol optical depth","[]",0,tau_col(1))
	2173	call WRITEDIAGFI(ngridmx,"tau_aero",
	2174	& "aerosol optical depth","[]",3,aerosol(1,1,1))
	2175	endif
	2176
	2177	! TRACERS 1D ONLY
	2178	if (tracer) then
	2179	do iq=1,nq
	2180	call WRITEDIAGFI(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
	2181	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_surf',
	2182	& trim(noms(iq))//'_surf',
	2183	& 'kg m^-2',2,qsurf(1,iq) )
	2184	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_col',
	2185	& trim(noms(iq))//'_col',
	2186	& 'kg m^-2',2,qcol(1,iq) )
	2187	enddo
	2188	endif
	2189
	2190	countG1D=1
	2191	else
	2192	countG1D=countG1D+1
	2193	endif ! if time to save
	2194
	2195	END IF ! if(ngrid.ne.1)
	2196
	2197	! countG3D=countG3D+1
	2198	icount=icount+1
	2199
	2200	RETURN
	2201	END

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

Download in other formats:

Original Format