Context Navigation

physiq_mod.F90 @ 3506

Last change on this file since 3506 was 3506, checked in by afalco, 2 weeks ago
Pluto: import write_output function from Mars. xios specific outputs. AF
Property svn:executable set to ``*
File size: 108.4 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
	102
	103	implicit none
	104
	105
	106	!==================================================================
[3195]	107	!
[3184]	108	! Purpose
	109	! -------
	110	! Central subroutine for all the physics parameterisations in the
	111	! universal model. Originally adapted from the Mars LMDZ model.
	112	!
	113	! The model can be run without or with tracer transport
	114	! depending on the value of "tracer" in file "callphys.def".
	115	!
	116	!
	117	! It includes:
	118	!
	119	! I. Initialization :
	120	! I.1 Firstcall initializations.
	121	! I.2 Initialization for every call to physiq.
	122	!
	123	! II. Compute radiative transfer tendencies (longwave and shortwave) :
	124	! II.a Option 1 : Call correlated-k radiative transfer scheme.
	125	! II.b Option 2 : Call Newtonian cooling scheme.
[3237]	126	! II.! Option 3 : Atmosphere has no radiative effect.
[3184]	127	!
	128	! III. Vertical diffusion (turbulent mixing) :
	129	!
	130	! IV. Convection :
	131	! IV.a Thermal plume model
	132	! IV.b Dry convective adjusment
	133	!
	134	! V. Condensation and sublimation of gases (currently just N2).
	135	!
	136	! VI. Tracers
	137	! VI.1. Water and water ice.
	138	! VI.2 Photochemistry
	139	! VI.3. Aerosols and particles.
	140	! VI.4. Updates (pressure variations, surface budget).
	141	! VI.5. Slab Ocean.
	142	! VI.6. Surface Tracer Update.
	143	!
	144	! VII. Surface and sub-surface soil temperature.
	145	!
	146	! VIII. Perform diagnostics and write output files.
	147	!
	148	!
	149	! arguments
	150	! ---------
	151	!
	152	! INPUT
	153	! -----
	154	!
	155	! ngrid Size of the horizontal grid.
	156	! nlayer Number of vertical layers.
	157	! nq Number of advected fields.
	158	!
	159	! firstcall True at the first call.
	160	! lastcall True at the last call.
	161	!
	162	! pday Number of days counted from the North. Spring equinoxe.
	163	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
	164	! ptimestep timestep (s).
	165	!
	166	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa).
	167	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
	168	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2.s-2).
	169	!
	170	! pu(ngrid,nlayer) u, zonal component of the wind (ms-1).
	171	! pv(ngrid,nlayer) v, meridional component of the wind (ms-1).
	172	!
	173	! pt(ngrid,nlayer) Temperature (K).
	174	!
	175	! pq(ngrid,nlayer,nq) Advected fields.
	176	!
	177	! pudyn(ngrid,nlayer) \
	178	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
	179	! ptdyn(ngrid,nlayer) / corresponding variables.
	180	! pqdyn(ngrid,nlayer,nq) /
	181	! flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
	182	!
	183	! OUTPUT
	184	! ------
	185	!
	186	! pdu(ngrid,nlayer) \
	187	! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding
	188	! pdt(ngrid,nlayer) / variables due to physical processes.
	189	! pdq(ngrid,nlayer) /
	190	! pdpsrf(ngrid) /
	191	!
	192	!
	193	! Authors
	194	! -------
	195	! Frederic Hourdin 15/10/93
	196	! Francois Forget 1994
[3195]	197	! Christophe Hourdin 02/1997
[3184]	198	! Subroutine completely rewritten by F. Forget (01/2000)
	199	! Water ice clouds: Franck Montmessin (update 06/2003)
	200	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
	201	! New correlated-k radiative scheme: R. Wordsworth (2009)
[3195]	202	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
[3184]	203	! Improved water cycle: R. Wordsworth / B. Charnay (2010)
	204	! To F90: R. Wordsworth (2010)
	205	! New turbulent diffusion scheme: J. Leconte (2012)
	206	! Loops converted to F90 matrix format: J. Leconte (2012)
[3195]	207	! No more ngrid/nq, F90 commons and adaptation to parallel: A. Spiga (2012)
[3184]	208	! Purge of the code : M. Turbet (2015)
	209	! Photochemical core developped by F. Lefevre: B. Charnay (2017)
	210	! Purge for Pluto model : A. Falco (2024)
	211	!==================================================================
	212
	213
	214	! 0. Declarations :
	215	! ------------------
	216
	217	include "netcdf.inc"
	218
	219	! Arguments :
	220	! -----------
	221
	222	! INPUTS:
	223	! -------
	224
	225	integer,intent(in) :: ngrid ! Number of atmospheric columns.
	226	integer,intent(in) :: nlayer ! Number of atmospheric layers.
	227	integer,intent(in) :: nq ! Number of tracers.
[3195]	228
[3184]	229	logical,intent(in) :: firstcall ! Signals first call to physics.
	230	logical,intent(in) :: lastcall ! Signals last call to physics.
[3195]	231
[3184]	232	real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0.
	233	real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
	234	real,intent(in) :: ptimestep ! Physics timestep (s).
	235	real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
	236	real,intent(in) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
	237	real,intent(in) :: pphi(ngrid,nlayer) ! Geopotential at mid-layer (m2s-2).
	238	real,intent(in) :: pu(ngrid,nlayer) ! Zonal wind component (m/s).
	239	real,intent(in) :: pv(ngrid,nlayer) ! Meridional wind component (m/s).
	240	real,intent(in) :: pt(ngrid,nlayer) ! Temperature (K).
	241	real,intent(in) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
	242	real,intent(in) :: flxw(ngrid,nlayer) ! Vertical mass flux (ks/s) at lower boundary of layer
	243
	244	! OUTPUTS:
	245	! --------
	246
	247	! Physical tendencies :
	248
	249	real,intent(out) :: pdu(ngrid,nlayer) ! Zonal wind tendencies (m/s/s).
	250	real,intent(out) :: pdv(ngrid,nlayer) ! Meridional wind tendencies (m/s/s).
	251	real,intent(out) :: pdt(ngrid,nlayer) ! Temperature tendencies (K/s).
	252	real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
	253	real,intent(out) :: pdpsrf(ngrid) ! Surface pressure tendency (Pa/s).
	254
	255	! Local saved variables:
	256	! ----------------------
	257	integer,save :: day_ini ! Initial date of the run (sol since Ls=0).
	258	integer,save :: icount ! Counter of calls to physiq during the run.
	259	!$OMP THREADPRIVATE(day_ini,icount)
	260
[3241]	261	!Pluto specific
[3237]	262	REAL,save :: acond,bcond
[3241]	263	REAL,save :: tcond1p4Pa
[3237]	264	DATA tcond1p4Pa/38/
	265
	266
[3241]	267
[3184]	268	! Local variables :
	269	! -----------------
[3241]	270	! Tendencies for the paleoclimate mode
	271	REAL qsurfyear(ngrid,nq) ! kg.m-2 averaged mass of ice lost/gained in the last Pluto year of the run
	272	REAL phisfinew(ngrid) ! geopotential of the bedrock (= phisfi-qsurf/1000*g)
	273	REAL qsurfpal(ngrid,nq) ! qsurf after a paleoclimate step : for physdem1 and restartfi
	274	REAL phisfipal(ngrid) ! geopotential after a paleoclimate step : for physdem1 and restartfi
	275	REAL oblipal ! change of obliquity
	276	REAL peri_daypal ! new periday
	277	REAL eccpal ! change of eccentricity
	278	REAL tpalnew ! change of time
	279	REAL adjustnew ! change in N2 ice albedo
	280	REAL pdaypal ! new pday = day_ini + step
	281	REAL zdt_tot ! time range corresponding to the flux of qsurfyear
	282	REAL massacc(nq) ! accumulated mass
	283	REAL masslost(nq) ! accumulated mass
[3184]	284
[3241]	285	REAL globave ! globalaverage 2D ps
	286	REAL globaveice(nq) ! globalaverage 2D ice
	287	REAL globavenewice(nq) ! globalaverage 2D ice
	288	INTEGER lecttsoil ! lecture of tsoil from proftsoil
	289	REAL qsurf1(ngrid,nq) ! saving qsurf to calculate flux over long timescales kg.m-2
	290	REAL flusurf(ngrid,nq) ! flux cond/sub kg.m-2.s-1
	291	REAL flusurfold(ngrid,nq) ! old flux cond/sub kg.m-2.s-1
[3243]	292	REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
[3241]	293
[3412]	294	REAL,SAVE :: ptime0 ! store the first time
	295	REAL dstep
	296	REAL,SAVE :: glastep=20 ! step in pluto day to spread glacier
[3241]	297
	298
[3195]	299	! Aerosol (dust or ice) extinction optical depth at reference wavelength
[3184]	300	! for the "naerkind" optically active aerosols:
	301
	302	real,save,allocatable :: aerosol(:,:,:) ! Aerosols
	303	!$OMP THREADPRIVATE(aerosol)
	304	real zh(ngrid,nlayer) ! Potential temperature (K).
	305	real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
	306	real omega(ngrid,nlayer) ! omega velocity (Pa/s, >0 when downward)
	307
[3275]	308	integer i,l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_gas, igcm_generic_ice
	309	logical call_ice_gas_generic ! to call only one time the ice/vap pair of a tracer
[3195]	310
[3184]	311	real zls ! Solar longitude (radians).
	312	real zlss ! Sub solar point longitude (radians).
	313	real zday ! Date (time since Ls=0, calculated in sols).
[3329]	314	REAL,save :: saveday ! saved date
	315	REAL,save :: savedeclin ! saved declin
[3184]	316	real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers.
	317	real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries.
	318
[3237]	319
[3184]	320	! VARIABLES for the thermal plume model (AF24: deleted)
[3195]	321
[3184]	322	! TENDENCIES due to various processes :
	323
	324	! For Surface Temperature : (K/s)
	325	real zdtsurf(ngrid) ! Cumulated tendencies.
	326	real zdtsurfmr(ngrid) ! Mass_redistribution routine.
	327	real zdtsurfc(ngrid) ! Condense_n2 routine.
	328	real zdtsdif(ngrid) ! Turbdiff/vdifc routines.
	329	! real zdtsurf_hyd(ngrid) ! Hydrol routine.
[3195]	330
	331	! For Atmospheric Temperatures : (K/s)
[3184]	332	real dtlscale(ngrid,nlayer) ! Largescale routine.
[3195]	333	real dt_generic_condensation(ngrid,nlayer) ! condensation_generic routine.
[3184]	334	real zdtc(ngrid,nlayer) ! Condense_n2 routine.
	335	real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	336	real zdtmr(ngrid,nlayer) ! Mass_redistribution routine.
	337	real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
	338	real zdtchim(ngrid,nlayer) ! Calchim routine.
[3195]	339
[3184]	340	! For Surface Tracers : (kg/m2/s)
	341	real dqsurf(ngrid,nq) ! Cumulated tendencies.
[3228]	342	!real zdqsurfc(ngrid) ! Condense_n2 routine.
[3195]	343	REAL zdqsc(ngrid,nq) ! Condense_n2 routine.
[3184]	344	real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines.
	345	real zdqssed(ngrid,nq) ! Callsedim routine.
	346	real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine.
[3195]	347
[3184]	348	! For Tracers : (kg/kg_of_air/s)
	349	real zdqc(ngrid,nlayer,nq) ! Condense_n2 routine.
	350	real zdqadj(ngrid,nlayer,nq) ! Convadj routine.
	351	real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines.
	352	real zdqevap(ngrid,nlayer) ! Turbdiff routine.
	353	real zdqsed(ngrid,nlayer,nq) ! Callsedim routine.
	354	real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine.
	355	real dqvaplscale_generic(ngrid,nlayer,nq) ! condensation_generic routine.
	356	real dqcldlscale_generic(ngrid,nlayer,nq) ! condensation_generic routine.
	357	REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine
	358	REAL,allocatable,save :: zdqschim(:,:) ! Calchim_asis routine
	359	!$OMP THREADPRIVATE(zdqchim,zdqschim)
	360
[3195]	361
	362	!! PLUTO variables
	363	REAL zdqch4cloud(ngrid,nlayer,nq)
	364	REAL zdqsch4cloud(ngrid,nq)
	365	REAL zdtch4cloud(ngrid,nlayer)
	366	REAL zdqcocloud(ngrid,nlayer,nq)
	367	REAL zdqscocloud(ngrid,nq)
	368	REAL zdtcocloud(ngrid,nlayer)
	369	REAL rice_ch4(ngrid,nlayer) ! Methane ice geometric mean radius (m)
	370	REAL rice_co(ngrid,nlayer) ! CO ice geometric mean radius (m)
	371
	372	REAL zdqsch4fast(ngrid) ! used only for fast model nogcm
	373	REAL zdqch4fast(ngrid) ! used only for fast model nogcm
	374	REAL zdqscofast(ngrid) ! used only for fast model nogcm
	375	REAL zdqcofast(ngrid) ! used only for fast model nogcm
	376	REAL zdqflow(ngrid,nq)
	377
[3455]	378	REAL zdtconduc(ngrid,nlayer) ! (K/s)
	379	REAL zdumolvis(ngrid,nlayer)
	380	REAL zdvmolvis(ngrid,nlayer)
	381	real zdqmoldiff(ngrid,nlayer,nq)
	382
[3195]	383	! Haze relatated tendancies
	384	REAL zdqhaze(ngrid,nlayer,nq)
	385	REAL zdqprodhaze(ngrid,nq)
	386	REAL zdqsprodhaze(ngrid)
	387	REAL zdqhaze_col(ngrid)
	388	REAL zdqphot_prec(ngrid,nlayer)
	389	REAL zdqphot_ch4(ngrid,nlayer)
	390	REAL zdqconv_prec(ngrid,nlayer)
	391	REAL zdq_source(ngrid,nlayer,nq)
	392	! Fast Haze relatated tendancies
	393	REAL fluxbot(ngrid)
	394	REAL gradflux(ngrid)
	395	REAL fluxlym_sol_bot(ngrid) ! Solar flux Lyman alpha ph.m-2.s-1 reaching the surface
	396	REAL fluxlym_ipm_bot(ngrid) ! IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 reaching the surface
	397	REAL flym_sol(ngrid) ! Incident Solar flux Lyman alpha ph.m-2.s-1
	398	REAL flym_ipm(ngrid) ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1
[3237]	399	REAL zfluxuv ! Lyman alpha flux at 1AU
[3195]	400
	401
	402
[3184]	403	REAL array_zero1(ngrid)
	404	REAL array_zero2(ngrid,nlayer)
[3195]	405
[3184]	406	! For Winds : (m/s/s)
	407	real zdvadj(ngrid,nlayer), zduadj(ngrid,nlayer) ! Convadj routine.
	408	real zdumr(ngrid,nlayer), zdvmr(ngrid,nlayer) ! Mass_redistribution routine.
	409	real zdvdif(ngrid,nlayer), zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	410	real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	411	real zdhadj(ngrid,nlayer) ! Convadj routine.
[3195]	412	REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer) ! condense_n2 routine.
	413
[3184]	414	! For Pressure and Mass :
	415	real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
	416	real zdmassmr_col(ngrid) ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
	417	real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s).
	418
[3195]	419
	420
[3184]	421	! Local variables for LOCAL CALCULATIONS:
	422	! ---------------------------------------
	423	real zflubid(ngrid)
	424	real zplanck(ngrid),zpopsk(ngrid,nlayer)
[3237]	425	REAL zdum1(ngrid,nlayer)
	426	REAL zdum2(ngrid,nlayer)
[3184]	427	real ztim1,ztim2,ztim3, z1,z2
	428	real ztime_fin
	429	real zdh(ngrid,nlayer)
	430	real gmplanet
	431	real taux(ngrid),tauy(ngrid)
	432
	433
	434	! local variables for DIAGNOSTICS : (diagfi & stat)
	435	! -------------------------------------------------
	436	real ps(ngrid) ! Surface Pressure.
	437	real zt(ngrid,nlayer) ! Atmospheric Temperature.
	438	real zu(ngrid,nlayer),zv(ngrid,nlayer) ! Zonal and Meridional Winds.
	439	real zq(ngrid,nlayer,nq) ! Atmospheric Tracers.
	440	real zdtadj(ngrid,nlayer) ! Convadj Diagnostic.
	441	real zdtdyn(ngrid,nlayer) ! Dynamical Heating (K/s).
	442	real zdudyn(ngrid,nlayer) ! Dynamical Zonal Wind tendency (m.s-2).
	443
	444	real reff(ngrid,nlayer) ! Effective dust radius (used if doubleq=T).
	445	real vmr(ngrid,nlayer) ! volume mixing ratio
	446	real time_phys
[3195]	447
[3184]	448	real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
[3195]	449
[3184]	450	real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
	451
[3195]	452	! Pluto adding variables
	453	real vmr_ch4(ngrid) ! vmr ch4
	454	real vmr_co(ngrid) ! vmr co
	455	real rho(ngrid,nlayer) ! density
	456	real zrho_ch4(ngrid,nlayer) ! density methane kg.m-3
	457	real zrho_co(ngrid,nlayer) ! density CO kg.m-3
	458	real zrho_haze(ngrid,nlayer) ! density haze kg.m-3
	459	real zdqrho_photprec(ngrid,nlayer) !photolysis rate kg.m-3.s-1
	460	real zq1temp_ch4(ngrid) !
	461	real qsat_ch4(ngrid) !
	462	real qsat_ch4_l1(ngrid) !
	463	! CHARACTER(LEN=20) :: txt ! to temporarly store text for eddy tracers
	464	real profmmr(ngrid,nlayer) ! fixed profile of haze if haze_proffix
	465	real sensiblehf1(ngrid) ! sensible heat flux
	466	real sensiblehf2(ngrid) ! sensible heat flux
	467
[3184]	468	! included by RW for H2O Manabe scheme
	469	real rneb_man(ngrid,nlayer) ! H2O cloud fraction (moistadj).
	470	real rneb_lsc(ngrid,nlayer) ! H2O cloud fraction (large scale).
	471
	472
	473	! to test energy conservation (RW+JL)
	474	real mass(ngrid,nlayer),massarea(ngrid,nlayer)
	475	real dEtot, dEtots, AtmToSurf_TurbFlux
	476	real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
	477	!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
[3195]	478
[3184]	479	!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
	480
[3195]	481	real dtmoist_max,dtmoist_min
[3184]	482	real dItot, dItot_tmp, dVtot, dVtot_tmp
	483
	484
	485	real dWtot, dWtot_tmp, dWtots, dWtots_tmp
	486
	487	real psat_tmp ! AF24: to remove?
	488
	489	real qsat_generic(ngrid,nlayer,nq) ! generic condensable tracers (GCS) specific concentration at saturation (kg/kg_of_air).
	490	real RH_generic(ngrid,nlayer,nq) ! generic condensable tracers (GCS) Relative humidity.
	491	real rneb_generic(ngrid,nlayer,nq) ! GCS cloud fraction (generic condensation).
	492	real psat_tmp_generic
	493	real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
	494	!$OMP THREADPRIVATE(metallicity)
	495
	496	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
	497
	498	! For Clear Sky Case. (AF24: deleted)
	499
	500	real nconsMAX, vdifcncons(ngrid), cadjncons(ngrid) ! Vdfic water conservation test. By RW
	501
	502	real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW
	503
	504	real,save,allocatable :: reffcol(:,:)
	505	!$OMP THREADPRIVATE(reffcol)
	506
	507	! Non-oro GW tendencies
	508	REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
	509	REAL d_t_hin(ngrid,nlayer)
	510	! Diagnostics 2D of gw_nonoro
	511	REAL zustrhi(ngrid), zvstrhi(ngrid)
	512
	513
	514	real :: tsurf2(ngrid)
	515	!! real :: flux_o(ngrid),flux_g(ngrid)
	516	real :: flux_g(ngrid)
	517	real :: flux_sens_lat(ngrid)
	518	real :: qsurfint(ngrid,nq)
[3195]	519	#ifdef MESOSCALE
[3184]	520	REAL :: lsf_dt(nlayer)
	521	REAL :: lsf_dq(nlayer)
	522	#endif
	523
[3483]	524	! local variables for skin depth check
	525	real :: inertia_min,inertia_max
	526	real :: diurnal_skin ! diurnal skin depth (m)
	527	real :: annual_skin ! anuual skin depth (m)
	528
[3500]	529	! when no startfi file is asked for init
	530	real alpha,lay1 ! coefficients for building layers
	531	integer iloop
[3483]	532
[3184]	533	! flags to trigger extra sanity checks
	534	logical, save :: check_physics_inputs=.false.
	535	logical, save :: check_physics_outputs=.false.
[3195]	536	!$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)
[3184]	537
	538	! Misc
	539	character*2 :: str2
[3195]	540	character(len=10) :: tmp1
[3184]	541	character(len=10) :: tmp2
	542	!==================================================================================================
	543
	544	! -----------------
	545	! I. INITIALISATION
	546	! -----------------
	547
	548	! --------------------------------
	549	! I.1 First Call Initialisation.
	550	! --------------------------------
	551	if (firstcall) then
	552	call getin_p("check_physics_inputs", check_physics_inputs)
	553	call getin_p("check_physics_outputs", check_physics_outputs)
	554
	555	! Allocate saved arrays (except for 1D model, where this has already
	556	! been done)
	557	#ifndef MESOSCALE
	558	if (ngrid>1) call phys_state_var_init(nq)
	559	#endif
	560
	561	! Variables set to 0
	562	! ~~~~~~~~~~~~~~~~~~
	563	dtrad(:,:) = 0.0
	564	fluxrad(:) = 0.0
	565	tau_col(:) = 0.0
	566	zdtsw(:,:) = 0.0
	567	zdtlw(:,:) = 0.0
	568
	569	! Initialize tracer names, indexes and properties.
	570	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	571	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
	572	if (tracer) then
	573	call initracer(ngrid,nq)
	574	! if(photochem) then !AF24: removed
	575	endif
	576	! Initialize aerosol indexes.
	577	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3195]	578	! call iniaerosol
[3184]	579	! allocate related local arrays
	580	! (need be allocated instead of automatic because of "naerkind")
	581	allocate(aerosol(ngrid,nlayer,naerkind))
	582	allocate(reffcol(ngrid,naerkind))
	583
	584	#ifdef CPP_XIOS
	585	! Initialize XIOS context
	586	write(,) "physiq: call wxios_context_init"
	587	CALL wxios_context_init
	588	#endif
	589
	590	! Read 'startfi.nc' file.
	591	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3195]	592	#ifndef MESOSCALE
[3184]	593	call phyetat0(startphy_file, &
	594	ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
[3483]	595	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,inertiedat)
[3184]	596
[3195]	597	#else
	598
[3184]	599	day_ini = pday
	600	#endif
	601
	602	#ifndef MESOSCALE
	603	if (.not.startphy_file) then
[3500]	604	! starting without startfi.nc and with callsoil
	605	! is not yet possible as soildepth default is not defined
	606	if (callsoil) then
	607	! default mlayer distribution, following a power law:
	608	! mlayer(k)=lay1alpha*(k-1/2)
	609	lay1=2.e-4
	610	alpha=2
	611	do iloop=0,nsoilmx-1
	612	mlayer(iloop)=lay1(alpha*(iloop-0.5))
	613	enddo
	614	lay1=sqrt(mlayer(0)*mlayer(1))
	615	alpha=mlayer(1)/mlayer(0)
	616	do iloop=1,nsoilmx
	617	layer(iloop)=lay1(alpha*(iloop-1))
	618	enddo
	619	endif
[3184]	620	! additionnal "academic" initialization of physics
	621	if (is_master) write(,) "Physiq: initializing tsurf(:) to pt(:,1) !!"
	622	tsurf(:)=pt(:,1)
	623	if (is_master) write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
	624	do isoil=1,nsoilmx
	625	tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
	626	enddo
[3232]	627	if (is_master) write(,) "Physiq: initializing day_ini to pday !"
[3184]	628	day_ini=pday
	629	endif
	630	#endif
	631	if (pday.ne.day_ini) then
	632	write(,) "ERROR in physiq.F90:"
	633	write(,) "bad synchronization between physics and dynamics"
	634	write(,) "dynamics day: ",pday
	635	write(,) "physics day: ",day_ini
	636	stop
	637	endif
	638
	639	write (,) 'In physiq day_ini =', day_ini
	640
[3195]	641	! Initialize albedo calculation.
[3184]	642	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	643	albedo(:,:)=0.0
	644	albedo_bareground(:)=0.0
	645	albedo_snow_SPECTV(:)=0.0
	646	albedo_n2_ice_SPECTV(:)=0.0
[3412]	647
	648	ptime0=ptime
	649	write (,) 'In physiq ptime0 =', ptime
	650
[3195]	651	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV)
[3184]	652
[3195]	653	! Initialize orbital calculation.
[3184]	654	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	655	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
	656
[3329]	657	savedeclin=0.
	658	saveday=pday
[3377]	659	adjust=0. ! albedo adjustment for convergeps
[3184]	660
	661	! Initialize soil.
	662	! ~~~~~~~~~~~~~~~~
	663	if (callsoil) then
	664	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
	665	ptimestep,tsurf,tsoil,capcal,fluxgrd)
	666	else ! else of 'callsoil'.
	667	print*,'WARNING! Thermal conduction in the soil turned off'
	668	capcal(:)=1.e6
	669	fluxgrd(:)=intheat
	670	print*,'Flux from ground = ',intheat,' W m^-2'
	671	endif ! end of 'callsoil'.
	672
	673	icount=1
	674
	675	! Initialize surface history variable.
	676	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	677	qsurf_hist(:,:)=qsurf(:,:)
	678
[3275]	679	!! call WriteField_phy("post_qsurf_hist_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
[3184]	680
	681	! Initialize variable for dynamical heating and zonal wind tendency diagnostic
	682	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	683	ztprevious(:,:)=pt(:,:)
	684	zuprevious(:,:)=pu(:,:)
	685
	686	! Set temperature just above condensation temperature (for Early Mars)
	687	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	688	if(nearn2cond) then
	689	write(,)' WARNING! Starting at Tcond+1K'
	690	do l=1, nlayer
	691	do ig=1,ngrid
	692	pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4 &
	693	-pt(ig,l)) / ptimestep
	694	enddo
	695	enddo
	696	endif
	697
[3195]	698	if(meanOLR)then
	699	call system('rm -f rad_bal.out') ! to record global radiative balance.
	700	call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
[3184]	701	call system('rm -f h2o_bal.out') ! to record global hydrological balance.
	702	endif
	703
	704
	705	!! Initialize variables for water cycle ! AF24: removed
[3195]	706
[3184]	707	! Set metallicity for GCS
	708	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	709	metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
	710	call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic
[3195]	711
[3184]	712	! Set some parameters for the thermal plume model !AF24: removed
	713	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3195]	714
[3184]	715	#ifndef MESOSCALE
	716	if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
	717	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
	718	ptimestep,pday+nday,time_phys,cell_area, &
	719	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
	720	endif
	721
[3275]	722	!! call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
[3195]	723	#endif
	724	if (corrk) then
	725	! We initialise the spectral grid here instead of
	726	! at firstcall of callcorrk so we can output XspecIR, XspecVI
[3184]	727	! when using Dynamico
[3504]	728	if (is_master) print*, "physiq_mod: Correlated-k data base folder:",trim(datadir)
[3184]	729	call getin_p("corrkdir",corrkdir)
[3504]	730	if (is_master) print*,"corrkdir = ", corrkdir
[3184]	731	write (tmp1, '(i4)') L_NSPECTI
	732	write (tmp2, '(i4)') L_NSPECTV
	733	banddir=trim(trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2)))
	734	banddir=trim(trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir)))
	735	call setspi !Basic infrared properties.
	736	call setspv ! Basic visible properties.
	737	call sugas_corrk ! Set up gaseous absorption properties.
[3195]	738	if (aerohaze) then
	739	call suaer_corrk ! Set up aerosol optical properties.
	740	endif
[3184]	741	endif
	742
[3275]	743	!! call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
[3184]	744	! XIOS outputs
	745	#ifdef CPP_XIOS
	746
[3506]	747	if (is_master) write(,) "physiq: call initialize_xios_output"
	748	call initialize_xios_output(pday,ptime,ptimestep,daysec,year_day, &
	749	presnivs,pseudoalt,mlayer,WNOI,WNOV)
[3184]	750	#endif
[3195]	751
[3275]	752	!! call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
[3184]	753
	754	write(,) "physiq: end of firstcall"
	755	endif ! end of 'firstcall'
	756
[3275]	757	!! call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
	758	!! call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
[3184]	759
	760	if (check_physics_inputs) then
	761	!check the validity of input fields coming from the dynamics
	762	call check_physics_fields("begin physiq:", pt, pu, pv, pplev, pq)
	763	endif
	764
	765	! call writediagfi(ngrid,"pre_physical_rnat"," "," ",2,rnat)
	766	! call writediagfi(ngrid,"pre_physical_capcal"," "," ",2,capcal)
	767
	768	! ------------------------------------------------------
	769	! I.2 Initializations done at every physical timestep:
	770	! ------------------------------------------------------
	771
[3195]	772	#ifdef CPP_XIOS
[3184]	773	! update XIOS time/calendar
	774	call update_xios_timestep
[3195]	775	#endif
[3184]	776
	777	! Initialize various variables
	778	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3195]	779
[3184]	780	if ( .not.nearn2cond ) then
	781	pdt(1:ngrid,1:nlayer) = 0.0
[3195]	782	endif
[3184]	783	zdtsurf(1:ngrid) = 0.0
	784	pdq(1:ngrid,1:nlayer,1:nq) = 0.0
	785	dqsurf(1:ngrid,1:nq)= 0.0
	786	pdu(1:ngrid,1:nlayer) = 0.0
	787	pdv(1:ngrid,1:nlayer) = 0.0
	788	pdpsrf(1:ngrid) = 0.0
[3195]	789	zflubid(1:ngrid) = 0.0
[3184]	790	flux_sens_lat(1:ngrid) = 0.0
	791	taux(1:ngrid) = 0.0
	792	tauy(1:ngrid) = 0.0
	793
[3412]	794	if (conservn2) then
	795	write(,) 'conservn2 iniloop'
	796	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
	797	endif
	798
[3184]	799	zday=pday+ptime ! Compute time, in sols (and fraction thereof).
	800
	801	! Compute Stellar Longitude (Ls), and orbital parameters.
	802	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	803	if (season) then
	804	call stellarlong(zday,zls)
	805	else
	806	call stellarlong(noseason_day,zls)
	807	end if
	808
[3195]	809
[3237]	810	! Get Lyman alpha flux at specific Ls
	811	if (haze) then
	812	call lymalpha(zls,zfluxuv)
	813	print*, 'Haze lyman-alpha zls,zfluxuv=',zls,zfluxuv
	814	end if
	815
	816
	817	IF (triton) then
	818	CALL orbitetriton(zls,zday,dist_star,declin)
	819	ELSE
	820	call orbite(zls,dist_star,declin,right_ascen)
	821	ENDIF
	822
	823
[3184]	824	if (diurnal) then
	825	zlss=-2.pi(zday-.5)
	826	else if(diurnal .eqv. .false.) then
	827	zlss=9999.
[3195]	828	endif
[3184]	829
	830
	831	glat(:) = g !AF24: removed oblateness
	832
	833	! Compute geopotential between layers.
	834	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	835	zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
[3195]	836	do l=1,nlayer
[3184]	837	zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
	838	enddo
	839
	840	zzlev(1:ngrid,1)=0.
	841
	842	do l=2,nlayer
	843	do ig=1,ngrid
	844	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
	845	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
	846	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	847	enddo
[3195]	848	enddo
[3184]	849
	850	!Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
[3195]	851
[3184]	852	zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
	853
	854	! Compute potential temperature
	855	! Note : Potential temperature calculation may not be the same in physiq and dynamic...
	856	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3195]	857	do l=1,nlayer
[3184]	858	do ig=1,ngrid
	859	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
	860	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
	861	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
	862	massarea(ig,l)=mass(ig,l)*cell_area(ig)
	863	enddo
	864	enddo
	865
	866	! Compute vertical velocity (m/s) from vertical mass flux
	867	! w = F / (rhoarea) and rho = P/(rT)
	868	! But first linearly interpolate mass flux to mid-layers
[3237]	869	if (.not.fast) then
[3228]	870	do l=1,nlayer-1
[3184]	871	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
[3228]	872	enddo
	873	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
	874	do l=1,nlayer
[3184]	875	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) / &
	876	(pplay(1:ngrid,l)*cell_area(1:ngrid))
[3228]	877	enddo
	878	! omega in Pa/s
	879	do l=1,nlayer-1
[3184]	880	omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
[3228]	881	enddo
	882	omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
	883	do l=1,nlayer
[3184]	884	omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid)
[3228]	885	enddo
	886	endif
[3412]	887
	888	if (conservn2) then
	889	write(,) 'conservn2 thermo'
	890	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
	891	endif
	892
[3455]	893
	894
	895	! --------------------------------------------------------
	896	! 1.3 thermosphere
	897	! --------------------------------------------------------
	898
	899	! ajout de la conduction depuis la thermosphere
	900	IF (callconduct) THEN
	901
	902	call conduction (ngrid,nlayer,ptimestep, &
	903	pplay,pt,zzlay,zzlev,zdtconduc,tsurf)
	904	DO l=1,nlayer
	905	DO ig=1,ngrid
	906	pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
	907	ENDDO
	908	ENDDO
	909
	910	ENDIF
	911
	912	! ajout de la viscosite moleculaire
	913	IF (callmolvis) THEN
	914	call molvis(ngrid,nlayer,ptimestep, &
	915	pplay,pt,zzlay,zzlev, &
	916	zdtconduc,pu,tsurf,zdumolvis)
	917	call molvis(ngrid,nlayer,ptimestep, &
	918	pplay,pt,zzlay,zzlev, &
	919	zdtconduc,pv,tsurf,zdvmolvis)
	920
	921	DO l=1,nlayer
	922	DO ig=1,ngrid
	923	! pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
	924	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
	925	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
	926	ENDDO
	927	ENDDO
	928	ENDIF
	929
	930	IF (callmoldiff) THEN
	931	call moldiff_red(ngrid,nlayer,nq, &
	932	pplay,pplev,pt,pdt,pq,pdq,ptimestep, &
	933	zzlay,zdtconduc,zdqmoldiff)
	934
	935	DO l=1,nlayer
	936	DO ig=1,ngrid
	937	DO iq=1, nq
	938	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
	939	ENDDO
	940	ENDDO
	941	ENDDO
	942	ENDIF
	943
	944	if (conservn2) then
	945	write(,) 'conservn2 thermo'
	946	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
	947	endif
	948
	949
[3184]	950	!---------------------------------
	951	! II. Compute radiative tendencies
	952	!---------------------------------
[3243]	953	! Saving qsurf to compute paleo flux condensation/sublimation
	954	DO iq=1, nq
	955	DO ig=1,ngrid
	956	IF (qsurf(ig,iq).lt.0.) then
	957	qsurf(ig,iq)=0.
	958	ENDIF
	959	qsurf1(ig,iq)=qsurf(ig,iq)
	960	ENDDO
	961	ENDDO
[3184]	962
[3243]	963
[3184]	964	! Compute local stellar zenith angles
	965	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3377]	966	fract = 0
[3184]	967	if (diurnal) then
	968	ztim1=SIN(declin)
	969	ztim2=COS(declin)COS(2.pi*(zday-.5))
	970	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
	971
	972	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
	973	ztim1,ztim2,ztim3,mu0,fract)
	974	else if(diurnal .eqv. .false.) then
	975
	976	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad)
	977	! WARNING: this function appears not to work in 1D
	978
	979	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
	980	mu0 = cos(pi*szangle/180.0)
[3329]	981	fract= 1/(4*mu0) ! AF24: from pluto.old
[3184]	982	endif
	983
[3195]	984	endif
[3184]	985
[3237]	986
[3329]	987	! Pluto albedo /IT changes depending on surface ices (only in 3D)
	988	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	989	if (ngrid.ne.1) then
	990
	991	!! Specific to change albedo of N2 so that Psurf
	992	!! converges toward 1.4 Pa in "1989" seasons for Triton
	993	!! converges toward 1.1 Pa in "2015" seasons for Pluto
	994	if (convergeps) then
	995	if (triton) then
	996	! 1989 declination
	997	if (declin180./pi.gt.-46..and.declin180./pi.lt.-45. &
	998	.and.zday.gt.saveday+1000. &
	999	.and.declin.lt.savedeclin) then
	1000	call globalaverage2d(ngrid,pplev(:,1),globave)
	1001	if (globave.gt.1.5) then
	1002	adjust=adjust+0.005
	1003	else if (globave.lt.1.3) then
	1004	adjust=adjust-0.005
	1005	endif
	1006	saveday=zday
	1007	endif
	1008	else
	1009	! Pluto : 2015 declination current epoch
	1010	if (declin180./pi.gt.50.and.declin180./pi.lt.51. &
	1011	.and.zday.gt.saveday+10000. &
	1012	.and.declin.gt.savedeclin) then
	1013	call globalaverage2d(ngrid,pplev(:,1),globave)
	1014	if (globave.gt.1.2) then
	1015	adjust=adjust+0.005
	1016	else if (globave.lt.1.) then
	1017	adjust=adjust-0.005
	1018	endif
	1019	saveday=zday
	1020	endif
	1021	endif
	1022	end if
	1023	end if ! if ngrid ne 1
	1024
[3483]	1025	call surfprop(ngrid,nq,fract,qsurf,tsurf, &
	1026	capcal,adjust,dist_star,flusurfold,ptimestep,zls,&
	1027	albedo,emis,inertiedat)
[3329]	1028	! do i=2,ngrid
	1029	! albedo(i,:) = albedo(1,:)
	1030	! enddo
	1031
[3500]	1032	if (firstcall.and.callsoil) then
[3483]	1033	! AF24 Originally in soil.F, but inertiedat is modified by surfprop
	1034	! Additional checks: is the vertical discretization sufficient
	1035	! to resolve diurnal and annual waves?
	1036	do ig=1,ngrid
	1037	! extreme inertia for this column
	1038	inertia_min=minval(inertiedat(ig,:))
	1039	inertia_max=maxval(inertiedat(ig,:))
[3500]	1040	print*, "inertia min max" , inertia_min, inertia_max
[3483]	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) )
	2336	call wstats(ngrid,trim(noms(iq))//'_reff',trim(noms(iq))//'_reff',&
	2337	'm',3,reffrad(1,1,iq))
	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
	2359	call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
	2360	call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
	2361	call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
	2362	call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
	2363	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	2364	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
[3453]	2365	call writediagfi(ngrid,"emis","Emissivity","",2,emis)
[3184]	2366
[3241]	2367	!! Pluto outputs
	2368	call writediagfi(ngrid,"dist_star","dist_star","AU",0,dist_star)
	2369
	2370	if (fast) then
	2371	call writediagfi(ngrid,"globave","surf press","Pa",0,globave)
[3453]	2372	call writediagfi(ngrid,"fluxrad","fluxrad","W m-2",2,fluxrad)
	2373	call writediagfi(ngrid,"fluxgrd","fluxgrd","W m-2",2,fluxgrd)
	2374	call writediagfi(ngrid,"capcal","capcal","W.s m-2 K-1",2,capcal)
	2375	! call writediagfi(ngrid,"dplanck","dplanck","W.s m-2 K-1",2,dplanck)
	2376	call writediagfi(ngrid,"tsoil","tsoil","K",3,tsoil)
[3241]	2377	else
[3243]	2378	if (check_physics_outputs) then
	2379	! Check the validity of updated fields at the end of the physics step
	2380	call check_physics_fields("HERE physiq:", zt, zu, zv, pplev, zq)
	2381	endif
	2382
[3453]	2383	call writediagfi(ngrid,"temp","temperature","K",3,zt)
	2384	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
	2385	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
	2386	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
	2387	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
	2388	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
[3228]	2389	endif
	2390
[3184]	2391	! Subsurface temperatures
	2392	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	2393	! call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
	2394
	2395	! Total energy balance diagnostics
	2396	if(callrad)then
[3195]	2397
[3198]	2398	call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
[3361]	2399	call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
[3184]	2400	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
	2401	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
	2402	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
	2403	call writediagfi(ngrid,"shad","rings"," ", 2, fract)
[3195]	2404
[3184]	2405	! call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
	2406	! call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
	2407
	2408	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
	2409	! endif
[3195]	2410
[3184]	2411	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
[3195]	2412
[3184]	2413	endif ! end of 'callrad'
[3195]	2414
[3184]	2415	if(enertest) then
	2416	if (calldifv) then
[3195]	2417
[3184]	2418	call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
	2419	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
[3195]	2420
[3184]	2421	! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
	2422	! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
	2423	! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
[3195]	2424
[3184]	2425	endif
[3195]	2426
[3184]	2427	if (corrk) then
	2428	call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
	2429	call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
	2430	endif
[3195]	2431
[3184]	2432	if (generic_condensation) then
	2433	call writediagfi(ngrid,"genericconddE","heat from generic condensation","W m-2",2,genericconddE)
	2434	call writediagfi(ngrid,"dt_generic_condensation","heating from generic condensation","K s-1",3,dt_generic_condensation)
	2435	endif
[3195]	2436
[3184]	2437	endif ! end of 'enertest'
	2438
[3501]	2439	! Diagnostics of optical thickness
	2440	! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
	2441	if (diagdtau) then
	2442	do nw=1,L_NSPECTV
	2443	write(str2,'(i2.2)') nw
	2444	call writediagfi(ngrid,'dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',1,int_dtauv(:,nlayer:1:-1,nw))
	2445	enddo
	2446	do nw=1,L_NSPECTI
	2447	write(str2,'(i2.2)') nw
	2448	call writediagfi(ngrid,'dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',1,int_dtaui(:,nlayer:1:-1,nw))
	2449	enddo
	2450	endif
[3184]	2451
[3501]	2452	! Temporary inclusions for heating diagnostics.
	2453	call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
	2454	call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
	2455	call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
	2456	call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
[3195]	2457
[3501]	2458	! For Debugging.
	2459	!call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
	2460	!call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
[3184]	2461
	2462	! Output tracers.
	2463	if (tracer) then
[3453]	2464	! call writediagfi(ngrid,"zdtc","tendancy T cond N2","K",3,zdtc)
[3195]	2465
[3184]	2466	do iq=1,nq
[3453]	2467	call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
	2468	! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
[3422]	2469	! 'kg m^-2',2,qsurf_hist(1,iq) )
[3412]	2470	call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
	2471	'kg m^-2',2,qcol(1,iq) )
[3453]	2472	call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
[3422]	2473	'kg m^-2',2,qsurf(1,iq) )
[3184]	2474
	2475	if(generic_condensation)then
	2476	call writediagfi(ngrid,"rneb_generic","GCS cloud fraction (generic condensation)"," ",3,rneb_generic)
	2477	call writediagfi(ngrid,"CLF","GCS cloud fraction"," ",3,cloudfrac)
	2478	call writediagfi(ngrid,"RH_generic","GCS relative humidity"," ",3,RH_generic)
	2479	endif
[3390]	2480	! call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
[3184]	2481	enddo ! end of 'nq' loop
[3195]	2482
[3241]	2483	!Pluto specific
[3422]	2484	call writediagfi(ngrid,'n2_iceflux','n2_iceflux',"kg m^-2 s^-1",2,flusurf(1,igcm_n2) )
[3453]	2485	if (haze_radproffix)then
[3444]	2486	call writediagfi(ngrid,'haze_reff','haze_reff','m',3,reffrad(1,1,1))
	2487	end if
[3241]	2488	if (methane) then
[3453]	2489	! call writediagfi(ngrid,"rice_ch4","ch4 ice mass mean radius","m",3,rice_ch4)
	2490	! call writediagfi(ngrid,"zq1temp_ch4"," "," ",2,zq1temp_ch4)
	2491	! call writediagfi(ngrid,"qsat_ch4"," "," ",2,qsat_ch4)
	2492	! call writediagfi(ngrid,"qsat_ch4_l1"," "," ",2,qsat_ch4_l1)
	2493
[3241]	2494	call writediagfi(ngrid,'ch4_iceflux','ch4_iceflux',&
	2495	"kg m^-2 s^-1",2,flusurf(1,igcm_ch4_ice) )
	2496	call writediagfi(ngrid,"vmr_ch4","vmr_ch4","%",2,vmr_ch4)
	2497	if (.not.fast) then
	2498	call writediagfi(ngrid,"zrho_ch4","zrho_ch4","kg.m-3",3,zrho_ch4(:,:))
	2499	endif
[3243]	2500
[3390]	2501	! Tendancies
[3453]	2502	call writediagfi(ngrid,"zdqcn2_ch4","zdq condn2 ch4","",&
	2503	3,zdqc(:,:,igcm_ch4_gas))
	2504	call writediagfi(ngrid,"zdqdif_ch4","zdqdif ch4","",&
	2505	3,zdqdif(:,:,igcm_ch4_gas))
	2506	call writediagfi(ngrid,"zdqsdif_ch4_ice","zdqsdif ch4","",&
	2507	2,zdqsdif(:,igcm_ch4_ice))
	2508	call writediagfi(ngrid,"zdqadj_ch4","zdqadj ch4","",&
	2509	3,zdqadj(:,:,igcm_ch4_gas))
[3501]	2510	if (sedimentation) then
[3243]	2511	call writediagfi(ngrid,"zdqsed_ch4","zdqsed ch4","",&
	2512	3,zdqsed(:,:,igcm_ch4_gas))
	2513	call writediagfi(ngrid,"zdqssed_ch4","zdqssed ch4","",&
	2514	2,zdqssed(:,igcm_ch4_gas))
[3501]	2515	endif
	2516	if (metcloud.and.(.not.fast)) then
[3243]	2517	call writediagfi(ngrid,"zdtch4cloud","ch4 cloud","T s-1",&
[3241]	2518	3,zdtch4cloud)
[3439]	2519	call writediagfi(ngrid,"zdqch4cloud","ch4 cloud","T s-1",&
	2520	3,zdqch4cloud(1,1,igcm_ch4_gas))
	2521	endif
[3243]	2522
[3241]	2523	endif
[3243]	2524
[3241]	2525	if (carbox) then
[3453]	2526	! call writediagfi(ngrid,"zdtcocloud","tendancy T cocloud","K",3,zdtcocloud)
[3241]	2527	call writediagfi(ngrid,'co_iceflux','co_iceflux',&
	2528	"kg m^-2 s^-1",2,flusurf(1,igcm_co_ice) )
	2529	call writediagfi(ngrid,"vmr_co","vmr_co","%",2,vmr_co)
	2530	if (.not.fast) THEN
	2531	call writediagfi(ngrid,"zrho_co","zrho_co","kg.m-3",3,zrho_co(:,:))
	2532	endif
	2533	endif
[3243]	2534
[3241]	2535	if (haze) then
	2536	! call writediagfi(ngrid,"zrho_haze","zrho_haze","kg.m-3",3,zrho_haze(:,:))
	2537	call writediagfi(ngrid,"zdqrho_photprec","zdqrho_photprec",&
	2538	"kg.m-3.s-1",3,zdqrho_photprec(:,:))
	2539	call writediagfi(ngrid,"zdqphot_prec","zdqphot_prec","",&
	2540	3,zdqphot_prec(:,:))
	2541	call writediagfi(ngrid,"zdqhaze_ch4","zdqhaze_ch4","",&
	2542	3,zdqhaze(:,:,igcm_ch4_gas))
	2543	call writediagfi(ngrid,"zdqhaze_prec","zdqhaze_prec","",&
	2544	3,zdqhaze(:,:,igcm_prec_haze))
	2545	if (igcm_haze.ne.0) then
[3243]	2546	call writediagfi(ngrid,"zdqhaze_haze","zdqhaze_haze","",&
	2547	3,zdqhaze(:,:,igcm_haze))
	2548	if (sedimentation) then
	2549	call writediagfi(ngrid,"zdqssed_haze","zdqssed haze",&
	2550	"kg/m2/s",2,zdqssed(:,igcm_haze))
	2551	endif
[3241]	2552	endif
	2553	call writediagfi(ngrid,"zdqphot_ch4","zdqphot_ch4","",&
	2554	3,zdqphot_ch4(:,:))
	2555	call writediagfi(ngrid,"zdqconv_prec","zdqconv_prec","",&
	2556	3,zdqconv_prec(:,:))
	2557	! call writediagfi(ngrid,"zdqhaze_col","zdqhaze col","kg/m2/s",
	2558	! & 2,zdqhaze_col(:))
	2559	endif
[3243]	2560
[3241]	2561	if (aerohaze) then
	2562	call writediagfi(ngrid,"tau_col",&
	2563	"Total aerosol optical depth","opacity",2,tau_col)
	2564	endif
[3184]	2565
[3241]	2566	endif ! end of 'tracer'
[3184]	2567
[3241]	2568
[3184]	2569	! Output spectrum.
[3195]	2570	if(specOLR.and.corrk)then
[3184]	2571	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
	2572	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
	2573	call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W/m^2/cm^-1",3,GSR_nu)
	2574	endif
	2575
	2576	#else
	2577	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
	2578	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
	2579	comm_ALBEQ(1:ngrid)=albedo_equivalent(1:ngrid)
	2580	if (.not.calldifv) comm_LATENT_HF(:)=0.0
[3195]	2581
[3184]	2582	if ((tracer).and.(generic_condensation)) then
	2583
	2584	! If you have set generic_condensation (and not water) and you have set several GCS
	2585	! then the outputs given to WRF will be only the ones for the last generic tracer
	2586	! (Because it is rewritten every tracer in the loop)
	2587	! WRF can take only one moist tracer
	2588
	2589	do iq=1,nq
[3275]	2590	call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
[3195]	2591
[3275]	2592	if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
[3184]	2593
	2594	reffrad_generic_zeros_for_wrf(:,:) = 1.
	2595
	2596	comm_CLOUDFRAC(1:ngrid,1:nlayer) = cloudfrac(1:ngrid,1:nlayer)
	2597	comm_TOTCLOUDFRAC(1:ngrid) = totcloudfrac(1:ngrid) !??????
	2598	comm_SURFRAIN(1:ngrid) = zdqsrain_generic(1:ngrid,iq)
[3275]	2599	comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_gas)
[3184]	2600	comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_generic_ice)
	2601	! comm_H2OICE_REFF(1:ngrid,1:nlayer) = reffrad_generic_zeros_for_wrf(1:ngrid,1:nlayer) ! for now zeros !
	2602	!comm_H2OICE_REFF(1:ngrid,1:nlayer) = 0*zdtrain_generic(1:ngrid,1:nlayer) ! for now zeros !
	2603	comm_REEVAP(1:ngrid) = reevap_precip_generic(1:ngrid,iq)
	2604	comm_DTRAIN(1:ngrid,1:nlayer) = zdtrain_generic(1:ngrid,1:nlayer)
	2605	comm_DTLSC(1:ngrid,1:nlayer) = dt_generic_condensation(1:ngrid,1:nlayer)
	2606	comm_RH(1:ngrid,1:nlayer) = RH_generic(1:ngrid,1:nlayer,iq)
	2607
[3195]	2608	endif
	2609	end do ! do iq=1,nq loop on tracers
[3184]	2610
	2611	else
	2612	comm_CLOUDFRAC(1:ngrid,1:nlayer)=0.
	2613	comm_TOTCLOUDFRAC(1:ngrid)=0.
	2614	comm_SURFRAIN(1:ngrid)=0.
	2615	comm_DQVAP(1:ngrid,1:nlayer)=0.
	2616	comm_DQICE(1:ngrid,1:nlayer)=0.
	2617	! comm_H2OICE_REFF(1:ngrid,1:nlayer)=0.
	2618	comm_REEVAP(1:ngrid)=0.
	2619	comm_DTRAIN(1:ngrid,1:nlayer)=0.
	2620	comm_DTLSC(1:ngrid,1:nlayer)=0.
	2621	comm_RH(1:ngrid,1:nlayer)=0.
	2622
	2623	endif ! if water, if generic_condensation, else
	2624
	2625	comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid)
	2626	comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid)
	2627	comm_FLUXTOP_LW(1:ngrid)=fluxtop_lw(1:ngrid)
	2628	comm_FLUXSURF_SW(1:ngrid)=fluxsurf_sw(1:ngrid)
	2629	comm_FLUXSURF_LW(1:ngrid)=fluxsurf_lw(1:ngrid)
	2630	comm_FLXGRD(1:ngrid)=fluxgrd(1:ngrid)
	2631	sensibFlux(1:ngrid) = zflubid(1:ngrid) - capcal(1:ngrid)*zdtsdif(1:ngrid) !!! ????
	2632	comm_HR_DYN(1:ngrid,1:nlayer) = zdtdyn(1:ngrid,1:nlayer)
	2633	#endif
	2634
	2635	! XIOS outputs
[3195]	2636	#ifdef CPP_XIOS
[3184]	2637	! Send fields to XIOS: (NB these fields must also be defined as
	2638	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
[3506]	2639	CALL send_xios_field("controle",tab_cntrl_mod,1)
[3195]	2640
[3506]	2641	CALL send_xios_field("ap",ap,1)
	2642	CALL send_xios_field("bp",bp,1)
	2643	CALL send_xios_field("aps",aps,1)
	2644	CALL send_xios_field("bps",bps,1)
[3195]	2645
[3184]	2646	if (lastcall.and.is_omp_master) then
	2647	write(,) "physiq: call xios_context_finalize"
	2648	call xios_context_finalize
	2649	endif
	2650	#endif
[3195]	2651
[3184]	2652	if (check_physics_outputs) then
	2653	! Check the validity of updated fields at the end of the physics step
	2654	call check_physics_fields("end of physiq:", zt, zu, zv, pplev, zq)
	2655	endif
	2656
	2657	icount=icount+1
[3195]	2658
[3184]	2659	end subroutine physiq
[3195]	2660
[3184]	2661	end module physiq_mod

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

Download in other formats:

Original Format