Context Navigation

physiq_mod.F90 @ 3356

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

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

Download in other formats:

Original Format