Context Navigation

physiq_mod.F90 @ 3508

Last change on this file since 3508 was 3508, checked in by afalco, 2 weeks ago
Pluto: physiq_mod to call write_output, which writes both via XIOS and diagfi routines. Included xml example files for XIOS. AF
Property svn:executable set to ``*
File size: 106.8 KB

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

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

Download in other formats:

Original Format