Context Navigation

physiq_mod.F90 @ 1862

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

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

Download in other formats:

Original Format