Context Navigation

physiq_mod.F90 @ 1915

Last change on this file since 1915 was 1915, checked in by jvatant, 7 years ago

Add main current send_xios_field outputs
--JVO

Property svn:executable set to *
Property svn:mergeinfo set to
/trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 merged eligible

File size: 63.2 KB

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

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

Download in other formats:

Original Format