Context Navigation

physiq_mod.F90 @ 3608

Last change on this file since 3608 was 3601, checked in by debatzbr, 2 weeks ago
Fixe bug in the microphysical haze production
Property svn:executable set to ``*
File size: 97.8 KB

Rev	Line
[3184]	1	module physiq_mod
[3195]	2
[3184]	3	implicit none
[3195]	4
[3184]	5	contains
[3195]	6
[3184]	7	subroutine physiq(ngrid,nlayer,nq, &
	8	firstcall,lastcall, &
	9	pday,ptime,ptimestep, &
	10	pplev,pplay,pphi, &
	11	pu,pv,pt,pq, &
	12	flxw, &
	13	pdu,pdv,pdt,pdq,pdpsrf)
	14
	15	!!
[3195]	16	use write_field_phy, only: Writefield_phy
[3184]	17	!!
[3195]	18	use ioipsl_getin_p_mod, only: getin_p
[3184]	19	use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind, corrkdir, banddir
	20	use gases_h, only: gnom, gfrac
	21	use radcommon_h, only: sigma, glat, grav, BWNV, WNOI, DWNI, DWNV, WNOV
	22	use suaer_corrk_mod, only: suaer_corrk
[3572]	23	use aerosol_mod, only: i_haze, haze_prof
[3184]	24	use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, &
	25	dryness
	26	use comdiurn_h, only: coslat, sinlat, coslon, sinlon
	27	use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
[3483]	28	use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
[3539]	29	use geometry_mod, only: latitude, longitude, cell_area, &
	30	cell_area_for_lonlat_outputs
[3184]	31	USE comgeomfi_h, only: totarea, totarea_planet
	32	USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
[3237]	33	igcm_n2,igcm_ch4_gas,igcm_ch4_ice,igcm_haze,&
[3421]	34	igcm_co_gas,igcm_co_ice,igcm_prec_haze,lw_n2,lw_ch4,lw_co,&
[3184]	35	alpha_lift, alpha_devil, qextrhor, &
[3585]	36	nesp, is_chim
[3184]	37	use time_phylmdz_mod, only: ecritphy, iphysiq, nday
[3506]	38	use phyetat0_mod, only: phyetat0,tab_cntrl_mod
[3184]	39	use wstats_mod, only: callstats, wstats, mkstats
	40	use phyredem, only: physdem0, physdem1
	41	use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
	42	use mod_phys_lmdz_para, only : is_master
	43	use planete_mod, only: apoastr, periastr, year_day, peri_day, &
[3539]	44	obliquit, z0, adjust, tpal
[3184]	45	use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
	46	use time_phylmdz_mod, only: daysec
[3455]	47	use callkeys_mod, only: albedo_spectral_mode, calladj, calldifv, &
	48	callrad, callsoil, nosurf, &
	49	callconduct,callmolvis,callmoldiff, &
[3572]	50	corrk, diagdtau, &
[3455]	51	diurnal, enertest, fat1au, &
	52	icetstep, intheat, iradia, kastprof, &
	53	lwrite, mass_redistrib, meanOLR, &
	54	fast,fasthaze,haze,metcloud,monoxcloud, &
[3539]	55	n2cond,noseason_day,conservn2,conservch4, &
[3412]	56	convergeps,kbo,triton,paleo,paleoyears,glaflow, &
[3455]	57	carbox, methane,condmetsurf,condcosurf, &
	58	oldplutovdifc,oldplutocorrk,oldplutosedim, &
[3572]	59	optichaze,haze_proffix,haze_radproffix, &
[3455]	60	source_haze, tsurfmax, albmin_ch4, &
[3572]	61	season, sedimentation, &
[3455]	62	specOLR, &
	63	startphy_file, testradtimes, &
	64	tracer, UseTurbDiff, &
[3557]	65	global1d, szangle, &
	66	callmufi
[3184]	67	use check_fields_mod, only: check_physics_fields
	68	use conc_mod, only: rnew, cpnew, ini_conc_mod
	69	use phys_state_var_mod
	70	use callcorrk_mod, only: callcorrk
[3329]	71	use callcorrk_pluto_mod, only: callcorrk_pluto
[3184]	72	use vdifc_mod, only: vdifc
[3237]	73	use vdifc_pluto_mod, only: vdifc_pluto
[3184]	74	use turbdiff_mod, only: turbdiff
	75	use turb_mod, only : q2,sensibFlux,turb_resolved
	76	use mass_redistribution_mod, only: mass_redistribution
	77	use datafile_mod, only: datadir
[3506]	78	USE vertical_layers_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt
[3184]	79	use mod_phys_lmdz_omp_data, ONLY: is_omp_master
[3539]	80	USE mod_grid_phy_lmdz, ONLY: regular_lonlat, grid_type, unstructured
[3557]	81	! Microphysical model (mp2m)
	82	use mp2m_calmufi
	83	use mp2m_diagnostics
[3184]	84
[3195]	85	#ifdef CPP_XIOS
[3184]	86	use xios_output_mod, only: initialize_xios_output, &
	87	update_xios_timestep, &
	88	send_xios_field
	89	use wxios, only: wxios_context_init, xios_context_finalize
	90	#endif
[3544]	91	USE mod_grid_phy_lmdz, ONLY: grid_type,unstructured,regular_lonlat
[3508]	92	use write_output_mod, only: write_output
[3184]	93
	94	implicit none
	95
	96
	97	!==================================================================
[3195]	98	!
[3184]	99	! Purpose
	100	! -------
	101	! Central subroutine for all the physics parameterisations in the
[3539]	102	! Pluto model. Originally adapted from the Mars LMDZ model.
[3184]	103	!
[3539]	104	! The model can be run with 1 (N2) or more tracer transport
[3184]	105	! depending on the value of "tracer" in file "callphys.def".
	106	!
[3557]	107	! It includes:
[3184]	108	!
	109	! I. Initialization :
	110	! I.1 Firstcall initializations.
	111	! I.2 Initialization for every call to physiq.
	112	!
[3539]	113	! II.1 Thermosphere
	114	! II.2 Compute radiative transfer tendencies (longwave and shortwave) :
[3557]	115	! II.2.a Option 1 : Call correlated-k radiative transfer scheme.
	116	! II.2.b Option 2 : Atmosphere has no radiative effect.
[3184]	117	!
[3539]	118	! III. Vertical diffusion (turbulent mixing)
[3184]	119	!
	120	! IV. Convection :
[3539]	121	! IV.a Dry convective adjusment
[3184]	122	!
[3539]	123	! V. Condensation and sublimation of gases.
[3184]	124	!
	125	! VI. Tracers
[3557]	126	! VI.1. Microphysics / Aerosols and particles.
[3539]	127	! VI.2. Updates (pressure variations, surface budget).
	128	! VI.3. Surface Tracer Update.
[3184]	129	!
	130	! VII. Surface and sub-surface soil temperature.
	131	!
	132	! VIII. Perform diagnostics and write output files.
	133	!
	134	!
	135	! arguments
	136	! ---------
	137	!
	138	! INPUT
	139	! -----
	140	!
	141	! ngrid Size of the horizontal grid.
	142	! nlayer Number of vertical layers.
	143	! nq Number of advected fields.
	144	!
	145	! firstcall True at the first call.
	146	! lastcall True at the last call.
	147	!
	148	! pday Number of days counted from the North. Spring equinoxe.
	149	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
	150	! ptimestep timestep (s).
	151	!
	152	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa).
	153	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
	154	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2.s-2).
	155	!
	156	! pu(ngrid,nlayer) u, zonal component of the wind (ms-1).
	157	! pv(ngrid,nlayer) v, meridional component of the wind (ms-1).
	158	!
	159	! pt(ngrid,nlayer) Temperature (K).
	160	!
	161	! pq(ngrid,nlayer,nq) Advected fields.
	162	!
	163	! pudyn(ngrid,nlayer) \
	164	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
	165	! ptdyn(ngrid,nlayer) / corresponding variables.
	166	! pqdyn(ngrid,nlayer,nq) /
	167	! flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
	168	!
	169	! OUTPUT
	170	! ------
	171	!
	172	! pdu(ngrid,nlayer) \
	173	! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding
	174	! pdt(ngrid,nlayer) / variables due to physical processes.
	175	! pdq(ngrid,nlayer) /
	176	! pdpsrf(ngrid) /
	177	!
	178	!
	179	! Authors
	180	! -------
	181	! Frederic Hourdin 15/10/93
	182	! Francois Forget 1994
[3195]	183	! Christophe Hourdin 02/1997
[3184]	184	! Subroutine completely rewritten by F. Forget (01/2000)
	185	! Water ice clouds: Franck Montmessin (update 06/2003)
	186	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
	187	! New correlated-k radiative scheme: R. Wordsworth (2009)
[3195]	188	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
[3184]	189	! Improved water cycle: R. Wordsworth / B. Charnay (2010)
	190	! To F90: R. Wordsworth (2010)
	191	! New turbulent diffusion scheme: J. Leconte (2012)
	192	! Loops converted to F90 matrix format: J. Leconte (2012)
[3195]	193	! No more ngrid/nq, F90 commons and adaptation to parallel: A. Spiga (2012)
[3184]	194	! Purge of the code : M. Turbet (2015)
	195	! Photochemical core developped by F. Lefevre: B. Charnay (2017)
	196	! Purge for Pluto model : A. Falco (2024)
[3539]	197	! Adapting to Pluto : A. Falco, T. Bertrand (2024)
[3557]	198	! Microphysical moment model: B. de Batz de Trenquelléon (2024)
[3184]	199	!==================================================================
	200
	201
	202	! 0. Declarations :
	203	! ------------------
	204
	205	include "netcdf.inc"
	206
	207	! Arguments :
	208	! -----------
	209
	210	! INPUTS:
	211	! -------
	212
	213	integer,intent(in) :: ngrid ! Number of atmospheric columns.
	214	integer,intent(in) :: nlayer ! Number of atmospheric layers.
	215	integer,intent(in) :: nq ! Number of tracers.
[3195]	216
[3184]	217	logical,intent(in) :: firstcall ! Signals first call to physics.
	218	logical,intent(in) :: lastcall ! Signals last call to physics.
[3195]	219
[3184]	220	real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0.
	221	real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
	222	real,intent(in) :: ptimestep ! Physics timestep (s).
	223	real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
	224	real,intent(in) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
	225	real,intent(in) :: pphi(ngrid,nlayer) ! Geopotential at mid-layer (m2s-2).
	226	real,intent(in) :: pu(ngrid,nlayer) ! Zonal wind component (m/s).
	227	real,intent(in) :: pv(ngrid,nlayer) ! Meridional wind component (m/s).
	228	real,intent(in) :: pt(ngrid,nlayer) ! Temperature (K).
	229	real,intent(in) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
	230	real,intent(in) :: flxw(ngrid,nlayer) ! Vertical mass flux (ks/s) at lower boundary of layer
	231
	232	! OUTPUTS:
	233	! --------
	234
	235	! Physical tendencies :
	236
	237	real,intent(out) :: pdu(ngrid,nlayer) ! Zonal wind tendencies (m/s/s).
	238	real,intent(out) :: pdv(ngrid,nlayer) ! Meridional wind tendencies (m/s/s).
	239	real,intent(out) :: pdt(ngrid,nlayer) ! Temperature tendencies (K/s).
	240	real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
	241	real,intent(out) :: pdpsrf(ngrid) ! Surface pressure tendency (Pa/s).
	242
	243	! Local saved variables:
	244	! ----------------------
	245	integer,save :: day_ini ! Initial date of the run (sol since Ls=0).
	246	integer,save :: icount ! Counter of calls to physiq during the run.
	247	!$OMP THREADPRIVATE(day_ini,icount)
	248
[3241]	249	!Pluto specific
[3237]	250	REAL,save :: acond,bcond
[3241]	251	REAL,save :: tcond1p4Pa
[3237]	252	DATA tcond1p4Pa/38/
	253
[3184]	254	! Local variables :
	255	! -----------------
[3241]	256	! Tendencies for the paleoclimate mode
	257	REAL qsurfyear(ngrid,nq) ! kg.m-2 averaged mass of ice lost/gained in the last Pluto year of the run
	258	REAL phisfinew(ngrid) ! geopotential of the bedrock (= phisfi-qsurf/1000*g)
	259	REAL qsurfpal(ngrid,nq) ! qsurf after a paleoclimate step : for physdem1 and restartfi
	260	REAL phisfipal(ngrid) ! geopotential after a paleoclimate step : for physdem1 and restartfi
	261	REAL oblipal ! change of obliquity
	262	REAL peri_daypal ! new periday
	263	REAL eccpal ! change of eccentricity
	264	REAL tpalnew ! change of time
	265	REAL adjustnew ! change in N2 ice albedo
	266	REAL pdaypal ! new pday = day_ini + step
	267	REAL zdt_tot ! time range corresponding to the flux of qsurfyear
	268	REAL massacc(nq) ! accumulated mass
	269	REAL masslost(nq) ! accumulated mass
[3184]	270
[3241]	271	REAL globave ! globalaverage 2D ps
	272	REAL globaveice(nq) ! globalaverage 2D ice
	273	REAL globavenewice(nq) ! globalaverage 2D ice
	274	INTEGER lecttsoil ! lecture of tsoil from proftsoil
	275	REAL qsurf1(ngrid,nq) ! saving qsurf to calculate flux over long timescales kg.m-2
	276	REAL flusurf(ngrid,nq) ! flux cond/sub kg.m-2.s-1
	277	REAL flusurfold(ngrid,nq) ! old flux cond/sub kg.m-2.s-1
[3243]	278	REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
[3241]	279
[3412]	280	REAL,SAVE :: ptime0 ! store the first time
	281	REAL dstep
	282	REAL,SAVE :: glastep=20 ! step in pluto day to spread glacier
[3241]	283
[3195]	284	! Aerosol (dust or ice) extinction optical depth at reference wavelength
[3184]	285	! for the "naerkind" optically active aerosols:
	286
[3585]	287	real,save,allocatable :: dtau_aer(:,:,:) ! Aerosols
	288	!$OMP THREADPRIVATE(dtau_aer)
[3184]	289	real zh(ngrid,nlayer) ! Potential temperature (K).
	290	real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
	291	real omega(ngrid,nlayer) ! omega velocity (Pa/s, >0 when downward)
	292
[3572]	293	integer i,l,ig,ierr,iq,nw,isoil,iesp
[3195]	294
[3184]	295	real zls ! Solar longitude (radians).
	296	real zlss ! Sub solar point longitude (radians).
	297	real zday ! Date (time since Ls=0, calculated in sols).
[3329]	298	REAL,save :: saveday ! saved date
	299	REAL,save :: savedeclin ! saved declin
[3184]	300	real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers.
	301	real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries.
	302
	303	! TENDENCIES due to various processes :
	304
	305	! For Surface Temperature : (K/s)
	306	real zdtsurf(ngrid) ! Cumulated tendencies.
	307	real zdtsurfmr(ngrid) ! Mass_redistribution routine.
	308	real zdtsurfc(ngrid) ! Condense_n2 routine.
	309	real zdtsdif(ngrid) ! Turbdiff/vdifc routines.
[3195]	310
	311	! For Atmospheric Temperatures : (K/s)
[3184]	312	real dtlscale(ngrid,nlayer) ! Largescale routine.
	313	real zdtc(ngrid,nlayer) ! Condense_n2 routine.
	314	real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	315	real zdtmr(ngrid,nlayer) ! Mass_redistribution routine.
	316	real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
	317	real zdtchim(ngrid,nlayer) ! Calchim routine.
[3195]	318
[3184]	319	! For Surface Tracers : (kg/m2/s)
	320	real dqsurf(ngrid,nq) ! Cumulated tendencies.
[3228]	321	!real zdqsurfc(ngrid) ! Condense_n2 routine.
[3195]	322	REAL zdqsc(ngrid,nq) ! Condense_n2 routine.
[3184]	323	real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines.
	324	real zdqssed(ngrid,nq) ! Callsedim routine.
	325	real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine.
[3195]	326
[3184]	327	! For Tracers : (kg/kg_of_air/s)
	328	real zdqc(ngrid,nlayer,nq) ! Condense_n2 routine.
	329	real zdqadj(ngrid,nlayer,nq) ! Convadj routine.
	330	real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines.
	331	real zdqevap(ngrid,nlayer) ! Turbdiff routine.
	332	real zdqsed(ngrid,nlayer,nq) ! Callsedim routine.
	333	real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine.
	334	REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine
	335	REAL,allocatable,save :: zdqschim(:,:) ! Calchim_asis routine
	336	!$OMP THREADPRIVATE(zdqchim,zdqschim)
	337
[3195]	338	!! PLUTO variables
	339	REAL zdqch4cloud(ngrid,nlayer,nq)
	340	REAL zdqsch4cloud(ngrid,nq)
	341	REAL zdtch4cloud(ngrid,nlayer)
	342	REAL zdqcocloud(ngrid,nlayer,nq)
	343	REAL zdqscocloud(ngrid,nq)
	344	REAL zdtcocloud(ngrid,nlayer)
	345	REAL rice_ch4(ngrid,nlayer) ! Methane ice geometric mean radius (m)
	346	REAL rice_co(ngrid,nlayer) ! CO ice geometric mean radius (m)
	347
	348	REAL zdqsch4fast(ngrid) ! used only for fast model nogcm
	349	REAL zdqch4fast(ngrid) ! used only for fast model nogcm
	350	REAL zdqscofast(ngrid) ! used only for fast model nogcm
	351	REAL zdqcofast(ngrid) ! used only for fast model nogcm
	352	REAL zdqflow(ngrid,nq)
	353
[3455]	354	REAL zdtconduc(ngrid,nlayer) ! (K/s)
	355	REAL zdumolvis(ngrid,nlayer)
	356	REAL zdvmolvis(ngrid,nlayer)
	357	real zdqmoldiff(ngrid,nlayer,nq)
	358
[3195]	359	! Haze relatated tendancies
	360	REAL zdqhaze(ngrid,nlayer,nq)
	361	REAL zdqprodhaze(ngrid,nq)
	362	REAL zdqsprodhaze(ngrid)
	363	REAL zdqhaze_col(ngrid)
	364	REAL zdqphot_prec(ngrid,nlayer)
	365	REAL zdqphot_ch4(ngrid,nlayer)
	366	REAL zdqconv_prec(ngrid,nlayer)
	367	REAL zdq_source(ngrid,nlayer,nq)
	368	! Fast Haze relatated tendancies
	369	REAL fluxbot(ngrid)
	370	REAL gradflux(ngrid)
	371	REAL fluxlym_sol_bot(ngrid) ! Solar flux Lyman alpha ph.m-2.s-1 reaching the surface
	372	REAL fluxlym_ipm_bot(ngrid) ! IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 reaching the surface
	373	REAL flym_sol(ngrid) ! Incident Solar flux Lyman alpha ph.m-2.s-1
	374	REAL flym_ipm(ngrid) ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1
[3237]	375	REAL zfluxuv ! Lyman alpha flux at 1AU
[3195]	376
[3184]	377	REAL array_zero1(ngrid)
	378	REAL array_zero2(ngrid,nlayer)
[3195]	379
[3184]	380	! For Winds : (m/s/s)
	381	real zdvadj(ngrid,nlayer), zduadj(ngrid,nlayer) ! Convadj routine.
	382	real zdumr(ngrid,nlayer), zdvmr(ngrid,nlayer) ! Mass_redistribution routine.
	383	real zdvdif(ngrid,nlayer), zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	384	real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	385	real zdhadj(ngrid,nlayer) ! Convadj routine.
[3195]	386	REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer) ! condense_n2 routine.
	387
[3184]	388	! For Pressure and Mass :
	389	real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
	390	real zdmassmr_col(ngrid) ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
	391	real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s).
	392
[3557]	393	! Local variables for MICROPHYSICS:
	394	! ---------------------------------
	395	real gzlat(ngrid,nlayer) ! Altitude-Latitude-dependent gravity (this should be stored elsewhere...).
	396	real pdqmufi(ngrid,nlayer,nq) ! Microphysical tendency (X/kg_of_air/s).
	397	real pdqmufi_prod(ngrid,nlayer,nq) ! Aerosols production tendency (kg/kg_of_air/s).
	398	real int2ext(ngrid,nlayer) ! Intensive to extensive factor (kg_air/m3: X/kg_air --> X/m3).
	399
[3184]	400	! Local variables for LOCAL CALCULATIONS:
	401	! ---------------------------------------
	402	real zflubid(ngrid)
	403	real zplanck(ngrid),zpopsk(ngrid,nlayer)
[3237]	404	REAL zdum1(ngrid,nlayer)
	405	REAL zdum2(ngrid,nlayer)
[3184]	406	real ztim1,ztim2,ztim3, z1,z2
	407	real ztime_fin
	408	real zdh(ngrid,nlayer)
	409	real gmplanet
	410	real taux(ngrid),tauy(ngrid)
	411
	412	! local variables for DIAGNOSTICS : (diagfi & stat)
	413	! -------------------------------------------------
	414	real ps(ngrid) ! Surface Pressure.
	415	real zt(ngrid,nlayer) ! Atmospheric Temperature.
	416	real zu(ngrid,nlayer),zv(ngrid,nlayer) ! Zonal and Meridional Winds.
	417	real zq(ngrid,nlayer,nq) ! Atmospheric Tracers.
	418	real zdtadj(ngrid,nlayer) ! Convadj Diagnostic.
	419	real zdtdyn(ngrid,nlayer) ! Dynamical Heating (K/s).
	420	real zdudyn(ngrid,nlayer) ! Dynamical Zonal Wind tendency (m.s-2).
	421
	422	real reff(ngrid,nlayer) ! Effective dust radius (used if doubleq=T).
	423	real vmr(ngrid,nlayer) ! volume mixing ratio
	424	real time_phys
[3195]	425
[3184]	426	real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
[3195]	427
[3184]	428	real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
	429
[3195]	430	! Pluto adding variables
	431	real vmr_ch4(ngrid) ! vmr ch4
	432	real vmr_co(ngrid) ! vmr co
	433	real rho(ngrid,nlayer) ! density
	434	real zrho_ch4(ngrid,nlayer) ! density methane kg.m-3
	435	real zrho_co(ngrid,nlayer) ! density CO kg.m-3
	436	real zrho_haze(ngrid,nlayer) ! density haze kg.m-3
	437	real zdqrho_photprec(ngrid,nlayer) !photolysis rate kg.m-3.s-1
	438	real zq1temp_ch4(ngrid) !
	439	real qsat_ch4(ngrid) !
	440	real qsat_ch4_l1(ngrid) !
	441	! CHARACTER(LEN=20) :: txt ! to temporarly store text for eddy tracers
	442	real profmmr(ngrid,nlayer) ! fixed profile of haze if haze_proffix
	443	real sensiblehf1(ngrid) ! sensible heat flux
	444	real sensiblehf2(ngrid) ! sensible heat flux
	445
[3184]	446	! included by RW for H2O Manabe scheme
	447	real rneb_man(ngrid,nlayer) ! H2O cloud fraction (moistadj).
	448	real rneb_lsc(ngrid,nlayer) ! H2O cloud fraction (large scale).
	449
	450	! to test energy conservation (RW+JL)
	451	real mass(ngrid,nlayer),massarea(ngrid,nlayer)
	452	real dEtot, dEtots, AtmToSurf_TurbFlux
	453	real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
	454	!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
[3195]	455
[3184]	456	!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
[3195]	457	real dtmoist_max,dtmoist_min
[3184]	458	real dItot, dItot_tmp, dVtot, dVtot_tmp
	459	real dWtot, dWtot_tmp, dWtots, dWtots_tmp
	460
	461	real nconsMAX, vdifcncons(ngrid), cadjncons(ngrid) ! Vdfic water conservation test. By RW
	462
	463	real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW
	464
	465	! Non-oro GW tendencies
	466	REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
	467	REAL d_t_hin(ngrid,nlayer)
	468	! Diagnostics 2D of gw_nonoro
	469	REAL zustrhi(ngrid), zvstrhi(ngrid)
	470
	471	real :: tsurf2(ngrid)
	472	!! real :: flux_o(ngrid),flux_g(ngrid)
	473	real :: flux_g(ngrid)
	474	real :: flux_sens_lat(ngrid)
	475	real :: qsurfint(ngrid,nq)
	476
[3483]	477	! local variables for skin depth check
[3539]	478	real :: therm_inertia(ngrid,nsoilmx)
[3483]	479	real :: inertia_min,inertia_max
	480	real :: diurnal_skin ! diurnal skin depth (m)
	481	real :: annual_skin ! anuual skin depth (m)
	482
[3500]	483	! when no startfi file is asked for init
	484	real alpha,lay1 ! coefficients for building layers
	485	integer iloop
[3483]	486
[3184]	487	! flags to trigger extra sanity checks
	488	logical, save :: check_physics_inputs=.false.
	489	logical, save :: check_physics_outputs=.false.
[3195]	490	!$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)
[3184]	491
	492	! Misc
	493	character*2 :: str2
[3195]	494	character(len=10) :: tmp1
[3184]	495	character(len=10) :: tmp2
	496	!==================================================================================================
	497
	498	! -----------------
	499	! I. INITIALISATION
	500	! -----------------
	501
	502	! --------------------------------
	503	! I.1 First Call Initialisation.
	504	! --------------------------------
	505	if (firstcall) then
	506	call getin_p("check_physics_inputs", check_physics_inputs)
	507	call getin_p("check_physics_outputs", check_physics_outputs)
	508
	509	! Allocate saved arrays (except for 1D model, where this has already
	510	! been done)
	511	if (ngrid>1) call phys_state_var_init(nq)
	512
	513	! Variables set to 0
	514	! ~~~~~~~~~~~~~~~~~~
	515	dtrad(:,:) = 0.0
	516	fluxrad(:) = 0.0
	517	tau_col(:) = 0.0
	518	zdtsw(:,:) = 0.0
	519	zdtlw(:,:) = 0.0
	520
	521	! Initialize tracer names, indexes and properties.
	522	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	523	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
	524	if (tracer) then
	525	call initracer(ngrid,nq)
	526	! if(photochem) then !AF24: removed
	527	endif
	528	! Initialize aerosol indexes.
	529	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3195]	530	! call iniaerosol
[3184]	531	! allocate related local arrays
	532	! (need be allocated instead of automatic because of "naerkind")
[3585]	533	allocate(dtau_aer(ngrid,nlayer,naerkind))
[3184]	534
	535	#ifdef CPP_XIOS
	536	! Initialize XIOS context
	537	write(,) "physiq: call wxios_context_init"
	538	CALL wxios_context_init
	539	#endif
	540
	541	! Read 'startfi.nc' file.
	542	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	543	call phyetat0(startphy_file, &
	544	ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
[3483]	545	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,inertiedat)
[3184]	546
	547	if (.not.startphy_file) then
[3500]	548	! starting without startfi.nc and with callsoil
	549	! is not yet possible as soildepth default is not defined
	550	if (callsoil) then
	551	! default mlayer distribution, following a power law:
	552	! mlayer(k)=lay1alpha*(k-1/2)
	553	lay1=2.e-4
	554	alpha=2
	555	do iloop=0,nsoilmx-1
[3539]	556	mlayer(iloop)=lay1(alpha*(iloop-0.5))
	557	enddo
[3500]	558	lay1=sqrt(mlayer(0)*mlayer(1))
	559	alpha=mlayer(1)/mlayer(0)
	560	do iloop=1,nsoilmx
	561	layer(iloop)=lay1(alpha*(iloop-1))
	562	enddo
	563	endif
[3184]	564	! additionnal "academic" initialization of physics
	565	if (is_master) write(,) "Physiq: initializing tsurf(:) to pt(:,1) !!"
	566	tsurf(:)=pt(:,1)
	567	if (is_master) write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
	568	do isoil=1,nsoilmx
	569	tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
	570	enddo
[3232]	571	if (is_master) write(,) "Physiq: initializing day_ini to pday !"
[3184]	572	day_ini=pday
	573	endif
[3539]	574
[3184]	575	if (pday.ne.day_ini) then
	576	write(,) "ERROR in physiq.F90:"
	577	write(,) "bad synchronization between physics and dynamics"
	578	write(,) "dynamics day: ",pday
	579	write(,) "physics day: ",day_ini
	580	stop
	581	endif
	582
	583	write (,) 'In physiq day_ini =', day_ini
	584
[3195]	585	! Initialize albedo calculation.
[3184]	586	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	587	albedo(:,:)=0.0
	588	albedo_bareground(:)=0.0
	589	albedo_snow_SPECTV(:)=0.0
	590	albedo_n2_ice_SPECTV(:)=0.0
[3412]	591
	592	ptime0=ptime
	593	write (,) 'In physiq ptime0 =', ptime
	594
[3195]	595	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV)
[3184]	596
[3195]	597	! Initialize orbital calculation.
[3184]	598	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	599	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
	600
[3329]	601	savedeclin=0.
	602	saveday=pday
[3377]	603	adjust=0. ! albedo adjustment for convergeps
[3184]	604
	605	! Initialize soil.
	606	! ~~~~~~~~~~~~~~~~
	607	if (callsoil) then
	608	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
	609	ptimestep,tsurf,tsoil,capcal,fluxgrd)
	610	else ! else of 'callsoil'.
	611	print*,'WARNING! Thermal conduction in the soil turned off'
	612	capcal(:)=1.e6
	613	fluxgrd(:)=intheat
	614	print*,'Flux from ground = ',intheat,' W m^-2'
	615	endif ! end of 'callsoil'.
	616
	617	icount=1
	618
	619	! Initialize variable for dynamical heating and zonal wind tendency diagnostic
	620	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	621	ztprevious(:,:)=pt(:,:)
	622	zuprevious(:,:)=pu(:,:)
	623
[3195]	624	if(meanOLR)then
	625	call system('rm -f rad_bal.out') ! to record global radiative balance.
	626	call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
[3184]	627	call system('rm -f h2o_bal.out') ! to record global hydrological balance.
	628	endif
	629
	630	if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
	631	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
	632	ptimestep,pday+nday,time_phys,cell_area, &
[3539]	633	albedo_bareground,zmea,zstd,zsig,zgam,zthe)
[3184]	634	endif
	635
[3557]	636	! Initialize correlated-k.
	637	! ~~~~~~~~~~~~~~~~~~~~~~~~
[3195]	638	if (corrk) then
	639	! We initialise the spectral grid here instead of
	640	! at firstcall of callcorrk so we can output XspecIR, XspecVI
[3184]	641	! when using Dynamico
[3504]	642	if (is_master) print*, "physiq_mod: Correlated-k data base folder:",trim(datadir)
[3184]	643	call getin_p("corrkdir",corrkdir)
[3504]	644	if (is_master) print*,"corrkdir = ", corrkdir
[3184]	645	write (tmp1, '(i4)') L_NSPECTI
	646	write (tmp2, '(i4)') L_NSPECTV
	647	banddir=trim(trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2)))
	648	banddir=trim(trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir)))
	649	call setspi !Basic infrared properties.
	650	call setspv ! Basic visible properties.
	651	call sugas_corrk ! Set up gaseous absorption properties.
[3572]	652	if (optichaze) then
[3195]	653	call suaer_corrk ! Set up aerosol optical properties.
	654	endif
[3184]	655	endif
	656
[3557]	657	! Initialize microphysics.
	658	! ~~~~~~~~~~~~~~~~~~~~~~~~
	659	IF (callmufi) THEN
	660	! Initialize microphysics arrays.
	661	call inimufi(ptimestep)
	662	ENDIF ! end callmufi
	663
[3275]	664	!! call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
[3184]	665	! XIOS outputs
	666	#ifdef CPP_XIOS
	667
[3506]	668	if (is_master) write(,) "physiq: call initialize_xios_output"
	669	call initialize_xios_output(pday,ptime,ptimestep,daysec,year_day, &
	670	presnivs,pseudoalt,mlayer,WNOI,WNOV)
[3184]	671	#endif
[3195]	672
[3275]	673	!! call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
[3184]	674
	675	write(,) "physiq: end of firstcall"
	676	endif ! end of 'firstcall'
	677
[3275]	678	!! call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
	679	!! call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
[3184]	680
	681	if (check_physics_inputs) then
	682	!check the validity of input fields coming from the dynamics
	683	call check_physics_fields("begin physiq:", pt, pu, pv, pplev, pq)
	684	endif
	685
	686	! call writediagfi(ngrid,"pre_physical_rnat"," "," ",2,rnat)
	687	! call writediagfi(ngrid,"pre_physical_capcal"," "," ",2,capcal)
	688
	689	! ------------------------------------------------------
	690	! I.2 Initializations done at every physical timestep:
	691	! ------------------------------------------------------
	692
[3195]	693	#ifdef CPP_XIOS
[3184]	694	! update XIOS time/calendar
	695	call update_xios_timestep
[3195]	696	#endif
[3184]	697
	698	! Initialize various variables
	699	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3195]	700
[3539]	701	pdt(1:ngrid,1:nlayer) = 0.0
[3184]	702	zdtsurf(1:ngrid) = 0.0
	703	pdq(1:ngrid,1:nlayer,1:nq) = 0.0
	704	dqsurf(1:ngrid,1:nq)= 0.0
	705	pdu(1:ngrid,1:nlayer) = 0.0
	706	pdv(1:ngrid,1:nlayer) = 0.0
	707	pdpsrf(1:ngrid) = 0.0
[3195]	708	zflubid(1:ngrid) = 0.0
[3184]	709	flux_sens_lat(1:ngrid) = 0.0
	710	taux(1:ngrid) = 0.0
	711	tauy(1:ngrid) = 0.0
	712
[3412]	713	if (conservn2) then
	714	write(,) 'conservn2 iniloop'
	715	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
	716	endif
	717
[3184]	718	zday=pday+ptime ! Compute time, in sols (and fraction thereof).
	719
	720	! Compute Stellar Longitude (Ls), and orbital parameters.
	721	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	722	if (season) then
	723	call stellarlong(zday,zls)
	724	else
	725	call stellarlong(noseason_day,zls)
	726	end if
	727
[3237]	728	! Get Lyman alpha flux at specific Ls
[3572]	729	if (callmufi.or.haze) then
[3237]	730	call lymalpha(zls,zfluxuv)
	731	print*, 'Haze lyman-alpha zls,zfluxuv=',zls,zfluxuv
	732	end if
	733
	734	IF (triton) then
	735	CALL orbitetriton(zls,zday,dist_star,declin)
	736	ELSE
	737	call orbite(zls,dist_star,declin,right_ascen)
	738	ENDIF
	739
[3184]	740	if (diurnal) then
	741	zlss=-2.pi(zday-.5)
	742	else if(diurnal .eqv. .false.) then
	743	zlss=9999.
[3195]	744	endif
[3184]	745
	746	glat(:) = g !AF24: removed oblateness
	747
	748	! Compute geopotential between layers.
	749	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	750	zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
[3195]	751	do l=1,nlayer
[3184]	752	zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
	753	enddo
	754
	755	zzlev(1:ngrid,1)=0.
	756
	757	do l=2,nlayer
	758	do ig=1,ngrid
	759	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
	760	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
	761	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	762	enddo
[3195]	763	enddo
[3184]	764
	765	!Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
[3195]	766
[3184]	767	zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
	768
	769	! Compute potential temperature
	770	! Note : Potential temperature calculation may not be the same in physiq and dynamic...
	771	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3195]	772	do l=1,nlayer
[3184]	773	do ig=1,ngrid
	774	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
	775	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
	776	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
	777	massarea(ig,l)=mass(ig,l)*cell_area(ig)
	778	enddo
	779	enddo
	780
	781	! Compute vertical velocity (m/s) from vertical mass flux
	782	! w = F / (rhoarea) and rho = P/(rT)
	783	! But first linearly interpolate mass flux to mid-layers
[3237]	784	if (.not.fast) then
[3228]	785	do l=1,nlayer-1
[3184]	786	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
[3228]	787	enddo
	788	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
	789	do l=1,nlayer
[3184]	790	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) / &
	791	(pplay(1:ngrid,l)*cell_area(1:ngrid))
[3228]	792	enddo
	793	! omega in Pa/s
	794	do l=1,nlayer-1
[3184]	795	omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
[3228]	796	enddo
	797	omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
	798	do l=1,nlayer
[3184]	799	omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid)
[3228]	800	enddo
	801	endif
[3412]	802
	803	if (conservn2) then
	804	write(,) 'conservn2 thermo'
	805	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
	806	endif
	807
[3557]	808	! Compute variations of g with latitude (to do).
	809	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	810	gzlat(:,:) = g
	811
	812	! Initialize microphysical diagnostics.
	813	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	814	IF (callmufi) THEN
	815	! Initialize intensive to extensive factor (kg_air/m3: X/kg_air --> X/m3).
	816	int2ext(:,:) = (pplev(:,1:nlayer)-pplev(:,2:nlayer+1)) / gzlat(:,1:nlayer) / (zzlev(:,2:nlayer+1)-zzlev(:,1:nlayer))
	817
	818	! Initialize microphysics diagnostics arrays.
	819	call inimufi_diag(ngrid,nlayer,nq,pq,int2ext)
	820	ENDIF ! end callmufi
	821
[3455]	822	! --------------------------------------------------------
[3539]	823	! II.1 Thermosphere
[3455]	824	! --------------------------------------------------------
	825
	826	! ajout de la conduction depuis la thermosphere
	827	IF (callconduct) THEN
	828
	829	call conduction (ngrid,nlayer,ptimestep, &
	830	pplay,pt,zzlay,zzlev,zdtconduc,tsurf)
	831	DO l=1,nlayer
	832	DO ig=1,ngrid
	833	pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
	834	ENDDO
	835	ENDDO
	836
	837	ENDIF
	838
	839	! ajout de la viscosite moleculaire
	840	IF (callmolvis) THEN
	841	call molvis(ngrid,nlayer,ptimestep, &
	842	pplay,pt,zzlay,zzlev, &
	843	zdtconduc,pu,tsurf,zdumolvis)
	844	call molvis(ngrid,nlayer,ptimestep, &
	845	pplay,pt,zzlay,zzlev, &
	846	zdtconduc,pv,tsurf,zdvmolvis)
	847
	848	DO l=1,nlayer
	849	DO ig=1,ngrid
	850	! pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
	851	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
	852	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
	853	ENDDO
	854	ENDDO
	855	ENDIF
	856
	857	IF (callmoldiff) THEN
	858	call moldiff_red(ngrid,nlayer,nq, &
	859	pplay,pplev,pt,pdt,pq,pdq,ptimestep, &
	860	zzlay,zdtconduc,zdqmoldiff)
	861
	862	DO l=1,nlayer
	863	DO ig=1,ngrid
	864	DO iq=1, nq
	865	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
	866	ENDDO
	867	ENDDO
	868	ENDDO
	869	ENDIF
	870
	871	if (conservn2) then
[3539]	872	write(,) 'conservn2 thermosphere'
[3455]	873	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
	874	endif
	875
	876
[3184]	877	!---------------------------------
[3539]	878	! II.2 Compute radiative tendencies
[3184]	879	!---------------------------------
[3243]	880	! Saving qsurf to compute paleo flux condensation/sublimation
	881	DO iq=1, nq
	882	DO ig=1,ngrid
	883	IF (qsurf(ig,iq).lt.0.) then
	884	qsurf(ig,iq)=0.
	885	ENDIF
	886	qsurf1(ig,iq)=qsurf(ig,iq)
	887	ENDDO
	888	ENDDO
[3184]	889
[3243]	890
[3184]	891	! Compute local stellar zenith angles
	892	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3377]	893	fract = 0
[3184]	894	if (diurnal) then
	895	ztim1=SIN(declin)
	896	ztim2=COS(declin)COS(2.pi*(zday-.5))
	897	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
	898
	899	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
	900	ztim1,ztim2,ztim3,mu0,fract)
	901	else if(diurnal .eqv. .false.) then
	902
	903	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad)
	904	! WARNING: this function appears not to work in 1D
	905
	906	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
	907	mu0 = cos(pi*szangle/180.0)
[3329]	908	fract= 1/(4*mu0) ! AF24: from pluto.old
[3184]	909	endif
	910
[3195]	911	endif
[3184]	912
[3237]	913
[3329]	914	! Pluto albedo /IT changes depending on surface ices (only in 3D)
	915	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	916	if (ngrid.ne.1) then
	917
	918	!! Specific to change albedo of N2 so that Psurf
	919	!! converges toward 1.4 Pa in "1989" seasons for Triton
	920	!! converges toward 1.1 Pa in "2015" seasons for Pluto
	921	if (convergeps) then
	922	if (triton) then
	923	! 1989 declination
	924	if (declin180./pi.gt.-46..and.declin180./pi.lt.-45. &
	925	.and.zday.gt.saveday+1000. &
	926	.and.declin.lt.savedeclin) then
[3539]	927	call planetwide_sumval(pplev(:,1)*cell_area(:)/totarea_planet,globave)
[3329]	928	if (globave.gt.1.5) then
	929	adjust=adjust+0.005
	930	else if (globave.lt.1.3) then
	931	adjust=adjust-0.005
	932	endif
	933	saveday=zday
	934	endif
	935	else
	936	! Pluto : 2015 declination current epoch
	937	if (declin180./pi.gt.50.and.declin180./pi.lt.51. &
	938	.and.zday.gt.saveday+10000. &
	939	.and.declin.gt.savedeclin) then
[3539]	940	call planetwide_sumval(pplev(:,1)*cell_area(:)/totarea_planet,globave)
[3329]	941	if (globave.gt.1.2) then
	942	adjust=adjust+0.005
	943	else if (globave.lt.1.) then
	944	adjust=adjust-0.005
	945	endif
	946	saveday=zday
	947	endif
	948	endif
	949	end if
	950	end if ! if ngrid ne 1
	951
[3483]	952	call surfprop(ngrid,nq,fract,qsurf,tsurf, &
	953	capcal,adjust,dist_star,flusurfold,ptimestep,zls,&
[3539]	954	albedo,emis,therm_inertia)
[3329]	955	! do i=2,ngrid
	956	! albedo(i,:) = albedo(1,:)
	957	! enddo
[3544]	958	! AF24: TODO check albedo has been initialized here
[3329]	959
[3500]	960	if (firstcall.and.callsoil) then
[3539]	961	! AF24 Originally in soil.F, but therm_inertia is modified by surfprop
[3483]	962	! Additional checks: is the vertical discretization sufficient
	963	! to resolve diurnal and annual waves?
	964	do ig=1,ngrid
	965	! extreme inertia for this column
[3539]	966	inertia_min=minval(therm_inertia(ig,:))
	967	inertia_max=maxval(therm_inertia(ig,:))
[3483]	968	! diurnal and annual skin depth
	969	diurnal_skin=(inertia_min/volcapa)*sqrt(daysec/pi)
	970	annual_skin=(inertia_max/volcapa)sqrt(year_daydaysec/pi)
	971	if (0.5*diurnal_skin<layer(1)) then
	972	! one should have the fist layer be at least half of diurnal skin depth
	973	write(,) "soil Error: grid point ig=",ig
	974	write(,) " longitude=",longitude(ig)*(180./pi)
	975	write(,) " latitude=",latitude(ig)*(180./pi)
	976	write(,) " first soil layer depth ",layer(1)
	977	write(,) " not small enough for a diurnal skin depth of ", &
	978	diurnal_skin
	979	write(,) " change soil layer distribution (comsoil_h.F90)"
	980	call abort_physic("soil","change soil layer distribution (comsoil_h.F90)",1)
	981	endif
	982	if (2.*annual_skin>layer(nsoilmx)) then
	983	! one should have the full soil be at least twice the diurnal skin depth
	984	write(,) "soil Error: grid point ig=",ig
	985	write(,) " longitude=",longitude(ig)*(180./pi)
	986	write(,) " latitude=",latitude(ig)*(180./pi)
	987	write(,) " total soil layer depth ",layer(nsoilmx)
	988	write(,) " not large enough for an annual skin depth of ", &
	989	annual_skin
	990	write(,) " change soil layer distribution (comsoil_h.F90)"
	991	call abort_physic("soil","change soil layer distribution (comsoil_h.F90)",1)
	992	endif
	993	enddo ! of do ig=1,ngrid
	994
	995	end if ! callsoil
	996
[3184]	997	if (callrad) then
	998	if( mod(icount-1,iradia).eq.0.or.lastcall) then
[3195]	999
[3184]	1000	! Eclipse incoming sunlight !AF24: removed
	1001
[3275]	1002	!! call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
	1003	!! call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
[3184]	1004
	1005
	1006	if (corrk) then
[3195]	1007
[3184]	1008	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1009	! II.a Call correlated-k radiative transfer scheme
	1010	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1011	if(kastprof)then
	1012	print*,'kastprof should not = true here'
	1013	call abort
	1014	endif
[3195]	1015
[3184]	1016	! standard callcorrk
[3329]	1017	if (oldplutocorrk) then
	1018	call callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf, &
[3585]	1019	albedo(:,1),emis,mu0,pplev,pplay,pt, &
	1020	zzlay,zzlev,tsurf,fract,dist_star,dtau_aer, &
[3329]	1021	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, &
	1022	fluxabs_sw,fluxtop_dn,reffrad,tau_col,ptime,pday, &
	1023	firstcall,lastcall)
[3353]	1024	albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
[3361]	1025	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+ &
	1026	fluxsurf_sw(1:ngrid)*(1.-albedo(1:ngrid,1))
[3329]	1027	else
[3572]	1028	muvar(1:ngrid,1:nlayer+1)=mugaz
[3329]	1029	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
[3184]	1030	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
[3585]	1031	zzlay,zzlev,tsurf,fract,dist_star,dtau_aer,muvar, &
[3184]	1032	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, &
	1033	fluxsurfabs_sw,fluxtop_lw, &
	1034	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,GSR_nu, &
	1035	int_dtaui,int_dtauv, &
[3572]	1036	tau_col,firstcall,lastcall)
[3361]	1037	! Radiative flux from the sky absorbed by the surface (W.m-2).
	1038	GSR=0.0
	1039	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+ &
	1040	fluxsurfabs_sw(1:ngrid)
[3329]	1041	endif ! oldplutocorrk
[3184]	1042	!GG (feb2021): Option to "artificially" decrease the raditive time scale in
	1043	!the deep atmosphere press > 0.1 bar. Suggested by MT.
[3195]	1044	!! COEFF_RAD to be "tuned" to facilitate convergence of tendency
	1045
[3184]	1046	!coeff_rad=0. ! 0 values, it doesn't accelerate the convergence
	1047	!coeff_rad=0.5
[3195]	1048	!coeff_rad=1.
[3184]	1049	!do l=1, nlayer
	1050	! do ig=1,ngrid
	1051	! if(pplay(ig,l).ge.1.d4) then
	1052	! zdtsw(ig,l)=zdtsw(ig,l)(pplay(ig,l)/1.d4)*coeff_rad
	1053	! zdtlw(ig,l)=zdtlw(ig,l)(pplay(ig,l)/1.d4)*coeff_rad
	1054	! endif
	1055	! enddo
	1056	!enddo
	1057
	1058	! AF24: removed CLFvarying Option
	1059
	1060
	1061	!if(noradsurf)then ! no lower surface; SW flux just disappears
	1062	! GSR = SUM(fluxsurf_sw(1:ngrid)*cell_area(1:ngrid))/totarea
	1063	! fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
	1064	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
	1065	!endif
	1066
	1067	! Net atmospheric radiative heating rate (K.s-1)
	1068	dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
[3195]	1069
	1070	! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that !
[3184]	1071	if (firstcall .and. albedo_spectral_mode) then
	1072	call spectral_albedo_calc(albedo_snow_SPECTV)
	1073	endif
	1074
	1075	else
	1076	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3539]	1077	! II.b Atmosphere has no radiative effect
[3184]	1078	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1079	fluxtop_dn(1:ngrid) = fract(1:ngrid)mu0(1:ngrid)Fat1AU/dist_star**2
	1080	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
	1081	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
	1082	endif
	1083	fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
	1084	print*,'------------WARNING---WARNING------------' ! by MT2015.
	1085	print*,'You are in corrk=false mode, '
[3195]	1086	print*,'and the surface albedo is taken equal to the first visible spectral value'
	1087
[3198]	1088	albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
[3184]	1089	fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
[3198]	1090	fluxabs_sw(1:ngrid)=fluxsurfabs_sw(1:ngrid)
[3184]	1091	fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid)
	1092	fluxtop_lw(1:ngrid) = emis(1:ngrid)sigmatsurf(1:ngrid)**4
	1093
	1094	dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating
	1095
	1096	endif ! end of corrk
	1097
	1098	endif ! of if(mod(icount-1,iradia).eq.0)
	1099
[3195]	1100
[3184]	1101	! Transformation of the radiative tendencies
	1102	! ------------------------------------------
	1103	zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
	1104	zplanck(1:ngrid)=emis(1:ngrid)sigmazplanck(1:ngrid)*zplanck(1:ngrid)
	1105	fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
	1106	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
[3195]	1107
[3184]	1108	! Test of energy conservation
	1109	!----------------------------
	1110	if(enertest)then
	1111	call planetwide_sumval(cppmassarea(:,:)zdtsw(:,:)/totarea_planet,dEtotSW)
	1112	call planetwide_sumval(cppmassarea(:,:)zdtlw(:,:)/totarea_planet,dEtotLW)
	1113	!call planetwide_sumval(fluxsurf_sw(:)(1.-albedo_equivalent(:))cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
	1114	call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
	1115	call planetwide_sumval((fluxsurf_lw(:)emis(:)-zplanck(:))cell_area(:)/totarea_planet,dEtotsLW)
	1116	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
	1117	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
	1118	if (is_master) then
	1119	print*,'---------------------------------------------------------------'
	1120	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
	1121	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
	1122	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
	1123	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
	1124	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
	1125	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
	1126	endif
	1127	endif ! end of 'enertest'
	1128
	1129	endif ! of if (callrad)
	1130
[3275]	1131	!! call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
	1132	!! call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
[3184]	1133
[3412]	1134	if (conservn2) then
	1135	write(,) 'conservn2 radiat'
	1136	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
	1137	endif
[3184]	1138
	1139	! --------------------------------------------
	1140	! III. Vertical diffusion (turbulent mixing) :
	1141	! --------------------------------------------
	1142
	1143	if (calldifv) then
[3195]	1144
[3184]	1145	zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
	1146
[3258]	1147	if (oldplutovdifc) then
[3244]	1148	zdum1(:,:) = 0
	1149	zdum2(:,:) = 0
	1150	zdh(:,:)=pdt(:,:)/zpopsk(:,:)
[3195]	1151
[3237]	1152	! Calling vdif (Martian version WITH N2 condensation)
	1153	CALL vdifc_pluto(ngrid,nlayer,nq,zpopsk, &
	1154	ptimestep,capcal,lwrite, &
	1155	pplay,pplev,zzlay,zzlev,z0, &
	1156	pu,pv,zh,pq,pt,tsurf,emis,qsurf, &
	1157	zdum1,zdum2,zdh,pdq,pdt,zflubid, &
	1158	zdudif,zdvdif,zdhdif,zdtsdif,q2, &
	1159	zdqdif,zdqsdif,qsat_ch4,qsat_ch4_l1) !,zq1temp_ch4,qsat_ch4)
	1160
	1161	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
	1162
	1163	bcond=1./tcond1p4Pa
	1164	acond=r/lw_n2
	1165
[3244]	1166	! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
	1167	else if (UseTurbDiff) then
[3237]	1168
[3244]	1169	call turbdiff(ngrid,nlayer,nq, &
	1170	ptimestep,capcal, &
	1171	pplay,pplev,zzlay,zzlev,z0, &
	1172	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
	1173	pdt,pdq,zflubid, &
	1174	zdudif,zdvdif,zdtdif,zdtsdif, &
	1175	sensibFlux,q2,zdqdif,zdqevap,zdqsdif, &
	1176	taux,tauy)
	1177
[3572]	1178	else ! if .not. (oldplutovdifc) .and. (UseTurbDiff)
[3244]	1179
[3184]	1180	zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
[3195]	1181
[3184]	1182	call vdifc(ngrid,nlayer,nq,zpopsk, &
	1183	ptimestep,capcal,lwrite, &
	1184	pplay,pplev,zzlay,zzlev,z0, &
	1185	pu,pv,zh,pq,tsurf,emis,qsurf, &
	1186	zdh,pdq,zflubid, &
	1187	zdudif,zdvdif,zdhdif,zdtsdif, &
	1188	sensibFlux,q2,zdqdif,zdqsdif)
	1189
	1190	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
	1191	zdqevap(1:ngrid,1:nlayer)=0.
	1192
	1193	end if !end of 'UseTurbDiff'
	1194
	1195	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
	1196
	1197	!!! this is always done, except for turbulence-resolving simulations
	1198	if (.not. turb_resolved) then
	1199	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
	1200	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
	1201	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
	1202	endif
	1203
	1204	! if(ok_slab_ocean)then !AF24: removed
	1205	! flux_sens_lat(1:ngrid)=(zdtsdif(1:ngrid)*capcal(1:ngrid)-fluxrad(1:ngrid))
	1206	! endif
	1207
[3275]	1208	!! call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_gas))
[3184]	1209
[3195]	1210	if (tracer) then
[3184]	1211	pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
	1212	dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
	1213	end if ! of if (tracer)
	1214
[3275]	1215	!! call writediagfi(ngrid,"vdifc_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
	1216	!! call writediagfi(ngrid,"vdifc_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
[3184]	1217
	1218	! test energy conservation
	1219	!-------------------------
	1220	if(enertest)then
[3195]	1221
[3184]	1222	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
	1223	do ig = 1, ngrid
	1224	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
	1225	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
	1226	enddo
[3195]	1227
[3184]	1228	call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
	1229	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
	1230	call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
	1231	call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
[3195]	1232
[3184]	1233	if (is_master) then
[3195]	1234
[3184]	1235	if (UseTurbDiff) then
	1236	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
	1237	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
	1238	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
	1239	else
	1240	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
	1241	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
	1242	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
	1243	end if
	1244	endif ! end of 'is_master'
[3195]	1245
[3184]	1246	! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
	1247	endif ! end of 'enertest'
	1248
[3539]	1249	else ! calldifv
[3184]	1250
[3539]	1251	ztim1=4.sigmaptimestep
	1252	DO ig=1,ngrid
	1253	ztim2=ztim1emis(ig)tsurf(ig)**3
	1254	z1=capcal(ig)*tsurf(ig)+ &
	1255	ztim2tsurf(ig)+ (fluxrad(ig)+fluxgrd(ig))ptimestep
	1256	z2= capcal(ig)+ztim2
	1257	zdtsurf(ig)=(z1/z2 - tsurf(ig))/ptimestep
[3184]	1258
[3539]	1259	! for output:
	1260	!dplanck(ig)=4.stephanptimestepemis(ig)tsurf(ig)**3
	1261	ENDDO
[3184]	1262
	1263	! if(.not.newtonian)then
[3539]	1264	!zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
[3184]	1265
[3421]	1266	! ------------------------------------------------------------------
	1267	! Methane surface sublimation and condensation in fast model (nogcm)
	1268	! ------------------------------------------------------------------
	1269	if ((methane).and.(fast).and.condmetsurf) THEN
	1270
[3539]	1271	call ch4surf(ngrid,nlayer,nq,ptimestep,capcal, &
[3421]	1272	tsurf,zdtsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf, &
	1273	zdqch4fast,zdqsch4fast)
	1274
	1275	dqsurf(1:ngrid,igcm_ch4_ice)= dqsurf(1:ngrid,igcm_ch4_ice) + &
	1276	zdqsch4fast(1:ngrid)
	1277	pdq(1:ngrid,1,igcm_ch4_gas)= pdq(1:ngrid,1,igcm_ch4_gas) + &
	1278	zdqch4fast(1:ngrid)
	1279	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+lw_ch4*zdqsch4fast(1:ngrid)/capcal(1:ngrid)
	1280	end if
	1281	! ------------------------------------------------------------------
	1282	! CO surface sublimation and condensation in fast model (nogcm)
	1283	! ------------------------------------------------------------------
	1284	if ((carbox).and.(fast).and.condcosurf) THEN
	1285
	1286	call cosurf(ngrid,nlayer,nq,ptimestep, &
	1287	tsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf, &
	1288	zdqcofast,zdqscofast)
	1289
	1290	dqsurf(1:ngrid,igcm_co_ice)= dqsurf(1:ngrid,igcm_co_ice) + &
	1291	zdqscofast(1:ngrid)
	1292	pdq(1:ngrid,1,igcm_co_gas)= pdq(1:ngrid,1,igcm_co_gas) + &
	1293	zdqcofast(1:ngrid)
	1294	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+lw_co*zdqscofast(1:ngrid)/capcal(1:ngrid)
	1295	end if
	1296
	1297
[3184]	1298	endif ! end of 'calldifv'
	1299
[3412]	1300	if (conservn2) then
[3539]	1301	write(,) 'conservn2 calldifv'
[3412]	1302	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)+ &
	1303	dqsurf(:,1)*ptimestep)
	1304	endif
[3539]	1305	if (methane.and.conservch4) then
	1306	write(,) 'conservch4 calldifv'
	1307	if (fast) then
	1308	call testconservfast(ngrid,nlayer,nq,pq(:,1,igcm_ch4_gas),pdq(:,1,igcm_ch4_gas), &
	1309	qsurf(:,igcm_ch4_ice),dqsurf(:,igcm_ch4_ice), &
	1310	ptimestep,pplev,zdqch4fast,zdqsch4fast,'CH4',' vdifc ')
	1311	else
	1312	call testconserv(ngrid,nlayer,nq,pq,pdq,qsurf,dqsurf, &
	1313	igcm_ch4_gas,igcm_ch4_ice, &
	1314	ptimestep,pplev,zdqdif,zdqsdif,'CH4',' vdifc ')
	1315	endif
	1316	endif
[3184]	1317
	1318	!-------------------
	1319	! IV. Convection :
	1320	!-------------------
[3195]	1321
[3184]	1322	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3539]	1323	! IV.a Dry convective adjustment :
[3184]	1324	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1325
	1326	if(calladj) then
	1327
	1328	zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
	1329	zduadj(1:ngrid,1:nlayer)=0.0
	1330	zdvadj(1:ngrid,1:nlayer)=0.0
	1331	zdhadj(1:ngrid,1:nlayer)=0.0
	1332
	1333
	1334	call convadj(ngrid,nlayer,nq,ptimestep, &
	1335	pplay,pplev,zpopsk, &
	1336	pu,pv,zh,pq, &
	1337	pdu,pdv,zdh,pdq, &
	1338	zduadj,zdvadj,zdhadj, &
	1339	zdqadj)
	1340
	1341	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
	1342	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer)
	1343	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
	1344	zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
	1345
[3195]	1346	if(tracer) then
	1347	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
[3184]	1348	end if
	1349
	1350	! Test energy conservation
	1351	if(enertest)then
	1352	call planetwide_sumval(cppmassarea(:,:)zdtadj(:,:)/totarea_planet,dEtot)
	1353	if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
	1354	endif
	1355
	1356	! ! Test water conservation !AF24: removed
[3195]	1357
[3184]	1358	endif ! end of 'calladj'
[3195]	1359
[3184]	1360	!-----------------------------------------------
[3195]	1361	! V. Nitrogen condensation-sublimation :
[3184]	1362	!-----------------------------------------------
	1363
	1364	if (n2cond) then
[3195]	1365
[3184]	1366	if (.not.tracer) then
	1367	print*,'We need a N2 ice tracer to condense N2'
	1368	call abort
	1369	endif
[3195]	1370
[3390]	1371	zdqc(:,:,:)=0.
	1372	zdqsc(:,:)=0.
[3184]	1373	call condense_n2(ngrid,nlayer,nq,ptimestep, &
[3195]	1374	capcal,pplay,pplev,tsurf,pt, &
	1375	pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, &
	1376	qsurf(1,igcm_n2),albedo,emis, &
	1377	zdtc,zdtsurfc,pdpsrf,zduc,zdvc, &
	1378	zdqc,zdqsc(1,igcm_n2))
[3184]	1379
	1380	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
[3228]	1381	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer)+zdvc(1:ngrid,1:nlayer)
	1382	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer)+zduc(1:ngrid,1:nlayer)
[3184]	1383	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
	1384
	1385	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq)
[3228]	1386	dqsurf(1:ngrid,igcm_n2) = dqsurf(1:ngrid,igcm_n2) + zdqsc(1:ngrid,igcm_n2)
[3184]	1387
[3275]	1388	!! call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
	1389	!! call writediagfi(ngrid,"condense_n2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
[3184]	1390
	1391	! test energy conservation
	1392	if(enertest)then
	1393	call planetwide_sumval(cppmassarea(:,:)zdtc(:,:)/totarea_planet,dEtot)
	1394	call planetwide_sumval(capcal(:)zdtsurfc(:)cell_area(:)/totarea_planet,dEtots)
	1395	if (is_master) then
	1396	print*,'In n2cloud atmospheric energy change =',dEtot,' W m-2'
	1397	print*,'In n2cloud surface energy change =',dEtots,' W m-2'
	1398	endif
	1399	endif
	1400
	1401	endif ! end of 'n2cond'
	1402
[3412]	1403	if (conservn2) then
	1404	write(,) 'conservn2 n2cond'
	1405	call testconservmass(ngrid,nlayer,pplev(:,1)+ &
	1406	pdpsrf(:)ptimestep,qsurf(:,1)+dqsurf(:,1)ptimestep)
	1407	endif
[3539]	1408	if (methane.and.conservch4) then
	1409	write(,) 'conservch4 n2cond'
	1410	if (fast) then
	1411	call testconservfast(ngrid,nlayer,nq,pq(:,1,igcm_ch4_gas),pdq(:,1,igcm_ch4_gas), &
	1412	qsurf(:,igcm_ch4_ice),dqsurf(:,igcm_ch4_ice), &
	1413	ptimestep,pplev,zdqch4fast,zdqsch4fast,'CH4',' n2cond')
	1414	else
	1415	call testconserv(ngrid,nlayer,nq,pq,pdq,qsurf,dqsurf, &
	1416	igcm_ch4_gas,igcm_ch4_ice, &
	1417	ptimestep,pplev,zdqc,zdqsc,'CH4',' n2cond')
	1418	endif
	1419	endif
[3184]	1420
	1421	!---------------------------------------------
[3195]	1422	! VI. Specific parameterizations for tracers
[3184]	1423	!---------------------------------------------
	1424
	1425	if (tracer) then
[3195]	1426
[3237]	1427	! ---------------------------------------
	1428	! Methane ice condensation in the atmosphere
	1429	! ----------------------------------------
[3374]	1430	rice_ch4(:,:)=0 ! initialization needed for callsedim
[3237]	1431	zdqch4cloud(:,:,:)=0.
	1432	if ((methane).and.(metcloud).and.(.not.fast)) THEN
[3585]	1433	call ch4cloud(ngrid,nlayer,ptimestep, &
[3237]	1434	pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, &
	1435	pq,pdq,zdqch4cloud,zdqsch4cloud,zdtch4cloud, &
	1436	nq,rice_ch4)
	1437
	1438	DO l=1,nlayer
	1439	DO ig=1,ngrid
	1440	pdq(ig,l,igcm_ch4_gas)=pdq(ig,l,igcm_ch4_gas)+ &
	1441	zdqch4cloud(ig,l,igcm_ch4_gas)
	1442	pdq(ig,l,igcm_ch4_ice)=pdq(ig,l,igcm_ch4_ice)+ &
	1443	zdqch4cloud(ig,l,igcm_ch4_ice)
	1444	ENDDO
	1445	ENDDO
	1446
	1447	! Increment methane ice surface tracer tendency
	1448	DO ig=1,ngrid
	1449	dqsurf(ig,igcm_ch4_ice)=dqsurf(ig,igcm_ch4_ice)+ &
	1450	zdqsch4cloud(ig,igcm_ch4_ice)
	1451	ENDDO
	1452
	1453	! update temperature tendancy
	1454	DO ig=1,ngrid
	1455	DO l=1,nlayer
	1456	pdt(ig,l)=pdt(ig,l)+zdtch4cloud(ig,l)
	1457	ENDDO
	1458	ENDDO
	1459	end if
	1460
	1461	! ---------------------------------------
	1462	! CO ice condensation in the atmosphere
	1463	! ----------------------------------------
	1464	zdqcocloud(:,:,:)=0.
	1465	IF ((carbox).and.(monoxcloud).and.(.not.fast)) THEN
[3585]	1466	call cocloud(ngrid,nlayer,ptimestep, &
[3237]	1467	pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, &
	1468	pq,pdq,zdqcocloud,zdqscocloud,zdtcocloud, &
	1469	nq,rice_co,qsurf(1,igcm_n2),dqsurf(1,igcm_n2))
	1470
	1471	DO l=1,nlayer
	1472	DO ig=1,ngrid
	1473	pdq(ig,l,igcm_co_gas)=pdq(ig,l,igcm_co_gas)+ &
	1474	zdqcocloud(ig,l,igcm_co_gas)
	1475	pdq(ig,l,igcm_co_ice)=pdq(ig,l,igcm_co_ice)+ &
	1476	zdqcocloud(ig,l,igcm_co_ice)
	1477	ENDDO
	1478	ENDDO
	1479
	1480	! Increment CO ice surface tracer tendency
	1481	DO ig=1,ngrid
	1482	dqsurf(ig,igcm_co_ice)=dqsurf(ig,igcm_co_ice)+ &
	1483	zdqscocloud(ig,igcm_co_ice)
	1484	ENDDO
	1485
	1486	! update temperature tendancy
	1487	DO ig=1,ngrid
	1488	DO l=1,nlayer
	1489	pdt(ig,l)=pdt(ig,l)+zdtcocloud(ig,l)
	1490	ENDDO
	1491	ENDDO
	1492	ELSE
	1493	rice_co(:,:)=0 ! initialization needed for callsedim
	1494	END IF ! of IF (carbox)
	1495
[3557]	1496	! ----------------------------------------
	1497	! VI.1. Microphysics / Aerosol particles
	1498	! ----------------------------------------
	1499	! Call of microphysics
	1500	IF (callmufi) THEN
[3237]	1501
[3572]	1502	! Production for microphysics
	1503	IF (call_haze_prod_pCH4) THEN
	1504	zdqphot_prec(:,:) = 0.
	1505	zdqphot_ch4(:,:) = 0.
	1506	pdqmufi_prod(:,:,:) = 0.
	1507	call hazecloud(ngrid,nlayer,nq,ptimestep, &
	1508	pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,pdqmufi_prod, &
	1509	zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
[3601]	1510	pdq(:,:,:) = pdq(:,:,:) + pdqmufi_prod(:,:,:) ! Should be updated
[3572]	1511	ENDIF ! end call_haze_prod_pCH4
	1512
[3557]	1513	pdqmufi(:,:,:) = 0.
	1514
	1515	call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,pdqmufi_prod,pdqmufi)
	1516
	1517	pdq(:,:,:) = pdq(:,:,:) + pdqmufi(:,:,:)
[3590]	1518
[3572]	1519	ELSE
	1520	IF (haze) THEN
	1521	zdqphot_prec(:,:) = 0.
	1522	zdqphot_ch4(:,:) = 0.
	1523	zdqhaze(:,:,:) = 0.
	1524	! Forcing to a fixed haze profile if haze_proffix
	1525	if (haze_proffix.and.i_haze.gt.0.) then
	1526	call haze_prof(ngrid,nlayer,zzlay,pplay,pt, &
	1527	reffrad,profmmr)
	1528	zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze))/ptimestep
	1529	else
	1530	call hazecloud(ngrid,nlayer,nq,ptimestep, &
	1531	pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,zdqhaze, &
	1532	zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
	1533	endif
	1534	pdq(:,:,:) = pdq(:,:,:) + zdqhaze(:,:,:) ! Should be updated
	1535	ENDIF ! end haze
[3237]	1536
[3572]	1537	IF (fast.and.fasthaze) THEN
	1538	call prodhaze(ngrid,nlayer,nq,ptimestep,pplev,pq,pdq,dist_star, &
	1539	mu0,declin,zdqprodhaze,zdqsprodhaze,gradflux,fluxbot, &
	1540	fluxlym_sol_bot,fluxlym_ipm_bot,flym_sol,flym_ipm)
	1541	DO ig=1,ngrid
	1542	pdq(ig,1,igcm_ch4_gas)=pdq(ig,1,igcm_ch4_gas)+ &
	1543	zdqprodhaze(ig,igcm_ch4_gas)
	1544	pdq(ig,1,igcm_prec_haze)=pdq(ig,1,igcm_prec_haze)+ &
	1545	zdqprodhaze(ig,igcm_prec_haze)
	1546	pdq(ig,1,igcm_haze)=abs(pdq(ig,1,igcm_haze)+ &
	1547	zdqprodhaze(ig,igcm_haze))
	1548	qsurf(ig,igcm_haze)= qsurf(ig,igcm_haze)+ &
	1549	zdqsprodhaze(ig)*ptimestep
	1550	ENDDO
	1551	ENDIF ! end fast.and.fasthaze
[3237]	1552
[3572]	1553	! Sedimentation.
	1554	if (sedimentation) then
	1555	zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
	1556	zdqssed(1:ngrid,1:nq) = 0.0
	1557	if (oldplutosedim)then
	1558	call callsedim_pluto(ngrid,nlayer,ptimestep, &
	1559	pplev,zzlev,pt,pdt,rice_ch4,rice_co, &
	1560	pq,pdq,zdqsed,zdqssed,nq,pphi)
	1561	else
	1562	call callsedim(ngrid,nlayer,ptimestep, &
	1563	pplev,zzlev,pt,pdt,pq,pdq, &
	1564	zdqsed,zdqssed,nq,pphi)
	1565	endif
	1566	! Whether it falls as rain or snow depends only on the surface temperature
	1567	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
	1568	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
	1569	end if ! end of 'sedimentation'
[3237]	1570
[3572]	1571	ENDIF ! end callmufi
[3237]	1572
[3184]	1573	! ---------------
[3557]	1574	! VI.2. Updates
[3184]	1575	! ---------------
	1576
	1577	! Updating Atmospheric Mass and Tracers budgets.
	1578	if(mass_redistrib) then
	1579
	1580	zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * 0
	1581	! ( zdqevap(1:ngrid,1:nlayer) &
[3275]	1582	! ! + zdqrain(1:ngrid,1:nlayer,igcm_h2o_gas) &
	1583	! ! + dqmoist(1:ngrid,1:nlayer,igcm_h2o_gas) &
[3184]	1584	! + dqvaplscale(1:ngrid,1:nlayer) )
	1585
	1586	do ig = 1, ngrid
	1587	zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
	1588	enddo
[3195]	1589
[3184]	1590	! call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
	1591	! call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
	1592	call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
	1593
	1594	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
	1595	capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
	1596	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
	1597	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
[3195]	1598
[3184]	1599	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
	1600	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
	1601	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
	1602	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
	1603	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
	1604	pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
	1605	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
[3195]	1606
[3184]	1607	endif
	1608
[3275]	1609	! call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
[3184]	1610
	1611	! -----------------------------
[3557]	1612	! VI.3. Surface Tracer Update
[3184]	1613	! -----------------------------
	1614
	1615	qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
	1616
	1617	endif! end of if 'tracer'
	1618
[3241]	1619	if (conservn2) then
	1620	write(,) 'conservn2 tracer'
[3412]	1621	call testconservmass(ngrid,nlayer,pplev(:,1)+ &
	1622	pdpsrf(:)*ptimestep,qsurf(:,1))
[3241]	1623	endif
[3184]	1624
[3241]	1625	DO ig=1,ngrid
	1626	flusurf(ig,igcm_n2)=(qsurf(ig,igcm_n2)- &
	1627	qsurf1(ig,igcm_n2))/ptimestep
	1628	flusurfold(ig,igcm_n2)=flusurf(ig,igcm_n2)
	1629	if (methane) then
	1630	flusurf(ig,igcm_ch4_ice)=(qsurf(ig,igcm_ch4_ice)- &
	1631	qsurf1(ig,igcm_ch4_ice))/ptimestep
	1632	flusurfold(ig,igcm_ch4_ice)=flusurf(ig,igcm_ch4_ice)
	1633	endif
	1634	if (carbox) then
	1635	flusurf(ig,igcm_co_ice)=(qsurf(ig,igcm_co_ice)- &
	1636	qsurf1(ig,igcm_co_ice))/ptimestep
	1637	!flusurfold(ig,igcm_co_ice)=flusurf(ig,igcm_co_ice)
	1638	endif
	1639	ENDDO
	1640
	1641	!! Special source of haze particle !
	1642	! todo: should be placed in haze and use tendency of n2 instead of flusurf
	1643	IF (source_haze) THEN
[3329]	1644	write(,) "Source haze not supported yet."
	1645	stop
[3241]	1646	! call hazesource(ngrid,nlayer,nq,ptimestep, &
	1647	! pplev,flusurf,mu0,zdq_source)
	1648
	1649	DO iq=1, nq
	1650	DO l=1,nlayer
	1651	DO ig=1,ngrid
	1652	pdq(ig,l,iq)=pdq(ig,l,iq)+zdq_source(ig,l,iq)
	1653	ENDDO
	1654	ENDDO
	1655	ENDDO
	1656	ENDIF
	1657
[3184]	1658	!------------------------------------------------
[3195]	1659	! VII. Surface and sub-surface soil temperature
[3184]	1660	!------------------------------------------------
	1661
[3241]	1662	! For diagnostic
	1663	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3539]	1664	if (.not.fast) then
	1665	DO ig=1,ngrid
[3241]	1666	rho(ig,1) = pplay(ig,1)/(r*pt(ig,1))
	1667	sensiblehf1(ig)=rho(ig,1)cpp(0.4/log(zzlay(ig,1)/z0))*2 &
	1668	(pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1))*0.5 &
	1669	(tsurf(ig)-pt(ig,1))
[3452]	1670	if (calldifv) then
	1671	sensiblehf2(ig)=zflubid(ig)-capcal(ig)*zdtsdif(ig)
	1672	end if
[3539]	1673	ENDDO
	1674	endif
[3184]	1675
[3241]	1676
[3412]	1677	! VII.1 Increment surface temperature
	1678	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1679	tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
[3184]	1680
[3412]	1681	! Prevent surface (.e.g. non volatile ch4) to exceed max temperature
	1682	! Lellouch et al., 2000,2011
	1683	IF (tsurfmax) THEN
	1684	DO ig=1,ngrid
	1685	if (albedo_equivalent(ig).gt.albmin_ch4.and. &
	1686	qsurf(ig,igcm_n2).eq.0.) then
	1687	tsurf(ig)=min(tsurf(ig),54.)
	1688	endif
	1689	ENDDO
	1690	ENDIF
	1691
	1692	! VII.2 Compute soil temperatures and subsurface heat flux.
	1693	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3184]	1694	if (callsoil) then
[3539]	1695	call soil(ngrid,nsoilmx,.false.,lastcall,therm_inertia, &
[3195]	1696	ptimestep,tsurf,tsoil,capcal,fluxgrd)
[3184]	1697	endif
	1698
[3412]	1699	! ! For output :
	1700	! tidat_out(:,:)=0.
	1701	! DO l=1,min(nlayermx,nsoilmx)
[3483]	1702	! tidat_out(:,l)=inertiedat(:,l)
[3412]	1703	! ENDDO
[3184]	1704
	1705	! Test energy conservation
	1706	if(enertest)then
[3195]	1707	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf(:)/totarea_planet,dEtots)
[3184]	1708	if (is_master) print*,'Surface energy change =',dEtots,' W m-2'
	1709	endif
	1710
	1711
[3412]	1712
	1713	! VII.3 multiply tendencies of cond/subli for paleo loop only in the
	1714	! last Pluto year of the simulation
	1715	! Year day must be adapted in the startfi for each object
	1716	! Paleo uses year_day to calculate the annual mean tendancies
	1717	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1718	IF (paleo) then
	1719	if (zday.gt.day_ini+ptime0+nday-year_day) then
	1720	DO iq=1,nq
	1721	DO ig=1,ngrid
	1722	qsurfyear(ig,iq)=qsurfyear(ig,iq)+ &
	1723	(qsurf(ig,iq)-qsurf1(ig,iq)) !kg m-2 !ptimestep
	1724	ENDDO
	1725	ENDDO
	1726	endif
	1727	endif
	1728
	1729	! VII.4 Glacial flow at each timestep glastep or at lastcall
	1730	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1731	IF (fast.and.glaflow) THEN
	1732	if ((mod(zday-day_ini-ptime0,glastep)).lt.1. &
	1733	.or.lastcall) then
	1734	IF (lastcall) then
	1735	dstep=mod(zday-day_ini-ptime0,glastep)*daysec
	1736	else
	1737	dstep=glastep*daysec
	1738	endif
	1739	zdqflow(:,:)=qsurf(:,:)
	1740	IF (paleo) then
	1741	call spreadglacier_paleo(ngrid,nq,qsurf, &
	1742	phisfinew,dstep,tsurf)
	1743	else
	1744	call spreadglacier_simple(ngrid,nq,qsurf,dstep)
	1745	endif
	1746	zdqflow(:,:)=(zdqflow(:,:)-qsurf(:,:))/dstep
	1747
	1748	if (conservn2) then
	1749	write(,) 'conservn2 glaflow'
	1750	call testconservmass(ngrid,nlayer,pplev(:,1)+ &
	1751	pdpsrf(:)*ptimestep,qsurf(:,1))
	1752	endif
	1753
	1754	endif
	1755	ENDIF
	1756
[3184]	1757	!---------------------------------------------------
	1758	! VIII. Perform diagnostics and write output files
	1759	!---------------------------------------------------
	1760
	1761	! Note : For output only: the actual model integration is performed in the dynamics.
	1762
	1763
	1764	! Temperature, zonal and meridional winds.
	1765	zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
	1766	zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
	1767	zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
[3195]	1768
[3184]	1769	!! Recast thermal plume vertical velocity array for outputs
	1770	!! AF24: removed
[3195]	1771
[3184]	1772	! Diagnostic.
	1773	zdtdyn(1:ngrid,1:nlayer) = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep
	1774	ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
	1775
	1776	zdudyn(1:ngrid,1:nlayer) = (pu(1:ngrid,1:nlayer)-zuprevious(1:ngrid,1:nlayer)) / ptimestep
	1777	zuprevious(1:ngrid,1:nlayer) = zu(1:ngrid,1:nlayer)
	1778
	1779	if(firstcall)then
	1780	zdtdyn(1:ngrid,1:nlayer)=0.0
	1781	zdudyn(1:ngrid,1:nlayer)=0.0
	1782	endif
	1783
[3539]	1784	! Dynamical heating diagnostic
	1785	fluxdyn(:)=0.0
	1786	if (.not.fast) then
	1787	do ig=1,ngrid
	1788	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp
	1789	enddo
	1790	endif
[3184]	1791
	1792	! Tracers.
	1793	zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
	1794
	1795	! Surface pressure.
	1796	ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
[3539]	1797	call planetwide_sumval(ps(:)*cell_area(:)/totarea_planet,globave)
[3184]	1798
[3243]	1799	! pressure density !pluto specific
	1800	IF (.not.fast) then !
	1801	do ig=1,ngrid
	1802	do l=1,nlayer
	1803	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
	1804	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
	1805	rho(ig,l) = zplay(ig,l)/(r*zt(ig,l))
	1806	enddo
	1807	zplev(ig,nlayer+1)=pplev(ig,nlayer+1)/pplev(ig,1)*ps(ig)
	1808	enddo
	1809	ENDIF
[3184]	1810
[3243]	1811
[3184]	1812	! Surface and soil temperature information
	1813	call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
	1814	call planetwide_minval(tsurf(:),Ts2)
	1815	call planetwide_maxval(tsurf(:),Ts3)
	1816	if(callsoil)then
	1817	TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
	1818	if (is_master) then
	1819	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
	1820	print*,Ts1,Ts2,Ts3,TsS
	1821	end if
	1822	else
	1823	if (is_master) then
	1824	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
	1825	print*,Ts1,Ts2,Ts3
	1826	endif
	1827	end if
	1828
	1829
	1830	! Check the energy balance of the simulation during the run
	1831	if(corrk)then
	1832
	1833	call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
	1834	call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
	1835	call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
	1836	call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
	1837	call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
	1838	do ig=1,ngrid
	1839	if(fluxtop_dn(ig).lt.0.0)then
	1840	print*,'fluxtop_dn has gone crazy'
	1841	print*,'fluxtop_dn=',fluxtop_dn(ig)
	1842	print*,'tau_col=',tau_col(ig)
[3585]	1843	print*,'dtau_aer=',dtau_aer(ig,:,:)
[3184]	1844	print*,'temp= ',pt(ig,:)
	1845	print*,'pplay= ',pplay(ig,:)
	1846	call abort
	1847	endif
	1848	end do
[3195]	1849
[3184]	1850	if(ngrid.eq.1)then
	1851	DYN=0.0
	1852	endif
[3195]	1853
[3184]	1854	if (is_master) then
	1855	print*,' ISR ASR OLR GND DYN [W m^-2]'
	1856	print*, ISR,ASR,OLR,GND,DYN
	1857	endif
	1858
	1859	if(enertest .and. is_master)then
	1860	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
	1861	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
	1862	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
	1863	endif
	1864
	1865	if(meanOLR .and. is_master)then
	1866	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
	1867	! to record global radiative balance
	1868	open(92,file="rad_bal.out",form='formatted',position='append')
	1869	write(92,*) zday,ISR,ASR,OLR
	1870	close(92)
	1871	open(93,file="tem_bal.out",form='formatted',position='append')
	1872	if(callsoil)then
	1873	write(93,*) zday,Ts1,Ts2,Ts3,TsS
	1874	else
	1875	write(93,*) zday,Ts1,Ts2,Ts3
	1876	endif
	1877	close(93)
	1878	endif
	1879	endif
	1880
	1881	endif ! end of 'corrk'
	1882
	1883
	1884	! Diagnostic to test radiative-convective timescales in code.
	1885	if(testradtimes)then
	1886	open(38,file="tau_phys.out",form='formatted',position='append')
	1887	ig=1
	1888	do l=1,nlayer
	1889	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
	1890	enddo
	1891	close(38)
	1892	print*,'As testradtimes enabled,'
	1893	print*,'exiting physics on first call'
	1894	call abort
	1895	endif
	1896
	1897
	1898	! Compute column amounts (kg m-2) if tracers are enabled.
	1899	if(tracer)then
	1900	qcol(1:ngrid,1:nq)=0.0
	1901	do iq=1,nq
	1902	do ig=1,ngrid
	1903	qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
	1904	enddo
	1905	enddo
	1906
	1907	endif ! end of 'tracer'
	1908
[3241]	1909	if (methane) then
	1910	IF (fast) then ! zq is the mixing ratio supposingly mixed in all atmosphere
	1911	DO ig=1,ngrid
	1912	vmr_ch4(ig)=zq(ig,1,igcm_ch4_gas)* &
	1913	mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.
	1914	ENDDO
	1915	ELSE
	1916	DO ig=1,ngrid
	1917	! compute vmr methane
	1918	vmr_ch4(ig)=qcol(ig,igcm_ch4_gas)* &
	1919	g/ps(ig)mmol(igcm_n2)/mmol(igcm_ch4_gas)100.
	1920	! compute density methane
	1921	DO l=1,nlayer
	1922	zrho_ch4(ig,l)=zq(ig,l,igcm_ch4_gas)*rho(ig,l)
	1923	ENDDO
	1924	ENDDO
	1925	ENDIF
[3501]	1926	endif
[3243]	1927
[3501]	1928	if (carbox) then
[3241]	1929	IF (fast) then
	1930	DO ig=1,ngrid
	1931	vmr_co(ig)=zq(ig,1,igcm_co_gas)* &
	1932	mmol(igcm_n2)/mmol(igcm_co_gas)*100.
	1933	ENDDO
	1934	ELSE
	1935	DO ig=1,ngrid
	1936	! compute vmr CO
	1937	vmr_co(ig)=qcol(ig,igcm_co_gas)* &
	1938	g/ps(ig)mmol(igcm_n2)/mmol(igcm_co_gas)100.
	1939	! compute density CO
	1940	DO l=1,nlayer
	1941	zrho_co(ig,l)=zq(ig,l,igcm_co_gas)*rho(ig,l)
	1942	ENDDO
	1943	ENDDO
	1944	ENDIF
[3501]	1945	endif
[3243]	1946
[3501]	1947	zrho_haze(:,:)=0.
	1948	zdqrho_photprec(:,:)=0.
[3572]	1949	IF (haze.and.optichaze) then
[3501]	1950	DO ig=1,ngrid
	1951	DO l=1,nlayer
	1952	zrho_haze(ig,l)=zq(ig,l,igcm_haze)*rho(ig,l)
	1953	zdqrho_photprec(ig,l)=zdqphot_prec(ig,l)*rho(ig,l)
	1954	ENDDO
[3241]	1955	ENDDO
	1956	ENDIF
[3243]	1957
[3501]	1958	IF (fasthaze) then
	1959	DO ig=1,ngrid
	1960	qcol(ig,igcm_haze)=zq(ig,1,igcm_haze)*pplev(ig,1)/g
	1961	qcol(ig,igcm_prec_haze)=zq(ig,1,igcm_prec_haze)*pplev(ig,1)/g
	1962	ENDDO
	1963	ENDIF
[3243]	1964
[3241]	1965	! Info about Ls, declin...
[3501]	1966	IF (fast) THEN
[3241]	1967	if (is_master) write(,),'Ls=',zls180./pi,' dec=',declin180./pi
	1968	if (is_master) write(,),'zday=',zday,' ps=',globave
[3501]	1969	IF (lastcall) then
	1970	if (is_master) write(,),'lastcall'
	1971	ENDIF
	1972	ELSE
[3241]	1973	if (is_master) write(,),'Ls=',zls180./pi,'decli=',declin180./pi,'zday=',zday
[3501]	1974	ENDIF
[3243]	1975
[3501]	1976	lecttsoil=0 ! default value for lecttsoil
	1977	call getin_p("lecttsoil",lecttsoil)
	1978	IF (lastcall.and.(ngrid.EQ.1).and.(lecttsoil.eq.1)) THEN
	1979	! save tsoil temperature profile for 1D profile
	1980	OPEN(13,file='proftsoil.out',form='formatted')
	1981	DO i=1,nsoilmx
	1982	write(13,*) tsoil(1,i)
	1983	ENDDO
	1984	CLOSE(13)
	1985	ENDIF
[3243]	1986
[3184]	1987	if (is_master) print,'--> Ls =',zls180./pi
[3195]	1988
[3184]	1989	if(lastcall) then
[3544]	1990	IF (grid_type==unstructured) THEN !IF DYNAMICO
	1991	! DYNAMICO: no need to add a dynamics time step to ztime_fin
	1992	ztime_fin = ptime
	1993	ELSE ! LMDZ
	1994	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
	1995	ENDIF ! of IF (grid_type==unstructured)
[3184]	1996	!! Update surface ice distribution to iterate to steady state if requested
	1997	!! AF24: removed
[3195]	1998
[3184]	1999	! endif
[3241]	2000	if (paleo) then
	2001	! time range for tendencies of ice flux qsurfyear
	2002	zdt_tot=year_day ! Last year of simulation
	2003
	2004	masslost(:)=0.
	2005	massacc(:)=0.
	2006
	2007	DO ig=1,ngrid
	2008	! update new reservoir of ice on the surface
	2009	DO iq=1,nq
	2010	! kg/m2 to be sublimed or condensed during paleoyears
	2011	qsurfyear(ig,iq)=qsurfyear(ig,iq)* &
	2012	paleoyears365.25/(zdt_totdaysec/86400.)
	2013
	2014	! special case if we sublime the entire reservoir
[3243]	2015	!! AF: TODO : fix following lines (real_area), using line below:
	2016	! call planetwide_sumval((-qsurfyear(:,iq)-qsurf(:,iq))*cell_area(:),masslost)
[3241]	2017
[3243]	2018	! IF (-qsurfyear(ig,iq).gt.qsurf(ig,iq)) THEN
	2019	! masslost(iq)=masslost(iq)+real_area(ig)* &
	2020	! (-qsurfyear(ig,iq)-qsurf(ig,iq))
	2021	! qsurfyear(ig,iq)=-qsurf(ig,iq)
	2022	! ENDIF
[3241]	2023
[3243]	2024	! IF (qsurfyear(ig,iq).gt.0.) THEN
	2025	! massacc(iq)=massacc(iq)+real_area(ig)*qsurfyear(ig,iq)
	2026	! ENDIF
	2027
	2028
[3241]	2029	ENDDO
	2030	ENDDO
	2031
	2032	DO ig=1,ngrid
	2033	DO iq=1,nq
	2034	qsurfpal(ig,iq)=qsurf(ig,iq)+qsurfyear(ig,iq)
	2035	IF (qsurfyear(ig,iq).gt.0.) THEN
	2036	qsurfpal(ig,iq)=qsurfpal(ig,iq)- &
	2037	qsurfyear(ig,iq)*masslost(iq)/massacc(iq)
	2038	ENDIF
	2039	ENDDO
	2040	ENDDO
	2041	! Finally ensure conservation of qsurf
	2042	DO iq=1,nq
[3539]	2043	call planetwide_sumval(qsurf(:,iq)*cell_area(:)/totarea_planet,globaveice(iq))
	2044	call planetwide_sumval(qsurfpal(:,iq)*cell_area(:)/totarea_planet,globavenewice(iq))
[3501]	2045	IF (globavenewice(iq).gt.0.) THEN
	2046	qsurfpal(:,iq)=qsurfpal(:,iq)* &
[3241]	2047	globaveice(iq)/globavenewice(iq)
[3501]	2048	ENDIF
[3241]	2049	ENDDO
	2050
	2051	! update new geopotential depending on the ice reservoir
	2052	phisfipal(:)=phisfinew(:)+qsurfpal(:,igcm_n2)*g/1000.
	2053	!phisfipal(ig)=phisfi(ig)
	2054
	2055	if (kbo.or.triton) then ! case of Triton : we do not change the orbital parameters
[3501]	2056	pdaypal=pday ! no increment of pdaypal to keep same evolution of the subsolar point
	2057	eccpal=1.-periastr/((periastr+apoastr)/2.) !no change of ecc
	2058	peri_daypal=peri_day ! no change
	2059	oblipal=obliquit ! no change
	2060	tpalnew=tpal
	2061	adjustnew=adjust
[3241]	2062
	2063	else ! Pluto
[3501]	2064	! update new pday and tpal (Myr) to be set in startfi controle
	2065	pdaypal=int(day_ini+paleoyears*365.25/6.3872)
	2066	tpalnew=tpal+paleoyears*1.e-6 ! Myrs
[3241]	2067
[3501]	2068	! update new N2 ice adjustment (not tested yet on Pluto)
	2069	adjustnew=adjust
[3241]	2070
[3501]	2071	! update milankovitch parameters : obliquity,Lsp,ecc
	2072	call calcmilank(tpalnew,oblipal,peri_daypal,eccpal)
	2073	!peri_daypal=peri_day
	2074	!eccpal=0.009
[3241]	2075	endif
	2076
	2077	if (is_master) write(,) "Paleo peri=",peri_daypal," tpal=",tpalnew
	2078	if (is_master) write(,) "Paleo eccpal=",eccpal," tpal=",tpalnew
	2079
[3501]	2080	!----------------------------------------------------------------------
	2081	! Writing NetCDF file "RESTARTFI" at the end of the run
	2082	!----------------------------------------------------------------------
	2083	! Note: 'restartfi' is stored just before dynamics are stored
	2084	! in 'restart'. Between now and the writing of 'restart',
	2085	! there will have been the itau=itau+1 instruction and
	2086	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	2087	! thus we store for time=time+dtvr
[3241]	2088
[3539]	2089	! create restartfi
[3501]	2090	if (ngrid.ne.1) then
	2091	print*, "physdem1pal not yet implemented"
	2092	stop
	2093	!TODO: import this routine from pluto.old
	2094	! call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq, &
	2095	! ptimestep,pdaypal, &
	2096	! ztime_fin,tsurf,tsoil,emis,q2,qsurfpal, &
[3539]	2097	! cell_area,albedodat,therm_inertia,zmea,zstd,zsig, &
[3501]	2098	! zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal, &
	2099	! peri_daypal)
	2100	endif
	2101	else ! 'paleo'
[3195]	2102
[3501]	2103	if (ngrid.ne.1) then
	2104	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
[3195]	2105
[3501]	2106	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
	2107	ptimestep,ztime_fin, &
[3539]	2108	tsurf,tsoil,therm_inertia,emis,q2,qsurf)
[3501]	2109	endif
[3539]	2110
[3501]	2111	endif ! end of 'paleo'
	2112	endif ! end of 'lastcall'
[3184]	2113
[3501]	2114	!------------------------------------------------------------------------------
	2115	! OUTPUT in netcdf file "DIAGFI.NC",
	2116	! containing any variable for diagnostic
[3184]	2117	!
	2118	! Note 1 : output with period "ecritphy", set in "run.def"
[3501]	2119	! Note 2 : writediagfi can also be called from any other subroutine
	2120	! for any variable, but its preferable to keep all the
	2121	! calls in one place ...
	2122	!------------------------------------------------------------------------------
[3184]	2123
[3539]	2124	!-------- General 1D variables
	2125
[3508]	2126	call write_output("Ls","solar longitude","deg",zls*180./pi)
	2127	! call write_output("Lss","sub solar longitude","deg",zlss*180./pi)
	2128	call write_output("RA","right ascension","deg",right_ascen*180./pi)
	2129	call write_output("Declin","solar declination","deg",declin*180./pi)
	2130	call write_output("dist_star","dist_star","AU",dist_star)
[3539]	2131	call write_output("globave","surf press","Pa",globave)
[3184]	2132
[3539]	2133	!-------- General 2D variables
	2134
[3508]	2135	call write_output("tsurf","Surface temperature","K",tsurf)
	2136	call write_output("ps","Surface pressure","Pa",ps)
	2137	call write_output("emis","Emissivity","",emis)
[3539]	2138	!if (grid_type == regular_lonlat) then
	2139	! call write_output("area","Mesh area","m2", &
	2140	! cell_area_for_lonlat_outputs)
	2141	! else ! unstructured grid (e.g. dynamico)
	2142	! call write_output("area","Mesh area","m2",cell_area)
	2143	!endif
[3241]	2144
	2145	if (fast) then
[3508]	2146	call write_output("fluxrad","fluxrad","W m-2",fluxrad)
	2147	call write_output("fluxgrd","fluxgrd","W m-2",fluxgrd)
[3539]	2148	! call write_output("dplanck","dplanck","W.s m-2 K-1",dplanck)
	2149	! "soil" variables
[3508]	2150	call write_output("capcal","capcal","W.s m-2 K-1",capcal)
	2151	call write_output("tsoil","tsoil","K",tsoil)
[3539]	2152	endif
	2153
	2154	! Total energy balance diagnostics
	2155	if(callrad)then
	2156	call write_output("ALB","Surface albedo"," ",albedo_equivalent)
	2157	call write_output("ASR","absorbed stellar rad.","W m-2",fluxabs_sw)
	2158	call write_output("ISR","incoming stellar rad.","W m-2",fluxtop_dn)
	2159	call write_output("OLR","outgoing longwave rad.","W m-2",fluxtop_lw)
	2160	call write_output("GND","heat flux from ground","W m-2",fluxgrd)
	2161	if (.not.fast) then
	2162	call write_output("DYN","dynamical heat input","W m-2",fluxdyn)
	2163	endif
	2164	endif ! end of 'callrad'
	2165
	2166	!-------- General 3D variables
	2167
	2168	if (.not.fast) then
[3243]	2169	if (check_physics_outputs) then
	2170	! Check the validity of updated fields at the end of the physics step
[3508]	2171	call check_physics_fields("physiq:", zt, zu, zv, pplev, zq)
[3243]	2172	endif
	2173
[3508]	2174	call write_output("zzlay","Midlayer altitude", "m",zzlay(:,:))
	2175	call write_output("zzlev","Interlayer altitude", "m",zzlev(:,1:nlayer))
	2176	!call write_output('pphi','Geopotential',' ',pphi)
	2177
	2178	call write_output("temperature","temperature","K",zt)
	2179	call write_output("teta","potential temperature","K",zh)
	2180	call write_output("u","Zonal wind","m.s-1",zu)
	2181	call write_output("v","Meridional wind","m.s-1",zv)
	2182	call write_output("w","Vertical wind","m.s-1",pw)
	2183	call write_output("p","Pressure","Pa",pplay)
[3522]	2184	call write_output("omega","omega","Pa/s",omega)
[3228]	2185	endif
	2186
[3184]	2187	if(enertest) then
	2188	if (calldifv) then
[3508]	2189	call write_output("q2","turbulent kinetic energy","J.kg^-1",q2)
	2190	call write_output("sensibFlux","sensible heat flux","w.m^-2",sensibFlux)
[3184]	2191	endif
[3195]	2192
[3184]	2193	if (corrk) then
[3508]	2194	call write_output("dEzradsw","radiative heating","w.m^-2",dEzradsw)
	2195	call write_output("dEzradlw","radiative heating","w.m^-2",dEzradlw)
[3184]	2196	endif
	2197	endif ! end of 'enertest'
	2198
[3539]	2199	! Diagnostics of optical thickness
	2200	! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
	2201	if (diagdtau) then
[3501]	2202	do nw=1,L_NSPECTV
	2203	write(str2,'(i2.2)') nw
[3508]	2204	call write_output('dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',int_dtauv(:,nlayer:1:-1,nw))
[3501]	2205	enddo
	2206	do nw=1,L_NSPECTI
	2207	write(str2,'(i2.2)') nw
[3508]	2208	call write_output('dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',int_dtaui(:,nlayer:1:-1,nw))
[3501]	2209	enddo
[3539]	2210	endif
[3184]	2211
[3539]	2212	! Temporary inclusions for heating diagnostics.
	2213	if (.not.fast) then
	2214	call write_output("zdtsw","SW heating","T s-1",zdtsw)
	2215	call write_output("zdtlw","LW heating","T s-1",zdtlw)
	2216	call write_output("dtrad","radiative heating","K s-1",dtrad)
	2217	call write_output("zdtdyn","Dyn. heating","T s-1",zdtdyn)
	2218	endif
[3195]	2219
[3539]	2220	! For Debugging.
	2221	!call write_output('rnat','Terrain type',' ',real(rnat))
[3184]	2222
	2223	! Output tracers.
	2224	if (tracer) then
[3195]	2225
[3184]	2226	do iq=1,nq
[3539]	2227	if (.not.fast) then
	2228	call write_output(noms(iq),noms(iq),'kg/kg',zq(:,:,iq))
	2229	endif
[3522]	2230	call write_output(trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
[3508]	2231	'kg m^-2',qcol(:,iq) )
	2232	call write_output(trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
	2233	'kg m^-2',qsurf(:,iq) )
[3184]	2234	enddo ! end of 'nq' loop
[3195]	2235
[3539]	2236	! N2 cycle
	2237	call write_output('n2_iceflux','n2_iceflux',"kg m^-2 s^-1",flusurf(:,igcm_n2) )
	2238	if (.not.fast) then
	2239	call write_output("zdtc","tendancy T cond N2","K",zdtc)
	2240	endif
	2241
	2242	! CH4 cycle
[3241]	2243	if (methane) then
[3453]	2244
[3508]	2245	call write_output('ch4_iceflux','ch4_iceflux',&
[3539]	2246	"kg m^-2 s^-1",flusurf(:,igcm_ch4_ice) )
[3508]	2247	call write_output("vmr_ch4","vmr_ch4","%",vmr_ch4)
[3539]	2248
[3241]	2249	if (.not.fast) then
[3508]	2250	call write_output("zrho_ch4","zrho_ch4","kg.m-3",zrho_ch4(:,:))
[3539]	2251	!call write_output("rice_ch4","ch4 ice mass mean radius","m",rice_ch4)
	2252	!call write_output("zq1temp_ch4"," "," ",zq1temp_ch4)
	2253	!call write_output("qsat_ch4"," "," ",qsat_ch4)
	2254	!call write_output("qsat_ch4_l1"," "," ",qsat_ch4_l1)
[3243]	2255
[3539]	2256	! 3D Tendancies
	2257	call write_output("zdqcn2_ch4","zdq condn2 ch4","",&
[3508]	2258	zdqc(:,:,igcm_ch4_gas))
[3539]	2259	call write_output("zdqdif_ch4","zdqdif ch4","",&
[3508]	2260	zdqdif(:,:,igcm_ch4_gas))
[3539]	2261	call write_output("zdqsdif_ch4_ice","zdqsdif ch4","",&
[3508]	2262	zdqsdif(:,igcm_ch4_ice))
[3539]	2263	call write_output("zdqadj_ch4","zdqadj ch4","",&
[3508]	2264	zdqadj(:,:,igcm_ch4_gas))
[3539]	2265	endif
	2266
[3501]	2267	if (sedimentation) then
[3508]	2268	call write_output("zdqsed_ch4","zdqsed ch4","",&
	2269	zdqsed(:,:,igcm_ch4_gas))
	2270	call write_output("zdqssed_ch4","zdqssed ch4","",&
	2271	zdqssed(:,igcm_ch4_gas))
[3501]	2272	endif
[3539]	2273
[3501]	2274	if (metcloud.and.(.not.fast)) then
[3508]	2275	call write_output("zdtch4cloud","ch4 cloud","T s-1",&
	2276	zdtch4cloud)
	2277	call write_output("zdqch4cloud","ch4 cloud","T s-1",&
[3539]	2278	zdqch4cloud(:,:,igcm_ch4_gas))
[3439]	2279	endif
[3243]	2280
[3241]	2281	endif
[3243]	2282
[3539]	2283	! CO cycle
[3241]	2284	if (carbox) then
[3508]	2285	! call write_output("zdtcocloud","tendancy T cocloud","K",zdtcocloud)
	2286	call write_output('co_iceflux','co_iceflux',&
[3539]	2287	"kg m^-2 s^-1",flusurf(:,igcm_co_ice) )
[3508]	2288	call write_output("vmr_co","vmr_co","%",vmr_co)
[3241]	2289	if (.not.fast) THEN
[3508]	2290	call write_output("zrho_co","zrho_co","kg.m-3",zrho_co(:,:))
[3241]	2291	endif
	2292	endif
[3243]	2293
[3539]	2294	! Haze
[3241]	2295	if (haze) then
[3539]	2296
	2297	if (haze_radproffix)then
	2298	call write_output('haze_reff','haze_reff','m',reffrad(:,:,1))
	2299	end if
	2300	!call write_output("zrho_haze","zrho_haze","kg.m-3",zrho_haze(:,:))
	2301	!call write_output("zdqhaze_col","zdqhaze col","kg/m2/s",&
	2302	! zdqhaze_col(:))
	2303
	2304	! 3D Tendencies
[3508]	2305	call write_output("zdqrho_photprec","zdqrho_photprec",&
	2306	"kg.m-3.s-1",zdqrho_photprec(:,:))
	2307	call write_output("zdqphot_prec","zdqphot_prec","",&
	2308	zdqphot_prec(:,:))
	2309	call write_output("zdqhaze_ch4","zdqhaze_ch4","",&
	2310	zdqhaze(:,:,igcm_ch4_gas))
	2311	call write_output("zdqhaze_prec","zdqhaze_prec","",&
	2312	zdqhaze(:,:,igcm_prec_haze))
[3539]	2313	call write_output("zdqphot_ch4","zdqphot_ch4","",&
	2314	zdqphot_ch4(:,:))
	2315	call write_output("zdqconv_prec","zdqconv_prec","",&
	2316	zdqconv_prec(:,:))
	2317
[3241]	2318	if (igcm_haze.ne.0) then
[3508]	2319	call write_output("zdqhaze_haze","zdqhaze_haze","",&
	2320	zdqhaze(:,:,igcm_haze))
[3243]	2321	if (sedimentation) then
[3508]	2322	call write_output("zdqssed_haze","zdqssed haze",&
	2323	"kg/m2/s",zdqssed(:,igcm_haze))
[3243]	2324	endif
[3241]	2325	endif
[3243]	2326
[3572]	2327	if (optichaze) then
[3539]	2328	call write_output("tau_col",&
	2329	"Total aerosol optical depth","opacity",tau_col)
	2330	endif
	2331
[3590]	2332	endif ! end haze
[3184]	2333
[3590]	2334	if (callmufi) then
	2335	! Tracers:
	2336	call write_output("m0as","Spherical mode 0th order moment","m-3",zq(:,:,micro_indx(1))*int2ext(:,:))
	2337	call write_output("m3as","Spherical mode 3rd order moment","m3.m-3",zq(:,:,micro_indx(2))*int2ext(:,:))
	2338	call write_output("m0af","Fractal mode 0th order moment","m-3",zq(:,:,micro_indx(3))*int2ext(:,:))
	2339	call write_output("m3af","Fractal mode 3rd order moment","m3.m-3",zq(:,:,micro_indx(4))*int2ext(:,:))
	2340
	2341	! Diagnostics:
	2342	call write_output("rcs","Spherical mode characteristic radius","m",mp2m_rc_sph(:,:))
	2343	call write_output("rcf","Fractal mode characteristic radius","m",mp2m_rc_fra(:,:))
	2344
	2345	if (optichaze) then
	2346	call write_output("tau_col",&
	2347	"Total aerosol optical depth","opacity",tau_col)
	2348	endif ! end optichaze
	2349	endif ! end callmufi
	2350
[3241]	2351	endif ! end of 'tracer'
[3184]	2352
	2353	! Output spectrum.
[3195]	2354	if(specOLR.and.corrk)then
[3184]	2355	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
	2356	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
	2357	call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W/m^2/cm^-1",3,GSR_nu)
	2358	endif
	2359
	2360	! XIOS outputs
[3195]	2361	#ifdef CPP_XIOS
[3184]	2362	! Send fields to XIOS: (NB these fields must also be defined as
	2363	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
[3506]	2364	CALL send_xios_field("controle",tab_cntrl_mod,1)
[3195]	2365
[3506]	2366	CALL send_xios_field("ap",ap,1)
	2367	CALL send_xios_field("bp",bp,1)
	2368	CALL send_xios_field("aps",aps,1)
	2369	CALL send_xios_field("bps",bps,1)
[3195]	2370
[3184]	2371	if (lastcall.and.is_omp_master) then
	2372	write(,) "physiq: call xios_context_finalize"
	2373	call xios_context_finalize
	2374	endif
	2375	#endif
[3195]	2376
[3184]	2377	if (check_physics_outputs) then
	2378	! Check the validity of updated fields at the end of the physics step
	2379	call check_physics_fields("end of physiq:", zt, zu, zv, pplev, zq)
	2380	endif
	2381
	2382	icount=icount+1
[3195]	2383
[3184]	2384	end subroutine physiq
[3195]	2385
[3184]	2386	end module physiq_mod

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

Download in other formats:

Original Format