Context Navigation

physiq_mod.F90 @ 1788

Last change on this file since 1788 was 1788, checked in by jvatant, 7 years ago
Get rid of all the old-generic dummy aerosol scheme ( just left scatterers_h for compilation ) as the new microphysics for Titan will be plugged in -> even removed sedimentation ( will be done in the microphysical model ) --JVO
Property svn:executable set to ``*
File size: 59.7 KB

Rev	Line
[1549]	1	module physiq_mod
	2
	3	implicit none
	4
	5	contains
	6
[751]	7	subroutine physiq(ngrid,nlayer,nq, &
[787]	8	nametrac, &
[253]	9	firstcall,lastcall, &
	10	pday,ptime,ptimestep, &
	11	pplev,pplay,pphi, &
	12	pu,pv,pt,pq, &
[1312]	13	flxw, &
[1576]	14	pdu,pdv,pdt,pdq,pdpsrf)
[253]	15
[1788]	16	use radinc_h, only : L_NSPECTI,L_NSPECTV
[1482]	17	use radcommon_h, only: sigma, glat, grav, BWNV
[1647]	18	use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe
[1216]	19	use comdiurn_h, only: coslat, sinlat, coslon, sinlon
[1327]	20	use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
[1216]	21	use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
[1543]	22	use geometry_mod, only: latitude, longitude, cell_area
[1542]	23	USE comgeomfi_h, only: totarea, totarea_planet
[1788]	24	USE tracer_h, only: noms, mmol
[1525]	25	use time_phylmdz_mod, only: ecritphy, iphysiq, nday
[1670]	26	use phyetat0_mod, only: phyetat0
[1216]	27	use phyredem, only: physdem0, physdem1
[1295]	28	use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
	29	use mod_phys_lmdz_para, only : is_master
[1308]	30	use planete_mod, only: apoastr, periastr, year_day, peri_day, &
	31	obliquit, nres, z0
[1524]	32	use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
	33	use time_phylmdz_mod, only: daysec
[1672]	34	use logic_mod, only: moyzon_ch
	35	use moyzon_mod, only: tmoy, playmoy, zphibar, zphisbar, zplevbar, &
	36	zplaybar, zzlevbar, zzlaybar, ztfibar
[1397]	37	use callkeys_mod
[1622]	38	use vertical_layers_mod, only: presnivs, pseudoalt
[1623]	39	#ifdef CPP_XIOS
[1622]	40	use xios_output_mod, only: initialize_xios_output, &
	41	update_xios_timestep, &
	42	send_xios_field
[1623]	43	#endif
[253]	44	implicit none
	45
	46
	47	!==================================================================
	48	!
	49	! Purpose
	50	! -------
	51	! Central subroutine for all the physics parameterisations in the
	52	! universal model. Originally adapted from the Mars LMDZ model.
	53	!
	54	! The model can be run without or with tracer transport
	55	! depending on the value of "tracer" in file "callphys.def".
	56	!
	57	!
	58	! It includes:
	59	!
[1477]	60	! I. Initialization :
	61	! I.1 Firstcall initializations.
	62	! I.2 Initialization for every call to physiq.
[253]	63	!
[1477]	64	! II. Compute radiative transfer tendencies (longwave and shortwave) :
	65	! II.a Option 1 : Call correlated-k radiative transfer scheme.
	66	! II.b Option 2 : Call Newtonian cooling scheme.
	67	! II.c Option 3 : Atmosphere has no radiative effect.
	68	!
	69	! III. Vertical diffusion (turbulent mixing) :
	70	!
	71	! IV. Dry Convective adjustment :
	72	!
[1647]	73	! V. Tracers
[1672]	74	! V.1. Chemistry
	75	! V.3. Updates (pressure variations, surface budget).
	76	! V.4. Surface Tracer Update.
[1477]	77	!
[1647]	78	! VI. Surface and sub-surface soil temperature.
[1477]	79	!
[1647]	80	! VII. Perform diagnostics and write output files.
[1477]	81	!
	82	!
[253]	83	! arguments
	84	! ---------
	85	!
[1477]	86	! INPUT
[253]	87	! -----
[1477]	88	!
[253]	89	! ngrid Size of the horizontal grid.
	90	! nlayer Number of vertical layers.
[1477]	91	! nq Number of advected fields.
	92	! nametrac Name of corresponding advected fields.
	93	!
	94	! firstcall True at the first call.
	95	! lastcall True at the last call.
	96	!
	97	! pday Number of days counted from the North. Spring equinoxe.
	98	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
	99	! ptimestep timestep (s).
	100	!
	101	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa).
	102	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
	103	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2.s-2).
	104	!
	105	! pu(ngrid,nlayer) u, zonal component of the wind (ms-1).
	106	! pv(ngrid,nlayer) v, meridional component of the wind (ms-1).
	107	!
	108	! pt(ngrid,nlayer) Temperature (K).
	109	!
	110	! pq(ngrid,nlayer,nq) Advected fields.
	111	!
[1216]	112	! pudyn(ngrid,nlayer) \
[253]	113	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
[1477]	114	! ptdyn(ngrid,nlayer) / corresponding variables.
[253]	115	! pqdyn(ngrid,nlayer,nq) /
[1312]	116	! flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
[253]	117	!
[1477]	118	! OUTPUT
[253]	119	! ------
	120	!
[1308]	121	! pdu(ngrid,nlayer) \
	122	! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding
	123	! pdt(ngrid,nlayer) / variables due to physical processes.
	124	! pdq(ngrid,nlayer) /
[253]	125	! pdpsrf(ngrid) /
	126	!
	127	!
	128	! Authors
	129	! -------
[1524]	130	! Frederic Hourdin 15/10/93
	131	! Francois Forget 1994
	132	! Christophe Hourdin 02/1997
[253]	133	! Subroutine completely rewritten by F. Forget (01/2000)
	134	! Water ice clouds: Franck Montmessin (update 06/2003)
	135	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
	136	! New correlated-k radiative scheme: R. Wordsworth (2009)
	137	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
	138	! Improved water cycle: R. Wordsworth / B. Charnay (2010)
	139	! To F90: R. Wordsworth (2010)
[594]	140	! New turbulent diffusion scheme: J. Leconte (2012)
[716]	141	! Loops converted to F90 matrix format: J. Leconte (2012)
[787]	142	! No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
[1477]	143	! Purge of the code : M. Turbet (2015)
[1672]	144	! Fork for Titan : J. Vatant d'Ollone (2017)
	145	! + clean of all too-generic (ocean, water, co2 ...) routines
	146	! + Titan's chemistry
[1647]	147	!============================================================================================
[253]	148
	149
	150	! 0. Declarations :
	151	! ------------------
	152
[1670]	153	include "netcdf.inc"
[253]	154
	155	! Arguments :
	156	! -----------
	157
[1477]	158	! INPUTS:
[253]	159	! -------
	160
[1477]	161	integer,intent(in) :: ngrid ! Number of atmospheric columns.
	162	integer,intent(in) :: nlayer ! Number of atmospheric layers.
	163	integer,intent(in) :: nq ! Number of tracers.
	164	character*20,intent(in) :: nametrac(nq) ! Names of the tracers taken from dynamics.
	165
	166	logical,intent(in) :: firstcall ! Signals first call to physics.
	167	logical,intent(in) :: lastcall ! Signals last call to physics.
	168
	169	real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0.
	170	real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
	171	real,intent(in) :: ptimestep ! Physics timestep (s).
	172	real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
	173	real,intent(in) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
	174	real,intent(in) :: pphi(ngrid,nlayer) ! Geopotential at mid-layer (m2s-2).
	175	real,intent(in) :: pu(ngrid,nlayer) ! Zonal wind component (m/s).
	176	real,intent(in) :: pv(ngrid,nlayer) ! Meridional wind component (m/s).
	177	real,intent(in) :: pt(ngrid,nlayer) ! Temperature (K).
	178	real,intent(in) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
	179	real,intent(in) :: flxw(ngrid,nlayer) ! Vertical mass flux (ks/s) at lower boundary of layer
[253]	180
[1477]	181	! OUTPUTS:
[253]	182	! --------
	183
[1477]	184	! Physical tendencies :
	185
	186	real,intent(out) :: pdu(ngrid,nlayer) ! Zonal wind tendencies (m/s/s).
	187	real,intent(out) :: pdv(ngrid,nlayer) ! Meridional wind tendencies (m/s/s).
	188	real,intent(out) :: pdt(ngrid,nlayer) ! Temperature tendencies (K/s).
	189	real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
	190	real,intent(out) :: pdpsrf(ngrid) ! Surface pressure tendency (Pa/s).
	191
[253]	192	! Local saved variables:
	193	! ----------------------
	194
[1622]	195	integer,save :: day_ini ! Initial date of the run (sol since Ls=0).
	196	integer,save :: icount ! Counter of calls to physiq during the run.
	197	!$OMP THREADPRIVATE(day_ini,icount)
	198
[1482]	199	real, dimension(:),allocatable,save :: tsurf ! Surface temperature (K).
	200	real, dimension(:,:),allocatable,save :: tsoil ! Sub-surface temperatures (K).
	201	real, dimension(:,:),allocatable,save :: albedo ! Surface Spectral albedo. By MT2015.
	202	real, dimension(:),allocatable,save :: albedo_equivalent ! Spectral Mean albedo.
[1477]	203
[1647]	204	!$OMP THREADPRIVATE(tsurf,tsoil,albedo,albedo_equivalent)
[253]	205
[1482]	206	real,dimension(:),allocatable,save :: albedo_bareground ! Bare Ground Albedo. By MT 2015.
[1477]	207
[1647]	208	!$OMP THREADPRIVATE(albedo_bareground)
[253]	209
[1477]	210	real,dimension(:),allocatable,save :: emis ! Thermal IR surface emissivity.
	211	real,dimension(:,:),allocatable,save :: dtrad ! Net atmospheric radiative heating rate (K.s-1).
	212	real,dimension(:),allocatable,save :: fluxrad_sky ! Radiative flux from sky absorbed by surface (W.m-2).
	213	real,dimension(:),allocatable,save :: fluxrad ! Net radiative surface flux (W.m-2).
	214	real,dimension(:),allocatable,save :: capcal ! Surface heat capacity (J m-2 K-1).
	215	real,dimension(:),allocatable,save :: fluxgrd ! Surface conduction flux (W.m-2).
	216	real,dimension(:,:),allocatable,save :: qsurf ! Tracer on surface (e.g. kg.m-2).
	217	real,dimension(:,:),allocatable,save :: q2 ! Turbulent Kinetic Energy.
	218
[1315]	219	!$OMP THREADPRIVATE(emis,dtrad,fluxrad_sky,fluxrad,capcal,fluxgrd,qsurf,q2)
[253]	220
	221
	222	! Local variables :
	223	! -----------------
	224
[1477]	225	real zh(ngrid,nlayer) ! Potential temperature (K).
	226	real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
	227
[1670]	228	integer l,ig,ierr,iq,nw,isoil
[1161]	229
[1477]	230	! FOR DIAGNOSTIC :
	231
[1482]	232	real,dimension(:),allocatable,save :: fluxsurf_lw ! Incident Long Wave (IR) surface flux (W.m-2).
	233	real,dimension(:),allocatable,save :: fluxsurf_sw ! Incident Short Wave (stellar) surface flux (W.m-2).
	234	real,dimension(:),allocatable,save :: fluxsurfabs_sw ! Absorbed Short Wave (stellar) flux by the surface (W.m-2).
	235	real,dimension(:),allocatable,save :: fluxtop_lw ! Outgoing LW (IR) flux to space (W.m-2).
	236	real,dimension(:),allocatable,save :: fluxabs_sw ! Absorbed SW (stellar) flux (W.m-2).
	237	real,dimension(:),allocatable,save :: fluxtop_dn ! Incoming SW (stellar) radiation at the top of the atmosphere (W.m-2).
	238	real,dimension(:),allocatable,save :: fluxdyn ! Horizontal heat transport by dynamics (W.m-2).
	239	real,dimension(:,:),allocatable,save :: OLR_nu ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1)).
	240	real,dimension(:,:),allocatable,save :: OSR_nu ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1)).
	241	real,dimension(:,:),allocatable,save :: zdtlw ! LW heating tendencies (K/s).
	242	real,dimension(:,:),allocatable,save :: zdtsw ! SW heating tendencies (K/s).
	243	real,dimension(:),allocatable,save :: sensibFlux ! Turbulent flux given by the atmosphere to the surface (W.m-2).
[1477]	244
[1482]	245	!$OMP THREADPRIVATE(fluxsurf_lw,fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,&
[1524]	246
[1477]	247	!$OMP zdtlw,zdtsw,sensibFlux)
[253]	248
[1477]	249	real zls ! Solar longitude (radians).
	250	real zlss ! Sub solar point longitude (radians).
	251	real zday ! Date (time since Ls=0, calculated in sols).
	252	real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers.
	253	real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries.
[253]	254
[1477]	255	! TENDENCIES due to various processes :
[253]	256
[1477]	257	! For Surface Temperature : (K/s)
	258	real zdtsurf(ngrid) ! Cumulated tendencies.
	259	real zdtsurfmr(ngrid) ! Mass_redistribution routine.
	260	real zdtsdif(ngrid) ! Turbdiff/vdifc routines.
	261
	262	! For Atmospheric Temperatures : (K/s)
[1647]	263	real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
[1477]	264	real zdtmr(ngrid,nlayer) ! Mass_redistribution routine.
	265	real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
	266
	267	! For Surface Tracers : (kg/m2/s)
	268	real dqsurf(ngrid,nq) ! Cumulated tendencies.
	269	real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines.
	270	real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine.
	271
	272	! For Tracers : (kg/kg_of_air/s)
	273	real zdqadj(ngrid,nlayer,nq) ! Convadj routine.
	274	real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines.
	275	real zdqevap(ngrid,nlayer) ! Turbdiff routine.
	276	real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine.
[1672]	277
	278	real zdqchi(ngrid,nlayer,nq) ! Chemical tendency ( chemistry routine ).
	279	real zdqmph(ngrid,nlayer,nq) ! Microphysical tendency ( condensation routine only for now, no microphysical routines ).
[1477]	280
	281	! For Winds : (m/s/s)
	282	real zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! Convadj routine.
	283	real zdumr(ngrid,nlayer),zdvmr(ngrid,nlayer) ! Mass_redistribution routine.
	284	real zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	285	real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	286	real zdhadj(ngrid,nlayer) ! Convadj routine.
	287
	288	! For Pressure and Mass :
	289	real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
	290	real zdmassmr_col(ngrid) ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
	291	real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s).
[253]	292
[1477]	293
	294
	295	! Local variables for LOCAL CALCULATIONS:
	296	! ---------------------------------------
[787]	297	real zflubid(ngrid)
[1308]	298	real zplanck(ngrid),zpopsk(ngrid,nlayer)
[253]	299	real ztim1,ztim2,ztim3, z1,z2
	300	real ztime_fin
[1308]	301	real zdh(ngrid,nlayer)
[1194]	302	real gmplanet
[1297]	303	real taux(ngrid),tauy(ngrid)
[1194]	304
[253]	305
[1477]	306	! local variables for DIAGNOSTICS : (diagfi & stat)
	307	! -------------------------------------------------
	308	real ps(ngrid) ! Surface Pressure.
	309	real zt(ngrid,nlayer) ! Atmospheric Temperature.
	310	real zu(ngrid,nlayer),zv(ngrid,nlayer) ! Zonal and Meridional Winds.
	311	real zq(ngrid,nlayer,nq) ! Atmospheric Tracers.
	312	real zdtadj(ngrid,nlayer) ! Convadj Diagnostic.
	313	real zdtdyn(ngrid,nlayer) ! Dynamical Heating (K/s).
[1637]	314	real zdudyn(ngrid,nlayer) ! Dynamical Zonal Wind tendency (m.s-2).
[1477]	315	real,allocatable,dimension(:,:),save :: ztprevious ! Previous loop Atmospheric Temperature (K) ! Useful for Dynamical Heating calculation.
[1637]	316	real,allocatable,dimension(:,:),save :: zuprevious ! Previous loop Zonal Wind (m.s-1) ! Useful for Zonal Wind tendency calculation.
	317	!$OMP THREADPRIVATE(ztprevious,zuprevious)
[253]	318
[1477]	319	real vmr(ngrid,nlayer) ! volume mixing ratio
[253]	320	real time_phys
[597]	321
[1477]	322	real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
	323
[594]	324	! to test energy conservation (RW+JL)
[1308]	325	real mass(ngrid,nlayer),massarea(ngrid,nlayer)
[651]	326	real dEtot, dEtots, AtmToSurf_TurbFlux
[959]	327	real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
[1315]	328	!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
[1308]	329	real dEzRadsw(ngrid,nlayer),dEzRadlw(ngrid,nlayer),dEzdiff(ngrid,nlayer)
[787]	330	real dEdiffs(ngrid),dEdiff(ngrid)
[1477]	331
[594]	332	!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
[1477]	333
[1295]	334	real dItot, dItot_tmp, dVtot, dVtot_tmp
[1647]	335
[1295]	336	real dWtot, dWtot_tmp, dWtots, dWtots_tmp
[1477]	337
	338
[1482]	339	! For Clear Sky Case.
	340	real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxsurfabs_sw1(ngrid) ! For SW/LW flux.
	341	real fluxtop_lw1(ngrid), fluxabs_sw1(ngrid) ! For SW/LW flux.
	342	real albedo_equivalent1(ngrid) ! For Equivalent albedo calculation.
[1647]	343	real tf, ntf
[253]	344
[787]	345	real,allocatable,dimension(:,:),save :: qsurf_hist
[1315]	346	!$OMP THREADPRIVATE(qsurf_hist)
[253]	347
[1672]	348	! Local variables for Titan chemistry and microphysics (JVO 2017)
	349	! ----------------------------------------------------------------------------
	350
	351	integer,parameter :: nmicro=0 ! Temporary ! To be put in start/def
[996]	352
[1672]	353	real ctimestep ! Chemistry timestep (s)
	354
	355	! Grandeurs en moyennes zonales ------------------------
	356	real temp_eq(nlayer), press_eq(nlayer)
	357	real zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
	358	! real zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)
	359	real ztemp(ngrid,nlayer)
	360
	361	real ychim(ngrid,nlayer,nq-nmicro)
	362
[1787]	363	real rat_mmol(nq) ! Molar fraction ratio
	364
[1672]	365	! 2D vmr tendencies ( chemistry and condensation )
	366	real,dimension(:,:,:),allocatable,save :: dycchi
	367	! Must be saved since chemistry is not called every step
	368
	369	real dycmph(ngrid,nlayer,nq-nmicro)
	370	! ----------------------------------------------------------
	371
	372	real,dimension(:,:),allocatable,save :: qysat
	373
	374	character*10,dimension(:),allocatable,save :: nomqy
	375	!$OMP THREADPRIVATE(dycchi,qysat,nomqy)
	376
	377	!-----------------------------------------------------------------------------
	378
[1477]	379	!==================================================================================================
[253]	380
	381	! -----------------
[1477]	382	! I. INITIALISATION
	383	! -----------------
[253]	384
[1477]	385	! --------------------------------
	386	! I.1 First Call Initialisation.
	387	! --------------------------------
[253]	388	if (firstcall) then
	389
[1477]	390	! Allocate saved arrays.
[858]	391	ALLOCATE(tsurf(ngrid))
	392	ALLOCATE(tsoil(ngrid,nsoilmx))
[1482]	393	ALLOCATE(albedo(ngrid,L_NSPECTV))
[1647]	394	ALLOCATE(albedo_equivalent(ngrid))
	395	ALLOCATE(albedo_bareground(ngrid))
[858]	396	ALLOCATE(emis(ngrid))
[1308]	397	ALLOCATE(dtrad(ngrid,nlayer))
[858]	398	ALLOCATE(fluxrad_sky(ngrid))
	399	ALLOCATE(fluxrad(ngrid))
	400	ALLOCATE(capcal(ngrid))
	401	ALLOCATE(fluxgrd(ngrid))
	402	ALLOCATE(qsurf(ngrid,nq))
[1308]	403	ALLOCATE(q2(ngrid,nlayer+1))
	404	ALLOCATE(ztprevious(ngrid,nlayer))
[1637]	405	ALLOCATE(zuprevious(ngrid,nlayer))
[858]	406	ALLOCATE(qsurf_hist(ngrid,nq))
	407	ALLOCATE(fluxsurf_lw(ngrid))
	408	ALLOCATE(fluxsurf_sw(ngrid))
[1482]	409	ALLOCATE(fluxsurfabs_sw(ngrid))
[858]	410	ALLOCATE(fluxtop_lw(ngrid))
	411	ALLOCATE(fluxabs_sw(ngrid))
	412	ALLOCATE(fluxtop_dn(ngrid))
	413	ALLOCATE(fluxdyn(ngrid))
	414	ALLOCATE(OLR_nu(ngrid,L_NSPECTI))
	415	ALLOCATE(OSR_nu(ngrid,L_NSPECTV))
	416	ALLOCATE(sensibFlux(ngrid))
[1308]	417	ALLOCATE(zdtlw(ngrid,nlayer))
	418	ALLOCATE(zdtsw(ngrid,nlayer))
[1672]	419	ALLOCATE(dycchi(ngrid,nlayer,nq-nmicro))
	420	ALLOCATE(qysat(nlayer,nq-nmicro))
	421	ALLOCATE(nomqy(nq-nmicro+1))
	422
[1477]	423	! This is defined in comsaison_h
[1161]	424	ALLOCATE(mu0(ngrid))
[1477]	425	ALLOCATE(fract(ngrid))
	426	! This is defined in radcommon_h
	427	ALLOCATE(glat(ngrid))
	428
[787]	429
[1477]	430	! Variables set to 0
[253]	431	! ~~~~~~~~~~~~~~~~~~
	432	dtrad(:,:) = 0.0
	433	fluxrad(:) = 0.0
	434	zdtsw(:,:) = 0.0
	435	zdtlw(:,:) = 0.0
[726]	436
[253]	437
[1477]	438	! Initialize tracer names, indexes and properties.
	439	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	440	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
[253]	441	if (tracer) then
[787]	442	call initracer(ngrid,nq,nametrac)
[1477]	443	endif
[253]	444
[1787]	445	rat_mmol(:) = mmol(:) / mugaz
[726]	446
[1477]	447	! Read 'startfi.nc' file.
[253]	448	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1672]	449	call phyetat0(startphy_file,ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
	450	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf)
[1670]	451	if (.not.startphy_file) then
	452	! additionnal "academic" initialization of physics
	453	if (is_master) write(,) "Physiq: initializing tsurf(:) to pt(:,1) !!"
	454	tsurf(:)=pt(:,1)
	455	if (is_master) write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
	456	do isoil=1,nsoilmx
	457	tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
	458	enddo
	459	if (is_master) write(,) "Physiq: initializing day_ini to pdat !"
	460	day_ini=pday
	461	endif
[253]	462
	463	if (pday.ne.day_ini) then
	464	write(,) "ERROR in physiq.F90:"
	465	write(,) "bad synchronization between physics and dynamics"
	466	write(,) "dynamics day: ",pday
	467	write(,) "physics day: ",day_ini
	468	stop
	469	endif
	470
	471	write (,) 'In physiq day_ini =', day_ini
	472
[1482]	473	! Initialize albedo calculation.
	474	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	475	albedo(:,:)=0.0
[1524]	476	albedo_bareground(:)=0.0
[1647]	477	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground)
[1482]	478
	479	! Initialize orbital calculation.
	480	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[253]	481	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
	482
	483
	484	if(tlocked)then
	485	print*,'Planet is tidally locked at resonance n=',nres
	486	print*,'Make sure you have the right rotation rate!!!'
	487	endif
	488
[1297]	489
[1477]	490	! Initialize soil.
	491	! ~~~~~~~~~~~~~~~~
[253]	492	if (callsoil) then
[1477]	493
[787]	494	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
[1477]	495	ptimestep,tsurf,tsoil,capcal,fluxgrd)
[1297]	496
[1477]	497	else ! else of 'callsoil'.
	498
[253]	499	print*,'WARNING! Thermal conduction in the soil turned off'
[918]	500	capcal(:)=1.e6
[952]	501	fluxgrd(:)=intheat
	502	print*,'Flux from ground = ',intheat,' W m^-2'
[1477]	503
	504	endif ! end of 'callsoil'.
	505
[253]	506	icount=1
[1477]	507
[253]	508
[1477]	509	! Initialize surface history variable.
[253]	510	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[728]	511	qsurf_hist(:,:)=qsurf(:,:)
[253]	512
[1637]	513	! Initialize variable for dynamical heating and zonal wind tendency diagnostic
	514	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[253]	515	ztprevious(:,:)=pt(:,:)
[1637]	516	zuprevious(:,:)=pu(:,:)
[253]	517
	518
[1477]	519	if(meanOLR)then
	520	call system('rm -f rad_bal.out') ! to record global radiative balance.
	521	call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
	522	call system('rm -f h2o_bal.out') ! to record global hydrological balance.
[253]	523	endif
	524
[1477]	525	if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
[1542]	526	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
	527	ptimestep,pday+nday,time_phys,cell_area, &
[1482]	528	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
[1216]	529	endif
[1672]	530
	531	! Sanity check for chemistry
	532	if ( ((moyzon_ch).and.(.not.callchim)) .or. ((.not.moyzon_ch).and.(callchim)) ) then
	533	print *, "moyzon_ch=",moyzon_ch," and callchim=",callchim
	534	print *, "This is not compatible..."
	535	stop
	536	endif
	537
	538	! Initialize names, timestep and saturation profiles for chemistry
	539
	540	if ( callchim .and. (nq.gt.nmicro) ) then
	541
	542	ctimestep = ptimestep*REAL(ichim)
	543
	544	do iq=nmicro+1,nq
	545	nomqy(iq-nmicro) = nametrac(iq)
	546	enddo
	547
	548	nomqy(nq-nmicro+1) = "HV"
	549
	550	! qysat is taken at the equator ( small variations of t,p)
	551	temp_eq(:) = tmoy(:)
	552	press_eq(:) = playmoy(:)/100. ! en mbar
	553
	554	call inicondens(nq-nmicro,press_eq,temp_eq,nomqy,qysat)
	555
	556	endif
[1216]	557
[1622]	558	! XIOS outputs
	559	#ifdef CPP_XIOS
	560
	561	write(,) "physiq: call initialize_xios_output"
	562	call initialize_xios_output(pday,ptime,ptimestep,daysec, &
	563	presnivs,pseudoalt)
	564	#endif
[1477]	565	endif ! end of 'firstcall'
[253]	566
[1477]	567	! ------------------------------------------------------
	568	! I.2 Initializations done at every physical timestep:
	569	! ------------------------------------------------------
	570
[1622]	571	#ifdef CPP_XIOS
	572	! update XIOS time/calendar
	573	call update_xios_timestep
	574	#endif
	575
[1477]	576	! Initialize various variables
	577	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1647]	578
	579	pdt(1:ngrid,1:nlayer) = 0.0
[1477]	580	zdtsurf(1:ngrid) = 0.0
[1308]	581	pdq(1:ngrid,1:nlayer,1:nq) = 0.0
[1477]	582	dqsurf(1:ngrid,1:nq)= 0.0
	583	pdu(1:ngrid,1:nlayer) = 0.0
	584	pdv(1:ngrid,1:nlayer) = 0.0
[787]	585	pdpsrf(1:ngrid) = 0.0
[1477]	586	zflubid(1:ngrid) = 0.0
[1297]	587	taux(1:ngrid) = 0.0
	588	tauy(1:ngrid) = 0.0
[253]	589
[1477]	590	zday=pday+ptime ! Compute time, in sols (and fraction thereof).
[1297]	591
[1477]	592	! Compute Stellar Longitude (Ls), and orbital parameters.
	593	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[253]	594	if (season) then
	595	call stellarlong(zday,zls)
	596	else
	597	call stellarlong(float(day_ini),zls)
	598	end if
	599
[1329]	600	call orbite(zls,dist_star,declin,right_ascen)
[1477]	601
[1329]	602	if (tlocked) then
	603	zlss=Mod(-(2.pi(zday/year_day)nres - right_ascen),2.pi)
	604	elseif (diurnal) then
[1524]	605	zlss=-2.pi(zday-.5)
[1329]	606	else if(diurnal .eqv. .false.) then
	607	zlss=9999.
	608	endif
[1194]	609
	610
[1477]	611	! Compute variations of g with latitude (oblate case).
	612	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	613	if (oblate .eqv. .false.) then
	614	glat(:) = g
	615	else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
	616	print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
	617	call abort
	618	else
	619	gmplanet = MassPlanetgrav1e24
	620	do ig=1,ngrid
[1542]	621	glat(ig)= gmplanet/(Rmean*2) (1.D0 + 0.75 J2 - 2.0flatten/3. + (2.flatten - 15./4. J2) * cos(2. * (pi/2. - latitude(ig))))
[1477]	622	end do
	623	endif
[1297]	624
[1477]	625	! Compute geopotential between layers.
	626	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1308]	627	zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
	628	do l=1,nlayer
[1477]	629	zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
[1194]	630	enddo
[728]	631
[787]	632	zzlev(1:ngrid,1)=0.
[1477]	633	zzlev(1:ngrid,nlayer+1)=1.e7 ! Dummy top of last layer above 10000 km...
[728]	634
[253]	635	do l=2,nlayer
	636	do ig=1,ngrid
	637	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
	638	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
	639	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	640	enddo
[1477]	641	enddo
[253]	642
[1672]	643	! -------------------------------Taken from old Titan --------------------------
	644	! zonal averages needed
	645	if (moyzon_ch) then
	646
	647	zzlaybar(1,:)=(zphibar(1,:)+zphisbar(1))/g
	648	! SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
	649	! zzlaybar(1,:)=RGRARA/(RG*RA-(zphibar(1,:)+zphisbar(1)))-RA
	650	zzlevbar(1,1)=zphisbar(1)/g
	651	DO l=2,nlayer
	652	z1=(zplaybar(1,l-1)+zplevbar(1,l))/(zplevbar(1,l-1)-zplevbar(1,l))
	653	z2=(zplevbar(1,l) +zplaybar(1,l))/(zplevbar(1,l) -zplaybar(1,l))
	654	zzlevbar(1,l)=(z1zzlaybar(1,l-1)+z2zzlaybar(1,l))/(z1+z2)
	655	ENDDO
	656	zzlevbar(1,nlayer+1)=zzlaybar(1,nlayer)+(zzlaybar(1,nlayer)-zzlevbar(1,nlayer))
	657
	658	DO ig=2,ngrid
	659	if (latitude(ig).ne.latitude(ig-1)) then
	660	DO l=1,nlayer
	661	zzlaybar(ig,l)=(zphibar(ig,l)+zphisbar(ig))/g
	662	! SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
	663	!zzlaybar(ig,l)=RGRARA/(RG*RA-(zphibar(ig,l)+zphisbar(ig)))-RA
	664	ENDDO
	665	zzlevbar(ig,1)=zphisbar(ig)/g
	666	DO l=2,nlayer
	667	z1=(zplaybar(ig,l-1)+zplevbar(ig,l))/ (zplevbar(ig,l-1)-zplevbar(ig,l))
	668	z2=(zplevbar(ig,l) +zplaybar(ig,l))/(zplevbar(ig,l) -zplaybar(ig,l))
	669	zzlevbar(ig,l)=(z1zzlaybar(ig,l-1)+z2zzlaybar(ig,l))/(z1+z2)
	670	ENDDO
	671	zzlevbar(ig,nlayer+1)=zzlaybar(ig,nlayer)+(zzlaybar(ig,nlayer)-zzlevbar(ig,nlayer))
	672	else
	673	zzlaybar(ig,:)=zzlaybar(ig-1,:)
	674	zzlevbar(ig,:)=zzlevbar(ig-1,:)
	675	endif
	676	ENDDO
	677
	678	endif ! moyzon
	679	! -------------------------------------------------------------------------------------
	680
[1477]	681	! Compute potential temperature
	682	! Note : Potential temperature calculation may not be the same in physiq and dynamic...
	683	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[597]	684	do l=1,nlayer
[787]	685	do ig=1,ngrid
[253]	686	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
[597]	687	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
[1194]	688	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
[1542]	689	massarea(ig,l)=mass(ig,l)*cell_area(ig)
[253]	690	enddo
	691	enddo
	692
[1312]	693	! Compute vertical velocity (m/s) from vertical mass flux
[1346]	694	! w = F / (rhoarea) and rho = P/(rT)
[1477]	695	! But first linearly interpolate mass flux to mid-layers
	696	do l=1,nlayer-1
	697	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
	698	enddo
	699	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
	700	do l=1,nlayer
	701	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) / &
[1542]	702	(pplay(1:ngrid,l)*cell_area(1:ngrid))
[1477]	703	enddo
[1194]	704
[1477]	705	!---------------------------------
	706	! II. Compute radiative tendencies
	707	!---------------------------------
[253]	708
	709	if (callrad) then
[526]	710	if( mod(icount-1,iradia).eq.0.or.lastcall) then
[253]	711
[1477]	712	! Compute local stellar zenith angles
	713	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	714	if (tlocked) then
	715	! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb
	716	ztim1=SIN(declin)
	717	ztim2=COS(declin)*COS(zlss)
	718	ztim3=COS(declin)*SIN(zlss)
[253]	719
[1477]	720	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
	721	ztim1,ztim2,ztim3,mu0,fract, flatten)
[253]	722
[1477]	723	elseif (diurnal) then
	724	ztim1=SIN(declin)
	725	ztim2=COS(declin)COS(2.pi*(zday-.5))
	726	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
[253]	727
[1477]	728	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
	729	ztim1,ztim2,ztim3,mu0,fract, flatten)
	730	else if(diurnal .eqv. .false.) then
[253]	731
[1542]	732	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad,flatten)
[1161]	733	! WARNING: this function appears not to work in 1D
[253]	734
[1477]	735	endif
[1161]	736
[1477]	737	! Eclipse incoming sunlight (e.g. Saturn ring shadowing).
[1429]	738	if(rings_shadow) then
	739	call call_rings(ngrid, ptime, pday, diurnal)
	740	endif
[1133]	741
[1329]	742
[1477]	743	if (corrk) then
	744
	745	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	746	! II.a Call correlated-k radiative transfer scheme
	747	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1297]	748
[1648]	749	call call_profilgases(nlayer)
	750
[1477]	751	! standard callcorrk
[1482]	752	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
	753	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
[1788]	754	tsurf,fract,dist_star, &
[1482]	755	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, &
	756	fluxsurfabs_sw,fluxtop_lw, &
	757	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, &
[1788]	758	firstcall,lastcall)
[1297]	759
[1482]	760	! Radiative flux from the sky absorbed by the surface (W.m-2).
[1477]	761	GSR=0.0
[1482]	762	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurfabs_sw(1:ngrid)
[253]	763
[1477]	764	!if(noradsurf)then ! no lower surface; SW flux just disappears
[1542]	765	! GSR = SUM(fluxsurf_sw(1:ngrid)*cell_area(1:ngrid))/totarea
[1477]	766	! fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
	767	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
	768	!endif
[253]	769
[1477]	770	! Net atmospheric radiative heating rate (K.s-1)
	771	dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
[1498]	772
[1477]	773	elseif(newtonian)then
[1482]	774
	775	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	776	! II.b Call Newtonian cooling scheme
	777	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1477]	778	call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
[253]	779
[1477]	780	zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep
	781	! e.g. surface becomes proxy for 1st atmospheric layer ?
[253]	782
[1477]	783	else
	784
	785	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	786	! II.c Atmosphere has no radiative effect
	787	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	788	fluxtop_dn(1:ngrid) = fract(1:ngrid)mu0(1:ngrid)Fat1AU/dist_star**2
	789	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
	790	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
	791	endif
	792	fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
[1482]	793	print*,'------------WARNING---WARNING------------' ! by MT2015.
	794	print*,'You are in corrk=false mode, '
[1498]	795	print*,'and the surface albedo is taken equal to the first visible spectral value'
	796
	797	fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
	798	fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid)
[1477]	799	fluxtop_lw(1:ngrid) = emis(1:ngrid)sigmatsurf(1:ngrid)**4
[253]	800
[1477]	801	dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating
[253]	802
[1477]	803	endif ! end of corrk
[253]	804
[1477]	805	endif ! of if(mod(icount-1,iradia).eq.0)
[787]	806
[253]	807
[1477]	808	! Transformation of the radiative tendencies
	809	! ------------------------------------------
	810	zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
	811	zplanck(1:ngrid)=emis(1:ngrid)sigmazplanck(1:ngrid)*zplanck(1:ngrid)
	812	fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
	813	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
	814
	815	! Test of energy conservation
	816	!----------------------------
[253]	817	if(enertest)then
[1524]	818	call planetwide_sumval(cppmassarea(:,:)zdtsw(:,:)/totarea_planet,dEtotSW)
	819	call planetwide_sumval(cppmassarea(:,:)zdtlw(:,:)/totarea_planet,dEtotLW)
[1542]	820	!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
	821	call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
	822	call planetwide_sumval((fluxsurf_lw(:)emis(:)-zplanck(:))cell_area(:)/totarea_planet,dEtotsLW)
[1524]	823	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
	824	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
	825	if (is_master) then
[1477]	826	print*,'---------------------------------------------------------------'
	827	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
	828	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
	829	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
	830	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
	831	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
	832	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
[1524]	833	endif
[1477]	834	endif ! end of 'enertest'
[253]	835
	836	endif ! of if (callrad)
	837
	838
[1477]	839
	840	! --------------------------------------------
	841	! III. Vertical diffusion (turbulent mixing) :
	842	! --------------------------------------------
	843
[253]	844	if (calldifv) then
[526]	845
[787]	846	zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
[253]	847
[1477]	848	! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
[1524]	849	if (UseTurbDiff) then
	850
[1647]	851	call turbdiff(ngrid,nlayer,nq, &
[1477]	852	ptimestep,capcal,lwrite, &
	853	pplay,pplev,zzlay,zzlev,z0, &
	854	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
	855	pdt,pdq,zflubid, &
	856	zdudif,zdvdif,zdtdif,zdtsdif, &
[1647]	857	sensibFlux,q2,zdqdif,zdqsdif, &
[1477]	858	taux,tauy,lastcall)
[594]	859
[1524]	860	else
	861
[1477]	862	zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
[594]	863
[1647]	864	call vdifc(ngrid,nlayer,nq,zpopsk, &
[1477]	865	ptimestep,capcal,lwrite, &
	866	pplay,pplev,zzlay,zzlev,z0, &
	867	pu,pv,zh,pq,tsurf,emis,qsurf, &
	868	zdh,pdq,zflubid, &
	869	zdudif,zdvdif,zdhdif,zdtsdif, &
[1524]	870	sensibFlux,q2,zdqdif,zdqsdif, &
[1477]	871	taux,tauy,lastcall)
[253]	872
[1477]	873	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
[1524]	874	zdqevap(1:ngrid,1:nlayer)=0.
[594]	875
[1477]	876	end if !end of 'UseTurbDiff'
[594]	877
[1477]	878
[1308]	879	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
	880	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
	881	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
[787]	882	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
[1297]	883
[253]	884	if (tracer) then
[1308]	885	pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
[787]	886	dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
[253]	887	end if ! of if (tracer)
	888
[1477]	889
	890	! test energy conservation
[253]	891	!-------------------------
	892	if(enertest)then
[1477]	893
[1524]	894	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
[253]	895	do ig = 1, ngrid
[1524]	896	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
	897	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
[253]	898	enddo
[1477]	899
[1542]	900	call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
[1524]	901	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
[1542]	902	call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
	903	call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
[1477]	904
[1524]	905	if (is_master) then
[1477]	906
	907	if (UseTurbDiff) then
[1524]	908	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
	909	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
[1477]	910	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
[1524]	911	else
	912	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
	913	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
	914	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
	915	end if
	916	endif ! end of 'is_master'
[1477]	917
	918	! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
	919	endif ! end of 'enertest'
[253]	920
[1477]	921	else ! calldifv
[253]	922
	923	if(.not.newtonian)then
	924
[787]	925	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
[253]	926
	927	endif
	928
[1477]	929	endif ! end of 'calldifv'
[253]	930
	931
[1477]	932	!----------------------------------
	933	! IV. Dry convective adjustment :
	934	!----------------------------------
[253]	935
	936	if(calladj) then
	937
[1308]	938	zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
	939	zduadj(1:ngrid,1:nlayer)=0.0
	940	zdvadj(1:ngrid,1:nlayer)=0.0
	941	zdhadj(1:ngrid,1:nlayer)=0.0
[253]	942
	943
[1477]	944	call convadj(ngrid,nlayer,nq,ptimestep, &
	945	pplay,pplev,zpopsk, &
	946	pu,pv,zh,pq, &
	947	pdu,pdv,zdh,pdq, &
	948	zduadj,zdvadj,zdhadj, &
	949	zdqadj)
[253]	950
[1308]	951	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
	952	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer)
	953	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
	954	zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
[1283]	955
[253]	956	if(tracer) then
[1308]	957	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
[253]	958	end if
	959
[1477]	960	! Test energy conservation
[253]	961	if(enertest)then
[1524]	962	call planetwide_sumval(cppmassarea(:,:)zdtadj(:,:)/totarea_planet,dEtot)
[1295]	963	if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
[253]	964	endif
	965
[787]	966
[1477]	967	endif ! end of 'calladj'
	968
[253]	969
[1477]	970	!---------------------------------------------
[1647]	971	! V. Specific parameterizations for tracers
[1477]	972	!---------------------------------------------
[253]	973
[1477]	974	if (tracer) then
	975
[1672]	976	! -------------------------
	977	! V.1. Chemistry
	978	! -------------------------
	979	if (callchim) then
	980
[1787]	981	! Utilisation de la moyenne zonale dans calchim
[1672]	982	zplev(:,:) = zplevbar(:,:)
	983	zplay(:,:) = zplaybar(:,:)
	984	zzlev(:,:) = zzlevbar(:,:)
	985	zzlay(:,:) = zzlaybar(:,:)
	986	ztemp(:,:) = ztfibar(:,:)
	987
	988	if (nq.gt.nmicro) then
	989	do iq = nmicro+1,nq
[1787]	990	ychim(:,:,iq-nmicro) = pq(:,:,iq) * rat_mmol(iq) ! convert to molar fraction
[1672]	991	enddo
	992	endif
	993
	994	! Condensation tendency after the transport
	995	do iq=1,nq-nmicro
	996	do l=1,nlayer
	997	do ig=1,ngrid
	998	if ( ychim(ig,l,iq).gt.qysat(l,iq) ) then
	999	dycmph(ig,l,nmicro+iq)= ( -ychim(ig,l,iq)+qysat(l,iq) ) / ptimestep
	1000	endif
	1001	enddo
	1002	enddo
	1003	enddo
	1004
	1005	if( mod(icount-1,ichim).eq.0.) then
	1006
	1007	print *, "On passe dans la chimie..."
	1008
	1009	call calchim(ngrid,nq-nmicro,ychim,nomqy,declin,zls,ctimestep, &
	1010	ztemp,zplay,zplev,zzlay,zzlev,dycchi,nlayer+70)
	1011
	1012	! JVO 2017 : NLEV = nlayer+70, en accord avec le C. Quid si nlay=/ 55 ?
	1013
	1014	endif
	1015
	1016	if (nq.gt.nmicro) then
[1787]	1017	! We convert tendencies back to mass mixing ratio
[1672]	1018	do iq=nmicro+1,nq
[1787]	1019	zdqchi(:,:,iq) = dycchi(:,:,iq-nmicro) / rat_mmol(iq)
	1020	zdqmph(:,:,iq) = dycmph(:,:,iq-nmicro) / rat_mmol(iq)
[1672]	1021	enddo
	1022
	1023	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + &
	1024	zdqchi(1:ngrid,1:nlayer,1:nq) + zdqmph(1:ngrid,1:nlayer,1:nq)
	1025
	1026	endif
	1027
	1028	endif
[1647]	1029
[1477]	1030	! ---------------
[1672]	1031	! V.3 Updates
[1477]	1032	! ---------------
[253]	1033
[1477]	1034	! Updating Atmospheric Mass and Tracers budgets.
[728]	1035	if(mass_redistrib) then
	1036
[1647]	1037	zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * zdqevap(1:ngrid,1:nlayer)
[863]	1038
	1039	do ig = 1, ngrid
[1524]	1040	zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
[863]	1041	enddo
[728]	1042
[1524]	1043	call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
	1044	call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
	1045	call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
[728]	1046
[1524]	1047	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
[1647]	1048	capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
[1524]	1049	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
	1050	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
	1051
[1308]	1052	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
[1477]	1053	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
	1054	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
	1055	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
	1056	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
[1524]	1057	pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
[1477]	1058	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
[1524]	1059
	1060	endif
[728]	1061
[1477]	1062	! -----------------------------
[1672]	1063	! V.4. Surface Tracer Update
[1477]	1064	! -----------------------------
[1297]	1065
[1647]	1066	qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
[253]	1067
[1477]	1068	! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
	1069	! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
[622]	1070	qsurf_hist(:,:) = qsurf(:,:)
[253]	1071
[1477]	1072	endif! end of if 'tracer'
[253]	1073
	1074
[1477]	1075	!------------------------------------------------
[1647]	1076	! VI. Surface and sub-surface soil temperature
[1477]	1077	!------------------------------------------------
[253]	1078
[1477]	1079
	1080	! Increment surface temperature
[1297]	1081
[1647]	1082	tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
[1477]	1083
	1084	! Compute soil temperatures and subsurface heat flux.
[253]	1085	if (callsoil) then
[787]	1086	call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, &
[1477]	1087	ptimestep,tsurf,tsoil,capcal,fluxgrd)
[253]	1088	endif
	1089
[1297]	1090
[1477]	1091	! Test energy conservation
[253]	1092	if(enertest)then
[1542]	1093	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf(:)/totarea_planet,dEtots)
[1524]	1094	if (is_master) print*,'Surface energy change =',dEtots,' W m-2'
[253]	1095	endif
	1096
	1097
[1477]	1098	!---------------------------------------------------
[1647]	1099	! VII. Perform diagnostics and write output files
[1477]	1100	!---------------------------------------------------
	1101
	1102	! Note : For output only: the actual model integration is performed in the dynamics.
	1103
	1104
[253]	1105
[1477]	1106	! Temperature, zonal and meridional winds.
[1308]	1107	zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
	1108	zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
	1109	zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
[253]	1110
[1477]	1111	! Diagnostic.
[1637]	1112	zdtdyn(1:ngrid,1:nlayer) = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep
[1308]	1113	ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
[253]	1114
[1637]	1115	zdudyn(1:ngrid,1:nlayer) = (pu(1:ngrid,1:nlayer)-zuprevious(1:ngrid,1:nlayer)) / ptimestep
	1116	zuprevious(1:ngrid,1:nlayer) = zu(1:ngrid,1:nlayer)
	1117
[253]	1118	if(firstcall)then
[1308]	1119	zdtdyn(1:ngrid,1:nlayer)=0.0
[1637]	1120	zdudyn(1:ngrid,1:nlayer)=0.0
[253]	1121	endif
	1122
[1477]	1123	! Dynamical heating diagnostic.
[253]	1124	do ig=1,ngrid
[1637]	1125	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp
[253]	1126	enddo
	1127
[1477]	1128	! Tracers.
[1308]	1129	zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
[253]	1130
[1477]	1131	! Surface pressure.
[787]	1132	ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
[253]	1133
	1134
	1135
[1477]	1136	! Surface and soil temperature information
[1542]	1137	call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
[1295]	1138	call planetwide_minval(tsurf(:),Ts2)
	1139	call planetwide_maxval(tsurf(:),Ts3)
[253]	1140	if(callsoil)then
[1542]	1141	TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
[1722]	1142	if (is_master) then
	1143	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
	1144	print*,Ts1,Ts2,Ts3,TsS
	1145	end if
[959]	1146	else
[1722]	1147	if (is_master) then
	1148	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
[1477]	1149	print*,Ts1,Ts2,Ts3
[1524]	1150	endif
[959]	1151	end if
[253]	1152
	1153
[1477]	1154	! Check the energy balance of the simulation during the run
[253]	1155	if(corrk)then
	1156
[1542]	1157	call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
	1158	call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
	1159	call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
	1160	call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
	1161	call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
[787]	1162	do ig=1,ngrid
[253]	1163	if(fluxtop_dn(ig).lt.0.0)then
	1164	print*,'fluxtop_dn has gone crazy'
	1165	print*,'fluxtop_dn=',fluxtop_dn(ig)
	1166	print*,'temp= ',pt(ig,:)
	1167	print*,'pplay= ',pplay(ig,:)
	1168	call abort
	1169	endif
	1170	end do
	1171
[787]	1172	if(ngrid.eq.1)then
[253]	1173	DYN=0.0
	1174	endif
[1524]	1175
	1176	if (is_master) then
[1477]	1177	print*,' ISR ASR OLR GND DYN [W m^-2]'
	1178	print*, ISR,ASR,OLR,GND,DYN
[1524]	1179	endif
[253]	1180
[1295]	1181	if(enertest .and. is_master)then
[651]	1182	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
	1183	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
	1184	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
[253]	1185	endif
	1186
[1295]	1187	if(meanOLR .and. is_master)then
[1216]	1188	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
[253]	1189	! to record global radiative balance
[588]	1190	open(92,file="rad_bal.out",form='formatted',position='append')
[651]	1191	write(92,*) zday,ISR,ASR,OLR
[253]	1192	close(92)
[588]	1193	open(93,file="tem_bal.out",form='formatted',position='append')
[1295]	1194	if(callsoil)then
[1524]	1195	write(93,*) zday,Ts1,Ts2,Ts3,TsS
	1196	else
	1197	write(93,*) zday,Ts1,Ts2,Ts3
	1198	endif
[253]	1199	close(93)
	1200	endif
	1201	endif
	1202
[1477]	1203	endif ! end of 'corrk'
[253]	1204
[651]	1205
[1477]	1206	! Diagnostic to test radiative-convective timescales in code.
[253]	1207	if(testradtimes)then
[588]	1208	open(38,file="tau_phys.out",form='formatted',position='append')
[253]	1209	ig=1
	1210	do l=1,nlayer
	1211	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
	1212	enddo
	1213	close(38)
[726]	1214	print*,'As testradtimes enabled,'
	1215	print*,'exiting physics on first call'
[253]	1216	call abort
	1217	endif
	1218
[1477]	1219
[1722]	1220	if (is_master) print,'--> Ls =',zls180./pi
[1477]	1221
	1222
	1223	!----------------------------------------------------------------------
[253]	1224	! Writing NetCDF file "RESTARTFI" at the end of the run
[1477]	1225	!----------------------------------------------------------------------
	1226
[253]	1227	! Note: 'restartfi' is stored just before dynamics are stored
	1228	! in 'restart'. Between now and the writting of 'restart',
	1229	! there will have been the itau=itau+1 instruction and
	1230	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	1231	! thus we store for time=time+dtvr
	1232
	1233
	1234
[1477]	1235	if(lastcall) then
	1236	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
[305]	1237
[1477]	1238	if (ngrid.ne.1) then
	1239	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
	1240
	1241	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
	1242	ptimestep,ztime_fin, &
[1647]	1243	tsurf,tsoil,emis,q2,qsurf_hist)
[1477]	1244	endif
[1216]	1245
[1477]	1246	endif ! end of 'lastcall'
[253]	1247
[861]	1248
[1477]	1249	!-----------------------------------
[253]	1250	! Saving statistics :
[1477]	1251	!-----------------------------------
[253]	1252
[1477]	1253	! Note :("stats" stores and accumulates 8 key variables in file "stats.nc"
	1254	! which can later be used to make the statistic files of the run:
	1255	! "stats") only possible in 3D runs !!!
	1256
[253]	1257
[1477]	1258	if (callstats) then
[253]	1259
[1477]	1260	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
	1261	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1262	call wstats(ngrid,"fluxsurf_lw", &
	1263	"Thermal IR radiative flux to surface","W.m-2",2, &
	1264	fluxsurf_lw)
	1265	call wstats(ngrid,"fluxtop_lw", &
	1266	"Thermal IR radiative flux to space","W.m-2",2, &
	1267	fluxtop_lw)
	1268
[253]	1269	! call wstats(ngrid,"fluxsurf_sw", &
	1270	! "Solar radiative flux to surface","W.m-2",2, &
[1477]	1271	! fluxsurf_sw_tot)
[253]	1272	! call wstats(ngrid,"fluxtop_sw", &
	1273	! "Solar radiative flux to space","W.m-2",2, &
	1274	! fluxtop_sw_tot)
[526]	1275
[253]	1276
[1477]	1277	call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
	1278	call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
	1279	call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
[1482]	1280	!call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
	1281	!call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
[1477]	1282	call wstats(ngrid,"p","Pressure","Pa",3,pplay)
	1283	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
	1284	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
	1285	call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
	1286	call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
	1287	call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
[526]	1288
[1477]	1289	if (tracer) then
	1290	do iq=1,nq
	1291	call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
	1292	call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
	1293	'kg m^-2',2,qsurf(1,iq) )
	1294
	1295	! call wstats(ngrid,trim(noms(iq))//'_reff', &
[726]	1296	! trim(noms(iq))//'_reff', &
	1297	! 'm',3,reffrad(1,1,iq))
[1477]	1298
	1299	end do
[253]	1300
[1477]	1301	endif ! end of 'tracer'
[253]	1302
[1477]	1303	if(lastcall) then
	1304	write (,) "Writing stats..."
	1305	call mkstats(ierr)
	1306	endif
	1307
	1308	endif ! end of 'callstats'
[253]	1309
	1310
[1477]	1311	!-----------------------------------------------------------------------------------------------------
	1312	! OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
	1313	!
	1314	! Note 1 : output with period "ecritphy", set in "run.def"
	1315	!
	1316	! Note 2 : writediagfi can also be called from any other subroutine for any variable,
	1317	! but its preferable to keep all the calls in one place ...
	1318	!-----------------------------------------------------------------------------------------------------
[253]	1319
	1320
[1477]	1321	call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
	1322	call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
	1323	call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
	1324	call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
	1325	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1326	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
	1327	call writediagfi(ngrid,"temp","temperature","K",3,zt)
	1328	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
	1329	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1330	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
	1331	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
	1332	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
	1333
[965]	1334	! Subsurface temperatures
[969]	1335	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1336	! call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
[965]	1337
[1477]	1338
[1297]	1339
[1477]	1340	! Total energy balance diagnostics
	1341	if(callrad.and.(.not.newtonian))then
	1342
[1482]	1343	!call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
	1344	!call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
[1477]	1345	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
	1346	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
	1347	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
	1348	call writediagfi(ngrid,"shad","rings"," ", 2, fract)
	1349
[1016]	1350	! call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
	1351	! call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
	1352	! call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
	1353	! call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
	1354	! call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
	1355	! call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
[1477]	1356
[1647]	1357
	1358	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
[1477]	1359
	1360	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
	1361
	1362	endif ! end of 'callrad'
[1524]	1363
[1477]	1364	if(enertest) then
	1365
[1524]	1366	if (calldifv) then
[1477]	1367
[1524]	1368	call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
[1477]	1369	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
	1370
[1524]	1371	! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
	1372	! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
	1373	! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
	1374
	1375	endif
[1477]	1376
[1524]	1377	if (corrk) then
	1378	call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
	1379	call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
	1380	endif
[1477]	1381
	1382	endif ! end of 'enertest'
[253]	1383
[1663]	1384	! Temporary inclusions for winds diagnostics.
	1385	call writediagfi(ngrid,"zdudif","Turbdiff tend. zon. wind","m s-2",3,zdudif)
	1386	call writediagfi(ngrid,"zdudyn","Dyn. tend. zon. wind","m s-2",3,zdudyn)
[1477]	1387
	1388	! Temporary inclusions for heating diagnostics.
[1663]	1389	call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
	1390	call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
	1391	call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
	1392	call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
[1477]	1393
	1394	! For Debugging.
[1663]	1395	call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
[1477]	1396
[253]	1397
[1477]	1398	! Output tracers.
	1399	if (tracer) then
	1400
	1401	do iq=1,nq
	1402	call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
	1403	call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
	1404	'kg m^-2',2,qsurf_hist(1,iq) )
	1405
[787]	1406	! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
[1477]	1407	! 'kg m^-2',2,qsurf(1,iq) )
[253]	1408
[1477]	1409	enddo ! end of 'nq' loop
	1410
	1411	endif ! end of 'tracer'
[253]	1412
[1477]	1413
	1414	! Output spectrum.
[526]	1415	if(specOLR.and.corrk)then
[728]	1416	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
	1417	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
[526]	1418	endif
[253]	1419
[1622]	1420	! XIOS outputs
	1421	#ifdef CPP_XIOS
	1422	! Send fields to XIOS: (NB these fields must also be defined as
	1423	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
[1626]	1424	CALL send_xios_field("ls",zls)
	1425
[1622]	1426	CALL send_xios_field("ps",ps)
	1427	CALL send_xios_field("area",cell_area)
	1428
	1429	CALL send_xios_field("temperature",zt)
	1430	CALL send_xios_field("u",zu)
	1431	CALL send_xios_field("v",zv)
	1432
	1433	#endif
	1434
[253]	1435	icount=icount+1
	1436
[751]	1437	end subroutine physiq
[1549]	1438
	1439	end module physiq_mod

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

Download in other formats:

Original Format