Context Navigation

physiq_mod.F90 @ 3595

Last change on this file since 3595 was 3497, checked in by debatzbr, 3 months ago

Add AC6H6 condensation in the microphysics

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

File size: 85.5 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
[3318]	17	use radcommon_h, only: sigma, gzlat, grav, BWNV, WAVEI, WAVEV
[1958]	18	use comchem_h, only: nkim, cnames, nlaykim_up, ykim_up, ykim_tot, botCH4
[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
[3083]	22	use datafile_mod, only: datadir, corrkdir, banddir, haze_opt_file, nudging_file
	23	use geometry_mod, only: latitude, latitude_deg, 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
[2366]	35	#ifndef MESOSCALE
[1914]	36	use logic_mod, only: moyzon_ch
[1966]	37	use moyzon_mod, only: zphibar, zphisbar, zplevbar, zplaybar, &
	38	zzlevbar, zzlaybar, ztfibar, zqfibar
[2366]	39	#endif
[1397]	40	use callkeys_mod
[2291]	41	use phys_state_var_mod
	42	use turb_mod, only : q2,sensibFlux,turb_resolved
	43	#ifndef MESOSCALE
[1622]	44	use vertical_layers_mod, only: presnivs, pseudoalt
[1896]	45	use mod_phys_lmdz_omp_data, ONLY: is_omp_master
[2291]	46	#else
	47	use comm_wrf, only : comm_HR_SW, comm_HR_LW, &
	48	comm_FLUXTOP_DN,comm_FLUXABS_SW,&
	49	comm_FLUXTOP_LW,comm_FLUXSURF_SW,&
	50	comm_FLUXSURF_LW,comm_FLXGRD
	51	#endif
[1623]	52	#ifdef CPP_XIOS
[1622]	53	use xios_output_mod, only: initialize_xios_output, &
	54	update_xios_timestep, &
	55	send_xios_field
[1896]	56	use wxios, only: wxios_context_init, xios_context_finalize
[1623]	57	#endif
[1926]	58	use muphy_diag
[253]	59	implicit none
	60
	61
	62	!==================================================================
	63	!
	64	! Purpose
	65	! -------
	66	! Central subroutine for all the physics parameterisations in the
	67	! universal model. Originally adapted from the Mars LMDZ model.
	68	!
	69	! The model can be run without or with tracer transport
	70	! depending on the value of "tracer" in file "callphys.def".
	71	!
	72	!
	73	! It includes:
	74	!
[1477]	75	! I. Initialization :
	76	! I.1 Firstcall initializations.
	77	! I.2 Initialization for every call to physiq.
[253]	78	!
[1477]	79	! II. Compute radiative transfer tendencies (longwave and shortwave) :
	80	! II.a Option 1 : Call correlated-k radiative transfer scheme.
	81	! II.b Option 2 : Call Newtonian cooling scheme.
	82	! II.c Option 3 : Atmosphere has no radiative effect.
	83	!
	84	! III. Vertical diffusion (turbulent mixing) :
	85	!
	86	! IV. Dry Convective adjustment :
	87	!
[1647]	88	! V. Tracers
[1926]	89	! V.1. Microphysics
	90	! V.2. Chemistry
[1672]	91	! V.3. Updates (pressure variations, surface budget).
	92	! V.4. Surface Tracer Update.
[1477]	93	!
[1647]	94	! VI. Surface and sub-surface soil temperature.
[1477]	95	!
[1647]	96	! VII. Perform diagnostics and write output files.
[1477]	97	!
	98	!
[253]	99	! arguments
	100	! ---------
	101	!
[1477]	102	! INPUT
[253]	103	! -----
[1477]	104	!
[253]	105	! ngrid Size of the horizontal grid.
	106	! nlayer Number of vertical layers.
[1477]	107	! nq Number of advected fields.
	108	! nametrac Name of corresponding advected fields.
	109	!
	110	! firstcall True at the first call.
	111	! lastcall True at the last call.
	112	!
	113	! pday Number of days counted from the North. Spring equinoxe.
	114	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
	115	! ptimestep timestep (s).
	116	!
	117	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa).
	118	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
	119	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2.s-2).
	120	!
	121	! pu(ngrid,nlayer) u, zonal component of the wind (ms-1).
	122	! pv(ngrid,nlayer) v, meridional component of the wind (ms-1).
	123	!
	124	! pt(ngrid,nlayer) Temperature (K).
	125	!
	126	! pq(ngrid,nlayer,nq) Advected fields.
	127	!
[1216]	128	! pudyn(ngrid,nlayer) \
[253]	129	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
[1477]	130	! ptdyn(ngrid,nlayer) / corresponding variables.
[253]	131	! pqdyn(ngrid,nlayer,nq) /
[1312]	132	! flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
[253]	133	!
[1477]	134	! OUTPUT
[253]	135	! ------
	136	!
[1308]	137	! pdu(ngrid,nlayer) \
	138	! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding
	139	! pdt(ngrid,nlayer) / variables due to physical processes.
	140	! pdq(ngrid,nlayer) /
[3340]	141	! pdpsrf(ngrid) /
[253]	142	!
	143	!
	144	! Authors
	145	! -------
[1524]	146	! Frederic Hourdin 15/10/93
	147	! Francois Forget 1994
	148	! Christophe Hourdin 02/1997
[253]	149	! Subroutine completely rewritten by F. Forget (01/2000)
	150	! Water ice clouds: Franck Montmessin (update 06/2003)
	151	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
	152	! New correlated-k radiative scheme: R. Wordsworth (2009)
	153	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
	154	! Improved water cycle: R. Wordsworth / B. Charnay (2010)
	155	! To F90: R. Wordsworth (2010)
[594]	156	! New turbulent diffusion scheme: J. Leconte (2012)
[716]	157	! Loops converted to F90 matrix format: J. Leconte (2012)
[787]	158	! No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
[1477]	159	! Purge of the code : M. Turbet (2015)
[1672]	160	! Fork for Titan : J. Vatant d'Ollone (2017)
	161	! + clean of all too-generic (ocean, water, co2 ...) routines
	162	! + Titan's chemistry
[3318]	163	! Microphysical moment model: J.Burgalat / B. de Batz de Trenquelléon (2022-2023)
	164	! Optics for haze and clouds: B. de Batz de Trenquelléon (2023)
[1647]	165	!============================================================================================
[253]	166
[3083]	167	! ---------------
	168	! 0. DECLARATIONS
	169	! ---------------
[253]	170
[1670]	171	include "netcdf.inc"
[253]	172
	173	! Arguments :
	174	! -----------
	175
[1477]	176	! INPUTS:
[253]	177	! -------
	178
[1477]	179	integer,intent(in) :: ngrid ! Number of atmospheric columns.
	180	integer,intent(in) :: nlayer ! Number of atmospheric layers.
	181	integer,intent(in) :: nq ! Number of tracers.
[1984]	182	character*30,intent(in) :: nametrac(nq) ! Names of the tracers taken from dynamics.
[1477]	183
	184	logical,intent(in) :: firstcall ! Signals first call to physics.
	185	logical,intent(in) :: lastcall ! Signals last call to physics.
	186
	187	real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0.
	188	real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
	189	real,intent(in) :: ptimestep ! Physics timestep (s).
	190	real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
	191	real,intent(in) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
	192	real,intent(in) :: pphi(ngrid,nlayer) ! Geopotential at mid-layer (m2s-2).
	193	real,intent(in) :: pu(ngrid,nlayer) ! Zonal wind component (m/s).
	194	real,intent(in) :: pv(ngrid,nlayer) ! Meridional wind component (m/s).
	195	real,intent(in) :: pt(ngrid,nlayer) ! Temperature (K).
	196	real,intent(in) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
	197	real,intent(in) :: flxw(ngrid,nlayer) ! Vertical mass flux (ks/s) at lower boundary of layer
[253]	198
[1477]	199	! OUTPUTS:
[253]	200	! --------
	201
[1477]	202	! Physical tendencies :
	203
	204	real,intent(out) :: pdu(ngrid,nlayer) ! Zonal wind tendencies (m/s/s).
	205	real,intent(out) :: pdv(ngrid,nlayer) ! Meridional wind tendencies (m/s/s).
	206	real,intent(out) :: pdt(ngrid,nlayer) ! Temperature tendencies (K/s).
	207	real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
	208	real,intent(out) :: pdpsrf(ngrid) ! Surface pressure tendency (Pa/s).
	209
[253]	210	! Local saved variables:
	211	! ----------------------
[1622]	212	integer,save :: day_ini ! Initial date of the run (sol since Ls=0).
	213	integer,save :: icount ! Counter of calls to physiq during the run.
[2366]	214	!$OMP THREADPRIVATE(day_ini,icount)
[1622]	215
[253]	216	! Local variables :
	217	! -----------------
	218
[1477]	219	real zh(ngrid,nlayer) ! Potential temperature (K).
	220	real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
	221
[3083]	222	integer l,ig,ierr,iq,nw,isoil,ilat,lat_idx,i,j
	223
[1477]	224	! FOR DIAGNOSTIC :
[253]	225
[1477]	226	real zls ! Solar longitude (radians).
	227	real zlss ! Sub solar point longitude (radians).
	228	real zday ! Date (time since Ls=0, calculated in sols).
[2097]	229	real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers (ref : local surf).
	230	real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries (ref : local surf).
	231	real zzlay_eff(ngrid,nlayer) ! Effective altitude at the middle of the atmospheric layers (ref : geoid ).
	232	real zzlev_eff(ngrid,nlayer+1) ! Effective altitude at the atmospheric layer boundaries ( ref : geoid ).
	233
[1477]	234	! TENDENCIES due to various processes :
[253]	235
[1477]	236	! For Surface Temperature : (K/s)
	237	real zdtsurf(ngrid) ! Cumulated tendencies.
	238	real zdtsurfmr(ngrid) ! Mass_redistribution routine.
[3083]	239	real zdtsurfevap(ngrid) ! Evaporation.
[1477]	240	real zdtsdif(ngrid) ! Turbdiff/vdifc routines.
	241
	242	! For Atmospheric Temperatures : (K/s)
[3083]	243	real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	244	real zdtmr(ngrid,nlayer) ! Mass_redistribution routine.
	245	real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
	246	real zdtlc(ngrid,nlayer) ! Condensation heating rate.
[1477]	247
	248	! For Surface Tracers : (kg/m2/s)
	249	real dqsurf(ngrid,nq) ! Cumulated tendencies.
	250	real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines.
	251	real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine.
	252
	253	! For Tracers : (kg/kg_of_air/s)
	254	real zdqadj(ngrid,nlayer,nq) ! Convadj routine.
	255	real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines.
	256	real zdqevap(ngrid,nlayer) ! Turbdiff routine.
	257	real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine.
[1672]	258
[3090]	259	real zdqchi(ngrid,nlayer,nq) ! Chemical tendency (chemistry routine).
[1795]	260
	261	real zdqmufi(ngrid,nlayer,nq) ! Microphysical tendency.
[1926]	262
	263	real zdqfibar(ngrid,nlayer,nq) ! For 2D chemistry
	264	real zdqmufibar(ngrid,nlayer,nq) ! For 2D chemistry
[1477]	265
	266	! For Winds : (m/s/s)
	267	real zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! Convadj routine.
	268	real zdumr(ngrid,nlayer),zdvmr(ngrid,nlayer) ! Mass_redistribution routine.
	269	real zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	270	real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	271	real zdhadj(ngrid,nlayer) ! Convadj routine.
[3090]	272	real zdundg(ngrid,nlayer) ! Nudging for zonal wind.
[1477]	273
	274	! For Pressure and Mass :
	275	real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
	276	real zdmassmr_col(ngrid) ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
	277	real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s).
	278
	279
	280	! Local variables for LOCAL CALCULATIONS:
	281	! ---------------------------------------
[787]	282	real zflubid(ngrid)
[1308]	283	real zplanck(ngrid),zpopsk(ngrid,nlayer)
[253]	284	real ztim1,ztim2,ztim3, z1,z2
	285	real ztime_fin
[1308]	286	real zdh(ngrid,nlayer)
[1194]	287	real gmplanet
[1297]	288	real taux(ngrid),tauy(ngrid)
[3340]	289	real factlat
	290	real zundg(nlayer)
[1194]	291
[253]	292
[1477]	293	! local variables for DIAGNOSTICS : (diagfi & stat)
	294	! -------------------------------------------------
	295	real ps(ngrid) ! Surface Pressure.
	296	real zt(ngrid,nlayer) ! Atmospheric Temperature.
	297	real zu(ngrid,nlayer),zv(ngrid,nlayer) ! Zonal and Meridional Winds.
	298	real zq(ngrid,nlayer,nq) ! Atmospheric Tracers.
	299	real zdtadj(ngrid,nlayer) ! Convadj Diagnostic.
	300	real zdtdyn(ngrid,nlayer) ! Dynamical Heating (K/s).
[1637]	301	real zdudyn(ngrid,nlayer) ! Dynamical Zonal Wind tendency (m.s-2).
[253]	302
[3090]	303	real zhorizwind(ngrid,nlayer) ! Horizontal Wind (sqrt(uu+vv))
[2002]	304
[1477]	305	real vmr(ngrid,nlayer) ! volume mixing ratio
[253]	306	real time_phys
[597]	307
[1477]	308	real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
	309
[3083]	310	! To test energy conservation (RW+JL)
[1308]	311	real mass(ngrid,nlayer),massarea(ngrid,nlayer)
[651]	312	real dEtot, dEtots, AtmToSurf_TurbFlux
[959]	313	real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
[1315]	314	!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
[1308]	315	real dEzRadsw(ngrid,nlayer),dEzRadlw(ngrid,nlayer),dEzdiff(ngrid,nlayer)
[787]	316	real dEdiffs(ngrid),dEdiff(ngrid)
[1477]	317
[3083]	318	! JL12 : conservation test for mean flow kinetic energy has been disabled temporarily
[1295]	319	real dItot, dItot_tmp, dVtot, dVtot_tmp
	320	real dWtot, dWtot_tmp, dWtots, dWtots_tmp
[1477]	321
[1482]	322	! For Clear Sky Case.
	323	real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxsurfabs_sw1(ngrid) ! For SW/LW flux.
	324	real fluxtop_lw1(ngrid), fluxabs_sw1(ngrid) ! For SW/LW flux.
	325	real albedo_equivalent1(ngrid) ! For Equivalent albedo calculation.
[1647]	326	real tf, ntf
[253]	327
[1822]	328	! Miscellaneous :
	329	character(len=10) :: tmp1
	330	character(len=10) :: tmp2
[2133]	331	character*2 :: str2
[253]	332
[2366]	333	#ifndef MESOSCALE
	334
[1903]	335	! Local variables for Titan chemistry and microphysics
	336	! ----------------------------------------------------
	337
[1914]	338	real,save :: ctimestep ! Chemistry timestep (s)
	339	!$OMP THREADPRIVATE(ctimestep)
[1672]	340
[3083]	341	! Chemical tracers in molar fraction [mol/mol]
	342	real, dimension(ngrid,nlayer,nkim) :: ychim
	343	real, dimension(ngrid,nlayer,nkim) :: ychimbar ! For 2D chemistry
[1672]	344
[3083]	345	! Molar fraction tendencies (chemistry, condensation and evaporation) for tracers [mol/mol/s]
	346	real, dimension(ngrid,nlayer,nq) :: dyccond ! Condensation rate. NB : for all tracers, as we want to use indx on it.
	347	real, dimension(ngrid,nlayer,size(ices_indx)) :: dmuficond ! Condensation rate from microphysics [kg/kg/s].
	348	real, dimension(ngrid,nlayer,nq) :: dyccondbar ! For 2D chemistry
	349	real, dimension(ngrid) :: dycevapCH4 ! Surface "pseudo-evaporation" rate for CH4.
[1672]	350
[3083]	351	! Saturation profiles [mol/mol]
	352	real, dimension(ngrid,nlayer,nkim) :: ysat
	353
[3090]	354	! Temporary fraction of CH4 [mol/mol]
[3083]	355	real, dimension(ngrid,nlayer) :: tpq_CH4
[1903]	356
[3083]	357	real :: i2e(ngrid,nlayer) ! int 2 ext factor (X.kg-1 -> X.m-3)
[1789]	358
[2366]	359	#ifdef USE_QTEST
[1897]	360	real,save,dimension(:,:,:), allocatable :: tpq ! Tracers for decoupled microphysical tests ( temporary in 01/18 )
	361	!$OMP THREADPRIVATE(tpq)
[1915]	362	real,dimension(ngrid,nlayer,nq) :: dtpq ! (temporary in 01/18)
[2366]	363	#endif
[1897]	364
[2242]	365	logical file_ok
[1897]	366
[1672]	367	!-----------------------------------------------------------------------------
[1795]	368	! Interface to calmufi
	369	! --> needed in order to pass assumed-shape arrays. Otherwise we must put calmufi in a module
[1897]	370	! (to have an explicit interface generated by the compiler).
[1795]	371	! Or one can put calmufi in MMP_GCM module (in muphytitan).
	372	INTERFACE
[1947]	373	SUBROUTINE calmufi(dt, plev, zlev, play, zlay, g3d, temp, pq, zdqfi, zdq)
[1902]	374	REAL(kind=8), INTENT(IN) :: dt !! Physics timestep (s).
[1795]	375	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: plev !! Pressure levels (Pa).
	376	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlev !! Altitude levels (m).
	377	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: play !! Pressure layers (Pa).
	378	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlay !! Altitude at the center of each layer (m).
[1947]	379	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: g3d !! Latitude-Altitude depending gravitational acceleration (m.s-2).
[1795]	380	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp !! Temperature at the center of each layer (K).
[2242]	381	REAL(kind=8), DIMENSION(:,:,:), INTENT(IN) :: pq !! Tracers ($X.kg^{-1}}$).
	382	REAL(kind=8), DIMENSION(:,:,:), INTENT(IN) :: zdqfi !! Tendency from former processes for tracers ($X.kg^{-1}}$).
	383	REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq !! Microphysical tendency for tracers ($X.kg^{-1}}$).
[1795]	384	END SUBROUTINE calmufi
	385	END INTERFACE
[2366]	386
	387	#endif
[1672]	388
[1477]	389	!==================================================================================================
[253]	390
	391	! -----------------
[1477]	392	! I. INITIALISATION
	393	! -----------------
[253]	394
[1477]	395	! --------------------------------
	396	! I.1 First Call Initialisation.
	397	! --------------------------------
[253]	398	if (firstcall) then
[2366]	399
	400	#ifdef USE_QTEST
[1897]	401	allocate(tpq(ngrid,nlayer,nq))
	402	tpq(:,:,:) = pq(:,:,:)
[2366]	403	#endif
[1844]	404	! Initialisation of nmicro as well as tracers names, indexes ...
[1903]	405	if (ngrid.ne.1) then ! Already done in rcm1d
[1926]	406	call initracer2(nq,nametrac) ! WARNING JB (27/03/2018): should be wrapped in an OMP SINGLE statement (see module notes)
[1844]	407	endif
[1795]	408
[2366]	409	! Allocate saved arrays (except for 1D model, where this has already been done)
	410	#ifndef MESOSCALE
	411	if (ngrid>1) call phys_state_var_init(nq)
[2291]	412	#endif
	413
[1477]	414	! Variables set to 0
[253]	415	! ~~~~~~~~~~~~~~~~~~
[2242]	416	dtrad(:,:) = 0.D0
[3083]	417	zdtlc(:,:) = 0.D0
[2242]	418	fluxrad(:) = 0.D0
	419	zdtsw(:,:) = 0.D0
	420	zdtlw(:,:) = 0.D0
[3083]	421	zpopthi(:,:,:,:) = 0.D0
	422	zpopthv(:,:,:,:) = 0.D0
[3318]	423	zpoptti(:,:,:,:) = 0.D0
	424	zpopttv(:,:,:,:) = 0.D0
[726]	425
[1822]	426	! Initialize setup for correlated-k radiative transfer
[2242]	427	! JVO 17 : Was in callcorrk firstcall, but we need spectral intervals for microphysics.
	428	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1822]	429
	430	if (corrk) then
	431
	432	call system('rm -f surf_vals_long.out')
	433
	434	write( tmp1, '(i3)' ) L_NSPECTI
	435	write( tmp2, '(i3)' ) L_NSPECTV
	436	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
	437	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
	438
	439	call setspi ! Basic infrared properties.
	440	call setspv ! Basic visible properties.
	441	call sugas_corrk ! Set up gaseous absorption properties.
	442
[2133]	443	OLR_nu(:,:) = 0.D0
	444	OSR_nu(:,:) = 0.D0
[1822]	445
[2133]	446	int_dtaui(:,:,:) = 0.D0
	447	int_dtauv(:,:,:) = 0.D0
[2372]	448
	449	#ifndef MESOSCALE
[2242]	450	IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
[3083]	451	haze_opt_file=trim(datadir)//'/optical_tables/HAZE_OPTIC_'//trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))//'.DAT'
[2242]	452	inquire(file=trim(haze_opt_file),exist=file_ok)
	453	if(.not.file_ok) then
	454	write(,) 'The file ',TRIM(haze_opt_file),' with the haze optical properties'
	455	write(,) 'was not found by optci.F90 ! Check in ', TRIM(datadir)
	456	write(,) 'that you have the one corresponding to the given spectral resolution !!'
	457	write(,) 'Meanwhile I abort ...'
	458	call abort
	459	endif
	460	ENDIF
[2372]	461	#endif
	462
[1822]	463	endif
	464
[2366]	465	#ifndef MESOSCALE
[1966]	466	! Initialize names and timestep for chemistry
	467	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1795]	468
[3090]	469	if (callchim) then
[1795]	470
[3090]	471	if (moyzon_ch .and. ngrid.eq.1) then
[1914]	472	print *, "moyzon_ch=",moyzon_ch," and ngrid=1"
	473	print *, "Please desactivate zonal mean for 1D !"
	474	stop
	475	endif
[2366]	476
[1903]	477	! Chemistry timestep
[1795]	478	ctimestep = ptimestep*REAL(ichim)
	479
[1477]	480	endif
[253]	481
[1795]	482	! Initialize microphysics.
	483	! ~~~~~~~~~~~~~~~~~~~~~~~~
[726]	484
[3090]	485	IF (callmufi) THEN
[2242]	486	! WARNING JB (27/03/2018): inimufi should be wrapped in an OMP SINGLE statement.
[1926]	487	call inimufi(ptimestep)
[1795]	488
[1926]	489	! initialize microphysics diagnostics arrays.
	490	call ini_diag_arrays(ngrid,nlayer,nice)
	491
[1795]	492	ENDIF
[2366]	493	#endif
[1795]	494
[1896]	495	#ifdef CPP_XIOS
	496	! Initialize XIOS context
	497	write(,) "physiq: call wxios_context_init"
	498	CALL wxios_context_init
	499	#endif
	500
[1477]	501	! Read 'startfi.nc' file.
[253]	502	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2291]	503	#ifndef MESOSCALE
[2366]	504	call phyetat0(startphy_file, &
	505	ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
[1789]	506	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,tankCH4)
[2291]	507	#else
	508	emis(:)=0.0
	509	q2(:,:)=0.0
	510	qsurf(:,:)=0.0
[2366]	511	tankCH4(:)=0.0
[2291]	512	day_ini = pday
	513	#endif
	514
	515	#ifndef MESOSCALE
[1670]	516	if (.not.startphy_file) then
	517	! additionnal "academic" initialization of physics
	518	if (is_master) write(,) "Physiq: initializing tsurf(:) to pt(:,1) !!"
	519	tsurf(:)=pt(:,1)
	520	if (is_master) write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
	521	do isoil=1,nsoilmx
[1914]	522	tsoil(:,isoil)=tsurf(:)
[1670]	523	enddo
	524	if (is_master) write(,) "Physiq: initializing day_ini to pdat !"
	525	day_ini=pday
	526	endif
[2291]	527	#endif
[253]	528
	529	if (pday.ne.day_ini) then
	530	write(,) "ERROR in physiq.F90:"
	531	write(,) "bad synchronization between physics and dynamics"
	532	write(,) "dynamics day: ",pday
	533	write(,) "physics day: ",day_ini
	534	stop
	535	endif
	536	write (,) 'In physiq day_ini =', day_ini
	537
[1482]	538	! Initialize albedo calculation.
	539	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2242]	540	albedo(:,:)=0.D0
	541	albedo_bareground(:)=0.D0
[1647]	542	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground)
[1482]	543
	544	! Initialize orbital calculation.
	545	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[253]	546	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
	547
[2366]	548
[253]	549	if(tlocked)then
	550	print*,'Planet is tidally locked at resonance n=',nres
	551	print*,'Make sure you have the right rotation rate!!!'
	552	endif
	553
[1297]	554
[1477]	555	! Initialize soil.
	556	! ~~~~~~~~~~~~~~~~
[253]	557	if (callsoil) then
[1477]	558
[787]	559	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
[1477]	560	ptimestep,tsurf,tsoil,capcal,fluxgrd)
[1297]	561
[1477]	562	else ! else of 'callsoil'.
	563
[253]	564	print*,'WARNING! Thermal conduction in the soil turned off'
[918]	565	capcal(:)=1.e6
[952]	566	fluxgrd(:)=intheat
	567	print*,'Flux from ground = ',intheat,' W m^-2'
[1477]	568
	569	endif ! end of 'callsoil'.
	570
[253]	571	icount=1
[1477]	572
[253]	573
[1477]	574	! Initialize surface history variable.
[253]	575	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[728]	576	qsurf_hist(:,:)=qsurf(:,:)
[253]	577
[1637]	578	! Initialize variable for dynamical heating and zonal wind tendency diagnostic
	579	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[253]	580	ztprevious(:,:)=pt(:,:)
[1637]	581	zuprevious(:,:)=pu(:,:)
[253]	582
	583
[1477]	584	if(meanOLR)then
	585	call system('rm -f rad_bal.out') ! to record global radiative balance.
	586	call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
	587	call system('rm -f h2o_bal.out') ! to record global hydrological balance.
[253]	588	endif
	589
[3318]	590	! Read NewNudging.dat and initialize the nudging for pu
	591	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3083]	592	if (nudging_u) then
[3318]	593	nudging_file=trim(datadir)//'/NewNudging.dat'
[3083]	594	inquire(file=trim(nudging_file),exist=file_ok)
	595	if(.not.file_ok) then
	596	write(,) 'ERROR : The file ',TRIM(nudging_file),' was not found by physiq_mod.F90 ! Check in ', TRIM(datadir)
	597	write(,) 'Meanwhile I abort ...'
	598	call abort
	599	endif
	600
	601	open(88,file=nudging_file,form='formatted')
	602	do ilat = 1, 49
	603	read(88,*) u_ref(ilat,:)
	604	enddo
	605	close(88)
	606	endif
	607
[2291]	608	#ifndef MESOSCALE
[2366]	609	if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
[1542]	610	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
[2366]	611	ptimestep,pday+nday,time_phys,cell_area, &
	612	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
[1904]	613	endif
[2366]	614	#endif
[1843]	615
[1622]	616	! XIOS outputs
	617	#ifdef CPP_XIOS
	618
	619	write(,) "physiq: call initialize_xios_output"
	620	call initialize_xios_output(pday,ptime,ptimestep,daysec, &
[3318]	621	presnivs,pseudoalt,wavei,wavev)
[1622]	622	#endif
[1896]	623	write(,) "physiq: end of firstcall"
[1477]	624	endif ! end of 'firstcall'
[253]	625
[1477]	626	! ------------------------------------------------------
	627	! I.2 Initializations done at every physical timestep:
	628	! ------------------------------------------------------
	629
[1622]	630	#ifdef CPP_XIOS
	631	! update XIOS time/calendar
	632	call update_xios_timestep
	633	#endif
	634
[1477]	635	! Initialize various variables
	636	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1647]	637
[3083]	638	pdt(:,:) = 0.D0
	639	zdtsurf(:) = 0.D0
	640	zdtsurfevap(:) = 0.D0
	641	pdq(:,:,:) = 0.D0
	642	dqsurf(:,:) = 0.D0
	643	pdu(:,:) = 0.D0
	644	pdv(:,:) = 0.D0
	645	pdpsrf(:) = 0.D0
	646	zflubid(:) = 0.D0
	647	taux(:) = 0.D0
	648	tauy(:) = 0.D0
[253]	649
[1477]	650	zday=pday+ptime ! Compute time, in sols (and fraction thereof).
[1297]	651
[1477]	652	! Compute Stellar Longitude (Ls), and orbital parameters.
	653	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[253]	654	if (season) then
	655	call stellarlong(zday,zls)
	656	else
[2366]	657	call stellarlong(noseason_day,zls)
[253]	658	end if
	659
[1329]	660	call orbite(zls,dist_star,declin,right_ascen)
[1477]	661
[1329]	662	if (tlocked) then
	663	zlss=Mod(-(2.pi(zday/year_day)nres - right_ascen),2.pi)
	664	elseif (diurnal) then
[1524]	665	zlss=-2.pi(zday-.5)
[1329]	666	else if(diurnal .eqv. .false.) then
	667	zlss=9999.
	668	endif
[1194]	669
	670
[1477]	671	! Compute variations of g with latitude (oblate case).
	672	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	673	if (oblate .eqv. .false.) then
[1947]	674	gzlat(:,:) = g
[1477]	675	else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
[1947]	676	print*,'I need values for flatten, J2, Rmean and MassPlanet to compute gzlat (else set oblate=.false.)'
[1477]	677	call abort
	678	else
	679	gmplanet = MassPlanetgrav1e24
	680	do ig=1,ngrid
[1947]	681	gzlat(ig,:)= gmplanet/(Rmean*2) (1.D0 + 0.75 J2 - 2.0flatten/3. + (2.flatten - 15./4. J2) * cos(2. * (pi/2. - latitude(ig))))
[1477]	682	end do
	683	endif
[1947]	684
	685	! Compute altitudes with the geopotential coming from the dynamics.
	686	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1297]	687
[1947]	688	if (eff_gz .eqv. .false.) then
	689
	690	do l=1,nlayer
[2097]	691	zzlay(:,l) = pphi(:,l) / gzlat(:,l) ! Reference = local surface
[1947]	692	enddo
	693
	694	else ! In this case we consider variations of g with altitude
	695
	696	do l=1,nlayer
[2097]	697	zzlay(:,l) = gradrad / ( g*rad - ( pphi(:,l) + phisfi(:) )) - rad
[1947]	698	gzlat(:,l) = gradrad / ( rad + zzlay(:,l) )**2
	699	end do
	700
	701	endif ! if eff_gz
[728]	702
[1914]	703	zzlev(:,1)=0.
	704	zzlev(:,nlayer+1)=1.e7 ! Dummy top of last layer above 10000 km...
[2097]	705	! JVO 19 : This altitude is indeed dummy for the GCM and fits ptop=0
	706	! but for upper chemistry that's a pb -> we anyway redefine it just after ..
[728]	707
[253]	708	do l=2,nlayer
	709	do ig=1,ngrid
	710	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
	711	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
	712	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	713	enddo
[1477]	714	enddo
[253]	715
[2097]	716	! Effective altitudes ( eg needed for chemistry ) with correct g, and with reference to the geoid
	717	! JVO 19 : We shall always have correct altitudes in chemistry no matter what's in physics
	718	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2366]	719	#ifndef MESOSCALE
[2097]	720	if (moyzon_ch) then ! Zonal averages
[1672]	721
[2097]	722	zzlaybar(1,:)=gradrad/(g*rad-(zphibar(1,:)+zphisbar(1)))-rad ! reference = geoid
[1947]	723	zzlevbar(1,1)=zphisbar(1)/g
[1672]	724	DO l=2,nlayer
[1947]	725	z1=(zplaybar(1,l-1)+zplevbar(1,l))/(zplaybar(1,l-1)-zplevbar(1,l))
[1672]	726	z2=(zplevbar(1,l) +zplaybar(1,l))/(zplevbar(1,l) -zplaybar(1,l))
	727	zzlevbar(1,l)=(z1zzlaybar(1,l-1)+z2zzlaybar(1,l))/(z1+z2)
	728	ENDDO
	729	zzlevbar(1,nlayer+1)=zzlaybar(1,nlayer)+(zzlaybar(1,nlayer)-zzlevbar(1,nlayer))
	730
	731	DO ig=2,ngrid
[1947]	732	if (latitude(ig).ne.latitude(ig-1)) then
	733	DO l=1,nlayer
[2097]	734	zzlaybar(ig,l)=gradrad/(g*rad-(zphibar(ig,l)+zphisbar(ig)))-rad
[1672]	735	ENDDO
	736	zzlevbar(ig,1)=zphisbar(ig)/g
	737	DO l=2,nlayer
[1947]	738	z1=(zplaybar(ig,l-1)+zplevbar(ig,l))/ (zplaybar(ig,l-1)-zplevbar(ig,l))
[1672]	739	z2=(zplevbar(ig,l) +zplaybar(ig,l))/(zplevbar(ig,l) -zplaybar(ig,l))
	740	zzlevbar(ig,l)=(z1zzlaybar(ig,l-1)+z2zzlaybar(ig,l))/(z1+z2)
	741	ENDDO
[1947]	742	zzlevbar(ig,nlayer+1)=zzlaybar(ig,nlayer)+(zzlaybar(ig,nlayer)-zzlevbar(ig,nlayer))
[1672]	743	else
	744	zzlaybar(ig,:)=zzlaybar(ig-1,:)
	745	zzlevbar(ig,:)=zzlevbar(ig-1,:)
[1947]	746	endif
[1672]	747	ENDDO
	748
[2097]	749	else ! if not moyzon
[2366]	750	#endif
[2097]	751
[2098]	752	DO ig=1,ngrid
	753	DO l=1,nlayer
	754	zzlay_eff(ig,l)=gradrad/(g*rad-(pphi(ig,l)+phisfi(ig)))-rad ! reference = geoid
	755	ENDDO
	756	zzlev_eff(ig,1)=phisfi(ig)/g
	757	DO l=2,nlayer
[2097]	758	z1=(pplay(ig,l-1)+pplev(ig,l))/ (pplay(ig,l-1)-pplev(ig,l))
	759	z2=(pplev(ig,l) +pplay(ig,l))/(pplev(ig,l) -pplay(ig,l))
	760	zzlev_eff(ig,l)=(z1zzlay_eff(ig,l-1)+z2zzlay_eff(ig,l))/(z1+z2)
[2098]	761	ENDDO
	762	zzlev_eff(ig,nlayer+1)=zzlay_eff(ig,nlayer)+(zzlay_eff(ig,nlayer)-zzlev_eff(ig,nlayer))
	763	ENDDO
[2097]	764
[2366]	765	#ifndef MESOSCALE
[1672]	766	endif ! moyzon
[2366]	767	#endif
[1795]	768
[1672]	769	! -------------------------------------------------------------------------------------
[1477]	770	! Compute potential temperature
	771	! Note : Potential temperature calculation may not be the same in physiq and dynamic...
	772	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[597]	773	do l=1,nlayer
[787]	774	do ig=1,ngrid
[253]	775	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
[597]	776	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
[1947]	777	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/gzlat(ig,l)
[1542]	778	massarea(ig,l)=mass(ig,l)*cell_area(ig)
[253]	779	enddo
	780	enddo
	781
[1312]	782	! Compute vertical velocity (m/s) from vertical mass flux
[1346]	783	! w = F / (rhoarea) and rho = P/(rT)
[1477]	784	! But first linearly interpolate mass flux to mid-layers
	785	do l=1,nlayer-1
[1914]	786	pw(:,l)=0.5*(flxw(:,l)+flxw(:,l+1))
[1477]	787	enddo
[1914]	788	pw(:,nlayer)=0.5*flxw(:,nlayer) ! since flxw(nlayer+1)=0
[1477]	789	do l=1,nlayer
[1914]	790	pw(:,l)=(pw(:,l)rpt(:,l)) / (pplay(:,l)*cell_area(:))
[1477]	791	enddo
[1194]	792
[1477]	793	!---------------------------------
	794	! II. Compute radiative tendencies
	795	!---------------------------------
[253]	796
	797	if (callrad) then
[526]	798	if( mod(icount-1,iradia).eq.0.or.lastcall) then
[253]	799
[1477]	800	! Compute local stellar zenith angles
	801	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	802	if (tlocked) then
	803	! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb
	804	ztim1=SIN(declin)
	805	ztim2=COS(declin)*COS(zlss)
	806	ztim3=COS(declin)*SIN(zlss)
[253]	807
[1477]	808	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
	809	ztim1,ztim2,ztim3,mu0,fract, flatten)
[253]	810
[1477]	811	elseif (diurnal) then
	812	ztim1=SIN(declin)
	813	ztim2=COS(declin)COS(2.pi*(zday-.5))
	814	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
[253]	815
[1477]	816	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
	817	ztim1,ztim2,ztim3,mu0,fract, flatten)
[3083]	818
[1477]	819	else if(diurnal .eqv. .false.) then
[253]	820
[1542]	821	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad,flatten)
[1161]	822	! WARNING: this function appears not to work in 1D
[253]	823
[1477]	824	endif
[1161]	825
[1477]	826	! Eclipse incoming sunlight (e.g. Saturn ring shadowing).
[1429]	827	if(rings_shadow) then
	828	call call_rings(ngrid, ptime, pday, diurnal)
	829	endif
[1133]	830
[1329]	831
[1477]	832	if (corrk) then
	833
	834	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	835	! II.a Call correlated-k radiative transfer scheme
	836	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1297]	837
[1648]	838	call call_profilgases(nlayer)
	839
[1477]	840	! standard callcorrk
[3083]	841	call callcorrk(ngrid,nlayer,pq,nq,qsurf,zday, &
	842	albedo,albedo_equivalent,emis,mu0,pplev,pplay,zzlev,&
	843	pt,tsurf,fract,dist_star, &
[1482]	844	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, &
	845	fluxsurfabs_sw,fluxtop_lw, &
	846	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, &
[3318]	847	int_dtaui,int_dtauv,zpopthi,zpopthv,zpoptti,zpopttv,&
[3083]	848	lastcall)
[1297]	849
[1482]	850	! Radiative flux from the sky absorbed by the surface (W.m-2).
[1477]	851	GSR=0.0
[1914]	852	fluxrad_sky(:)=emis(:)*fluxsurf_lw(:)+fluxsurfabs_sw(:)
[253]	853
[1477]	854	!if(noradsurf)then ! no lower surface; SW flux just disappears
[1914]	855	! GSR = SUM(fluxsurf_sw(:)*cell_area(:))/totarea
	856	! fluxrad_sky(:)=emis(:)*fluxsurf_lw(:)
[1477]	857	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
	858	!endif
[253]	859
[1477]	860	! Net atmospheric radiative heating rate (K.s-1)
[1914]	861	dtrad(:,:)=zdtsw(:,:)+zdtlw(:,:)
[1498]	862
[1477]	863	elseif(newtonian)then
[1482]	864
	865	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	866	! II.b Call Newtonian cooling scheme
	867	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1477]	868	call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
[253]	869
[1914]	870	zdtsurf(:) = +(pt(:,1)-tsurf(:))/ptimestep
[1477]	871	! e.g. surface becomes proxy for 1st atmospheric layer ?
[253]	872
[1477]	873	else
	874
	875	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	876	! II.c Atmosphere has no radiative effect
	877	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1914]	878	fluxtop_dn(:) = fract(:)mu0(:)Fat1AU/dist_star**2
[1477]	879	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
	880	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
	881	endif
[1914]	882	fluxsurf_sw(:) = fluxtop_dn(:)
[1482]	883	print*,'------------WARNING---WARNING------------' ! by MT2015.
	884	print*,'You are in corrk=false mode, '
[1498]	885	print*,'and the surface albedo is taken equal to the first visible spectral value'
	886
[1914]	887	fluxsurfabs_sw(:) = fluxtop_dn(:)*(1.-albedo(:,1))
	888	fluxrad_sky(:) = fluxsurfabs_sw(:)
	889	fluxtop_lw(:) = emis(:)sigmatsurf(:)**4
[253]	890
[2242]	891	dtrad(:,:)=0.D0 ! no atmospheric radiative heating
[253]	892
[1477]	893	endif ! end of corrk
[253]	894
[1477]	895	endif ! of if(mod(icount-1,iradia).eq.0)
[787]	896
[253]	897
[1477]	898	! Transformation of the radiative tendencies
	899	! ------------------------------------------
[1914]	900	zplanck(:)=tsurf(:)*tsurf(:)
	901	zplanck(:)=emis(:)sigmazplanck(:)*zplanck(:)
	902	fluxrad(:)=fluxrad_sky(:)-zplanck(:)
	903	pdt(:,:)=pdt(:,:)+dtrad(:,:)
[1477]	904
	905	! Test of energy conservation
	906	!----------------------------
[253]	907	if(enertest)then
[1524]	908	call planetwide_sumval(cppmassarea(:,:)zdtsw(:,:)/totarea_planet,dEtotSW)
	909	call planetwide_sumval(cppmassarea(:,:)zdtlw(:,:)/totarea_planet,dEtotLW)
[1542]	910	!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
	911	call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
	912	call planetwide_sumval((fluxsurf_lw(:)emis(:)-zplanck(:))cell_area(:)/totarea_planet,dEtotsLW)
[1524]	913	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
	914	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
	915	if (is_master) then
[1477]	916	print*,'---------------------------------------------------------------'
	917	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
	918	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
	919	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
	920	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
	921	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
	922	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
[1524]	923	endif
[1477]	924	endif ! end of 'enertest'
[253]	925
	926	endif ! of if (callrad)
	927
	928
[1477]	929
	930	! --------------------------------------------
	931	! III. Vertical diffusion (turbulent mixing) :
	932	! --------------------------------------------
[253]	933	if (calldifv) then
[526]	934
[1914]	935	zflubid(:)=fluxrad(:)+fluxgrd(:)
[253]	936
[1477]	937	! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
[1524]	938	if (UseTurbDiff) then
	939
[1647]	940	call turbdiff(ngrid,nlayer,nq, &
[1477]	941	ptimestep,capcal,lwrite, &
	942	pplay,pplev,zzlay,zzlev,z0, &
	943	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
	944	pdt,pdq,zflubid, &
	945	zdudif,zdvdif,zdtdif,zdtsdif, &
[1647]	946	sensibFlux,q2,zdqdif,zdqsdif, &
[1477]	947	taux,tauy,lastcall)
[594]	948
[1524]	949	else
	950
[1914]	951	zdh(:,:)=pdt(:,:)/zpopsk(:,:)
[594]	952
[1647]	953	call vdifc(ngrid,nlayer,nq,zpopsk, &
[1477]	954	ptimestep,capcal,lwrite, &
	955	pplay,pplev,zzlay,zzlev,z0, &
	956	pu,pv,zh,pq,tsurf,emis,qsurf, &
	957	zdh,pdq,zflubid, &
	958	zdudif,zdvdif,zdhdif,zdtsdif, &
[1524]	959	sensibFlux,q2,zdqdif,zdqsdif, &
[1477]	960	taux,tauy,lastcall)
[253]	961
[1914]	962	zdtdif(:,:)=zdhdif(:,:)*zpopsk(:,:) ! for diagnostic only
	963	zdqevap(:,:)=0.
[594]	964
[1477]	965	end if !end of 'UseTurbDiff'
[594]	966
[2291]	967	if (.not. turb_resolved) then
	968	pdv(:,:)=pdv(:,:)+zdvdif(:,:)
	969	pdu(:,:)=pdu(:,:)+zdudif(:,:)
	970	pdt(:,:)=pdt(:,:)+zdtdif(:,:)
	971	endif
[3083]	972
	973	zdtsurf(:)=zdtsurf(:)+zdtsdif(:)
[1477]	974
[253]	975	if (tracer) then
[1914]	976	pdq(:,:,:)=pdq(:,:,:)+ zdqdif(:,:,:)
	977	dqsurf(:,:)=dqsurf(:,:) + zdqsdif(:,:)
[253]	978	end if ! of if (tracer)
	979
[1477]	980
	981	! test energy conservation
[253]	982	!-------------------------
	983	if(enertest)then
[1477]	984
[1524]	985	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
[253]	986	do ig = 1, ngrid
[1524]	987	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
	988	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
[253]	989	enddo
[1477]	990
[1542]	991	call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
[1524]	992	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
[1542]	993	call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
	994	call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
[1477]	995
[1524]	996	if (is_master) then
[1477]	997
	998	if (UseTurbDiff) then
[1524]	999	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
	1000	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
[1477]	1001	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
[1524]	1002	else
	1003	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
	1004	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
	1005	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
	1006	end if
	1007	endif ! end of 'is_master'
[1477]	1008
	1009	! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
	1010	endif ! end of 'enertest'
[253]	1011
[1477]	1012	else ! calldifv
[253]	1013
	1014	if(.not.newtonian)then
	1015
[1914]	1016	zdtsurf(:) = zdtsurf(:) + (fluxrad(:) + fluxgrd(:))/capcal(:)
[253]	1017
	1018	endif
	1019
[1477]	1020	endif ! end of 'calldifv'
[253]	1021
	1022
[1477]	1023	!----------------------------------
	1024	! IV. Dry convective adjustment :
	1025	!----------------------------------
[253]	1026
	1027	if(calladj) then
	1028
[1914]	1029	zdh(:,:) = pdt(:,:)/zpopsk(:,:)
[2242]	1030	zduadj(:,:)=0.D0
	1031	zdvadj(:,:)=0.D0
	1032	zdhadj(:,:)=0.D0
[253]	1033
	1034
[1477]	1035	call convadj(ngrid,nlayer,nq,ptimestep, &
	1036	pplay,pplev,zpopsk, &
	1037	pu,pv,zh,pq, &
	1038	pdu,pdv,zdh,pdq, &
	1039	zduadj,zdvadj,zdhadj, &
	1040	zdqadj)
[253]	1041
[1914]	1042	pdu(:,:) = pdu(:,:) + zduadj(:,:)
	1043	pdv(:,:) = pdv(:,:) + zdvadj(:,:)
	1044	pdt(:,:) = pdt(:,:) + zdhadj(:,:)*zpopsk(:,:)
	1045	zdtadj(:,:) = zdhadj(:,:)*zpopsk(:,:) ! for diagnostic only
[1283]	1046
[253]	1047	if(tracer) then
[1914]	1048	pdq(:,:,:) = pdq(:,:,:) + zdqadj(:,:,:)
[253]	1049	end if
	1050
[1477]	1051	! Test energy conservation
[253]	1052	if(enertest)then
[1524]	1053	call planetwide_sumval(cppmassarea(:,:)zdtadj(:,:)/totarea_planet,dEtot)
[1295]	1054	if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
[253]	1055	endif
	1056
[787]	1057
[1477]	1058	endif ! end of 'calladj'
	1059
[253]	1060
[1477]	1061	!---------------------------------------------
[1647]	1062	! V. Specific parameterizations for tracers
[1477]	1063	!---------------------------------------------
[253]	1064
[1477]	1065	if (tracer) then
[1914]	1066
[2366]	1067	#ifndef MESOSCALE
[3090]	1068	! -------------------
	1069	! V.1 Microphysics
	1070	! -------------------
[2366]	1071
[3090]	1072	! We must call microphysics before chemistry, for condensation !
[1926]	1073	if (callmufi) then
[2242]	1074
	1075	zzlev(:,nlayer+1)=zzlay(:,nlayer)+(zzlay(:,nlayer)-zzlev(:,nlayer)) ! JVO 19 : We assume zzlev isn't reused later on (could be done cleaner)
	1076
[1926]	1077	#ifdef USE_QTEST
[2242]	1078	dtpq(:,:,:) = 0.D0 ! we want tpq to go only through mufi
[1947]	1079	call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,tpq,dtpq,zdqmufi)
[1926]	1080	tpq(:,:,:) = tpq(:,:,:) + zdqmufi(:,:,:)ptimestep ! only manipulation of tpq->ptimestep here
[3090]	1081
	1082	#else
[2109]	1083	call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,zdqmufi)
[1926]	1084	pdq(:,:,:) = pdq(:,:,:) + zdqmufi(:,:,:)
[2242]	1085
[3090]	1086	! Sanity check (way safer to be done here rather than within YAMMS)
[2242]	1087	! Important : the sanity check intentionally include the former processes tendency !
	1088	! NB : Despite this sanity check there might be still some unphysical values going through :
	1089	! - Negatives, but harmless as it will be only for the output files
	1090	! just remove them in post-proc.
	1091	! - Weird unphysical ratio of m0 and m3, ok for now, but take care of them if
	1092	! you want to compute optics from radii.
[3083]	1093	WHERE ((pq(:,:,1)+pdq(:,:,1)ptimestep < 0.D0) .OR. (pq(:,:,2)+pdq(:,:,2)ptimestep < 0.D0))
	1094	pdq(:,:,1) = (epsilon(1.0)-1.D0)*pq(:,:,1)/ptimestep
	1095	pdq(:,:,2) = (epsilon(1.0)-1.D0)*pq(:,:,2)/ptimestep
[2242]	1096	ENDWHERE
[3083]	1097	WHERE ((pq(:,:,3)+pdq(:,:,3)ptimestep < 0.D0) .OR. (pq(:,:,4)+pdq(:,:,4)ptimestep < 0.D0))
	1098	pdq(:,:,3) = (epsilon(1.0)-1.D0)*pq(:,:,3)/ptimestep
	1099	pdq(:,:,4) = (epsilon(1.0)-1.D0)*pq(:,:,4)/ptimestep
[2242]	1100	ENDWHERE
[3083]	1101	IF (callclouds) THEN
	1102	WHERE ((pq(:,:,5)+pdq(:,:,5)ptimestep < 0.D0) .OR. (pq(:,:,6)+pdq(:,:,6)ptimestep < 0.D0))
	1103	pdq(:,:,5) = (epsilon(1.0)-1.D0)*pq(:,:,5)/ptimestep
	1104	pdq(:,:,6) = (epsilon(1.0)-1.D0)*pq(:,:,6)/ptimestep
	1105	ENDWHERE
	1106	DO iq = 1, size(ices_indx)
	1107	! For ices :
	1108	WHERE (pq(:,:,ices_indx(iq))+pdq(:,:,ices_indx(iq))*ptimestep < 0.D0)
	1109	pdq(:,:,ices_indx(iq)) = (epsilon(1.0)-1.D0)*pq(:,:,ices_indx(iq))/ptimestep
	1110	ENDWHERE
	1111	! For gases :
	1112	WHERE (pq(:,:,gazs_indx(iq))+pdq(:,:,gazs_indx(iq))*ptimestep < 0.D0)
	1113	pdq(:,:,gazs_indx(iq)) = (epsilon(1.0)-1.D0)*pq(:,:,gazs_indx(iq))/ptimestep
	1114	ENDWHERE
	1115	ENDDO
	1116	ENDIF
[3497]	1117	! In case there is no clouds, in the troposphere the moments are fixed to evacuate all aerosols
[3340]	1118	IF (.NOT. callclouds) THEN
[3497]	1119	WHERE (pplay(:,:) .gt. 1000.)
	1120	pdq(:,:,1) = 0.
	1121	pdq(:,:,2) = 0.
	1122	ENDWHERE
	1123	WHERE (pplay(:,:) .gt. 5000.)
	1124	pdq(:,:,3) = 0.
	1125	pdq(:,:,4) = 0.
	1126	ENDWHERE
[3340]	1127	ENDIF
[1926]	1128	#endif
	1129
	1130	! Microphysics condensation for 2D fields to sent non-saturated fields to photochem
[3090]	1131	if (callclouds .and. moyzon_ch .and. mod(icount-1,ichim).eq.0) then
	1132	zdqfibar(:,:,:) = 0.D0 ! We work in zonal average -> forget processes other than condensation
	1133	call calmufi(ptimestep,zplevbar,zzlevbar,zplaybar,zzlaybar, &
	1134	gzlat,ztfibar,zqfibar,zdqfibar,zdqmufibar)
	1135	! TODO : Add a sanity check here !
[1926]	1136	endif
[3083]	1137
[3090]	1138	! Condensation heating rate :
	1139	if (callclouds) then
	1140	! Default value -> no condensation [kg/kg_air/s] :
	1141	dmuficond(:,:,:) = 0.D0
	1142	do iq = 1, size(ices_indx)
	1143	dmuficond(:,:,iq) = zdqmufi(:,:,gazs_indx(iq))
	1144	enddo
	1145	call cond_muphy(ngrid,nlayer,pt,dmuficond,zdtlc)
	1146	!pdt(:,:) = pdt(:,:) + zdtlc(:,:)
	1147	endif
	1148	endif ! callmufi
[1926]	1149
	1150	! -----------------
	1151	! V.2. Chemistry
	1152	! -----------------
	1153	! NB : Must be call last ( brings fields back to an equilibrium )
[1672]	1154	if (callchim) then
	1155
[1914]	1156	! o. Convert updated tracers to molar fraction
	1157	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1158	do iq = 1,nkim
[3090]	1159	ychim(:,:,iq) = (pq(:,:,iq+nmicro) + pdq(:,:,iq+nmicro)*ptimestep) / rat_mmol(iq+nmicro)
[1914]	1160	enddo
[1672]	1161
[1926]	1162	! JVO 05/18 : We update zonally averaged fields with condensation
	1163	! as it is compulsory to have correct photochem production. But for other
	1164	! processes ( convadj ... ) we miss them in any case as we work in zonally/diurnal
	1165	! mean -> no fine diurnal/short time evolution, only seasonal evolution only.
[3090]	1166	if (moyzon_ch .and. mod(icount-1,ichim).eq. 0) then
[1926]	1167	do iq = 1,nkim
	1168	ychimbar(:,:,iq) = zqfibar(:,:,iq+nmicro) / rat_mmol(iq+nmicro)
	1169	if ( callclouds ) then
[3083]	1170	ychimbar(:,:,iq) = ychimbar(:,:,iq) + ( zdqmufibar(:,:,iq+nmicro)*ptimestep / rat_mmol(iq+nmicro) )
[1926]	1171	endif
	1172	enddo
	1173	endif
	1174
[1958]	1175	! i. Condensation of the 3D tracers after the transport
	1176	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1966]	1177
[3090]	1178	call calc_ysat(ngrid,nlayer,pplay/100.0,pt,ysat) ! Compute saturation profiles for every grid point (!!p in mbar!!)
[1966]	1179
[2117]	1180	dyccond(:,:,:) = 0.D0 ! Default value -> no condensation
[1966]	1181
[1903]	1182	do iq=1,nkim
[3090]	1183	where (ychim(:,:,iq).gt.ysat(:,:,iq))
	1184	dyccond(:,:,iq+nmicro) = (-ychim(:,:,iq)+ysat(:,:,iq)) / ptimestep
	1185	endwhere
[1672]	1186	enddo
	1187
[3083]	1188	if (callclouds) then
	1189	do iq = 1, size(ices_indx)
	1190	dyccond(:,:,gazs_indx(iq)) = 0.D0 ! Condensation have been calculated in the cloud microphysics
	1191	enddo
	1192	endif
[1795]	1193
[1914]	1194	do iq=1,nkim
[3083]	1195	ychim(:,:,iq) = ychim(:,:,iq) + dyccond(:,:,iq+nmicro)*ptimestep ! update molar ychim for following calchim
[1958]	1196	pdq(:,:,iq+nmicro) = pdq(:,:,iq+nmicro) + dyccond(:,:,iq+nmicro)*rat_mmol(iq+nmicro) ! convert tendencies to mass mixing ratio
[1914]	1197	enddo
	1198
	1199
[1958]	1200	! ii. 2D zonally averaged fields need to condense and evap before photochem
	1201	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3090]	1202	if (moyzon_ch .and. mod(icount-1,ichim).eq. 0) then
[1958]	1203
[3090]	1204	call calc_ysat(ngrid,nlayer,zplaybar/100.0,ztfibar,ysat) ! Compute saturation profiles for every grid point for the zon-ave fields
[1966]	1205
[3090]	1206	dyccondbar(:,:,:) = 0.D0 ! Default value -> no condensation
	1207
	1208	do iq = 1,nkim
	1209	where ( ychimbar(:,:,iq).gt.ysat(:,:,iq) )
	1210	dyccondbar(:,:,iq+nmicro) = ( -ychimbar(:,:,iq)+ysat(:,:,iq) ) / ptimestep
	1211	endwhere
	1212	enddo
[1926]	1213
[3083]	1214	if (callclouds) then
	1215	do iq = 1, size(ices_indx)
	1216	dyccondbar(:,:,gazs_indx(iq)) = 0.D0 ! Condensation have been calculated in the cloud microphysics
	1217	enddo
	1218	endif
[1926]	1219
[3090]	1220	do iq=1,nkim
	1221	ychimbar(:,:,iq) = ychimbar(:,:,iq) + dyccondbar(:,:,iq+nmicro)*ptimestep
	1222	enddo
[1926]	1223
[3090]	1224	endif ! if ( moyzon_ch .and. mod(icount-1,ichim).eq.0 )
[1958]	1225
	1226	! iii. Photochemistry ( must be call after condensation (and evap of 2D) )
	1227	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1795]	1228	if( mod(icount-1,ichim).eq.0. ) then
[3090]	1229	print *, "We enter in the photochemistry ..."
[1672]	1230
[1914]	1231	if (moyzon_ch) then ! 2D zonally averaged chemistry
[1926]	1232	! Here we send zonal average fields ( corrected with cond ) from dynamics to chem. module
[2097]	1233	call calchim(ngrid,ychimbar,declin,ctimestep,ztfibar,zphibar,zphisbar, &
[1958]	1234	zplaybar,zplevbar,zzlaybar,zzlevbar,dycchi)
[3090]	1235
[1904]	1236	else ! 3D chemistry (or 1D run)
[2097]	1237	call calchim(ngrid,ychim,declin,ctimestep,pt,pphi,phisfi, &
	1238	pplay,pplev,zzlay_eff,zzlev_eff,dycchi)
[1914]	1239	endif ! if moyzon
[3090]	1240	endif ! if (mod(icount-1,ichim).eq.0)
[1672]	1241
[1958]	1242	! iv. Surface pseudo-evaporation
	1243	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3318]	1244	! Infinite tank of CH4
[3083]	1245	if (resCH4_inf) then
[3090]	1246	do ig=1,ngrid
	1247	if ((ychim(ig,1,7)+dycchi(ig,1,7)*ptimestep) .lt. botCH4) then ! + dycchi because ychim not yet updated
	1248	dycevapCH4(ig) = ( -ychim(ig,1,7)+botCH4 ) / ptimestep - dycchi(ig,1,7)
	1249	else
	1250	dycevapCH4(ig) = 0.D0
	1251	endif
	1252	enddo
[1958]	1253
[3083]	1254	else
[3318]	1255	! Fill lakes with precipitation :
[3497]	1256	!if (REAL(latitude_deg(ig)) .ge. 70.0) then
	1257	! tankCH4(ig) = 200.0 ! [m]
	1258	!else if (REAL(latitude_deg(ig)) .le. -70.0) then
	1259	! tankCH4(ig) = 50.0 ! [m]
	1260	!else
	1261	! tankCH4(ig) = 0.0 ! [m]
	1262	!endif
	1263	tankCH4(:) = tankCH4(:) + (mmd_ice_prec(:,1) / 422. * ptimestep) ! [m]
[3318]	1264
	1265	! Evaporation of lakes :
[3090]	1266	if (moyzon_ch) then
	1267	tpq_CH4(:,:) = ychimbar(:,:,7) + dycchi(:,:,7)*ptimestep ! + dycchi because ychim not yet updated [mol/mol]
	1268	else
	1269	tpq_CH4(:,:) = ychim(:,:,7) + dycchi(:,:,7)*ptimestep ! + dycchi because ychim not yet updated [mol/mol]
	1270	endif
	1271	call evapCH4(ngrid,nlayer,ptimestep,pplev,zzlay,zzlev,&
	1272	pu,pv,tsurf,tpq_CH4,tankCH4,dycevapCH4,zdtsurfevap)
	1273	zdtsurf(:) = zdtsurf(:) + zdtsurfevap(:)
[3083]	1274	endif
	1275
[3318]	1276
	1277	pdq(:,1,7+nmicro) = pdq(:,1,7+nmicro) + dycevapCH4(:)*rat_mmol(7+nmicro) ! convert tendencies to mass mixing ratio
[1956]	1278
[1958]	1279	! v. Updates and positivity check
	1280	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2242]	1281	zdqchi(:,:,:) = 0.D0 ! -> dycchi is saved but for the nmicro tracers we must update to 0 at each step
[2117]	1282
[1903]	1283	do iq=1,nkim
[1926]	1284	zdqchi(:,:,iq+nmicro) = dycchi(:,:,iq)*rat_mmol(iq+nmicro) ! convert tendencies to mass mixing ratio
[1672]	1285
[1926]	1286	where( (pq(:,:,iq+nmicro) + ( pdq(:,:,iq+nmicro)+zdqchi(:,:,iq+nmicro) )*ptimestep ) .LT. 0.) & ! When using zonal means we set the same tendency
[2117]	1287	zdqchi(:,:,iq+nmicro) = 1.D-30 - pdq(:,:,iq+nmicro) - pq(:,:,iq+nmicro)/ptimestep ! everywhere in longitude -> could lead to negs !
[1926]	1288	enddo
	1289
[1914]	1290	pdq(:,:,:) = pdq(:,:,:) + zdqchi(:,:,:)
[3083]	1291
[1795]	1292	endif ! end of 'callchim'
	1293
[3090]	1294	! END MESOSCALE
[2366]	1295	#endif
	1296
[1477]	1297	! ---------------
[1672]	1298	! V.3 Updates
[1477]	1299	! ---------------
[253]	1300
[1477]	1301	! Updating Atmospheric Mass and Tracers budgets.
[728]	1302	if(mass_redistrib) then
	1303
[1914]	1304	zdmassmr(:,:) = mass(:,:) * zdqevap(:,:)
[863]	1305
	1306	do ig = 1, ngrid
[1914]	1307	zdmassmr_col(ig)=SUM(zdmassmr(ig,:))
[863]	1308	enddo
[728]	1309
[1524]	1310	call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
	1311	call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
	1312	call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
[728]	1313
[1524]	1314	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
[1647]	1315	capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
[1524]	1316	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
	1317	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
	1318
[1914]	1319	pdq(:,:,:) = pdq(:,:,:) + zdqmr(:,:,:)
	1320	dqsurf(:,:) = dqsurf(:,:) + zdqsurfmr(:,:)
	1321	pdt(:,:) = pdt(:,:) + zdtmr(:,:)
	1322	pdu(:,:) = pdu(:,:) + zdumr(:,:)
	1323	pdv(:,:) = pdv(:,:) + zdvmr(:,:)
	1324	pdpsrf(:) = pdpsrf(:) + zdpsrfmr(:)
[3090]	1325	zdtsurf(:) = zdtsurf(:) + zdtsurfmr(:)
[1524]	1326	endif
[728]	1327
[1477]	1328	! -----------------------------
[1672]	1329	! V.4. Surface Tracer Update
[1477]	1330	! -----------------------------
[1297]	1331
[1914]	1332	qsurf(:,:) = qsurf(:,:) + ptimestep*dqsurf(:,:)
[253]	1333
[1477]	1334	! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
	1335	! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
[622]	1336	qsurf_hist(:,:) = qsurf(:,:)
[253]	1337
[1477]	1338	endif! end of if 'tracer'
[253]	1339
	1340
[1477]	1341	!------------------------------------------------
[1647]	1342	! VI. Surface and sub-surface soil temperature
[1477]	1343	!------------------------------------------------
[253]	1344
[1477]	1345
	1346	! Increment surface temperature
[1297]	1347
[1914]	1348	tsurf(:)=tsurf(:)+ptimestep*zdtsurf(:)
[1477]	1349
	1350	! Compute soil temperatures and subsurface heat flux.
[253]	1351	if (callsoil) then
[787]	1352	call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, &
[1477]	1353	ptimestep,tsurf,tsoil,capcal,fluxgrd)
[253]	1354	endif
	1355
[1297]	1356
[1477]	1357	! Test energy conservation
[253]	1358	if(enertest)then
[1542]	1359	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf(:)/totarea_planet,dEtots)
[1524]	1360	if (is_master) print*,'Surface energy change =',dEtots,' W m-2'
[253]	1361	endif
	1362
	1363
[1477]	1364	!---------------------------------------------------
[1647]	1365	! VII. Perform diagnostics and write output files
[1477]	1366	!---------------------------------------------------
[3083]	1367
	1368	! Nudging of zonal wind !
	1369	! ~~~~~~~~~~~~~~~~~~~~~~~
	1370	if (nudging_u) then
	1371	zdundg(:,:) = 0.D0
[3340]	1372	j=1
[3083]	1373	! boucle sur les points de grille :
	1374	do i = 1, ngrid
[3340]	1375	! interpolation linéaire des données dans le fichier lu (sur 49 latitudes)
	1376	do while ((u_ref(j,1).ge.latitude_deg(i)).and.(j.lt.49))
	1377	j=j+1
	1378	enddo
	1379	factlat = (latitude_deg(i)-u_ref(j,1))/(u_ref(j-1,1)-u_ref(j,1))
	1380	! Nudging of the first 23 layers only !!
	1381	! IF CHANGE IN VERTICAL RESOLUTION IN THE FIRST 23 LEVELS, IT DOES NOT WORK !!!
	1382	zundg(1:23) = factlatu_ref(j-1,2:24)+(1-factlat)u_ref(j,2:24)
	1383	zdundg(i,1:23) = (zundg(1:23) - (pu(i,1:23)+pdu(i,1:23)*ptimestep)) / nudging_dt
[3083]	1384	enddo
[3340]	1385
[3083]	1386	pdu(:,:) = pdu(:,:) + zdundg(:,:)
	1387	endif
	1388
[1477]	1389	! Note : For output only: the actual model integration is performed in the dynamics.
[253]	1390
[1477]	1391	! Temperature, zonal and meridional winds.
[1914]	1392	zt(:,:) = pt(:,:) + pdt(:,:)*ptimestep
	1393	zu(:,:) = pu(:,:) + pdu(:,:)*ptimestep
	1394	zv(:,:) = pv(:,:) + pdv(:,:)*ptimestep
[253]	1395
[1477]	1396	! Diagnostic.
[1914]	1397	zdtdyn(:,:) = (pt(:,:)-ztprevious(:,:)) / ptimestep
	1398	ztprevious(:,:) = zt(:,:)
[253]	1399
[1914]	1400	zdudyn(:,:) = (pu(:,:)-zuprevious(:,:)) / ptimestep
	1401	zuprevious(:,:) = zu(:,:)
[1637]	1402
[253]	1403	if(firstcall)then
[2242]	1404	zdtdyn(:,:)=0.D0
	1405	zdudyn(:,:)=0.D0
[253]	1406	endif
	1407
[2002]	1408	! Horizotal wind
	1409	zhorizwind(:,:) = sqrt( zu(:,:)zu(:,:) + zv(:,:)zv(:,:) )
	1410
[1477]	1411	! Dynamical heating diagnostic.
[253]	1412	do ig=1,ngrid
[1637]	1413	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp
[253]	1414	enddo
	1415
[3497]	1416	! [Forcage de la photochimie pour les nuages]
	1417	if (callclouds) then
	1418	do ig = 1, ngrid
	1419	do iq = 1, size(ices_indx)
	1420	! C2H2 :
	1421	!-------
	1422	if(trim(nameOfTracer(gazs_indx(iq))) .eq. "C2H2") then
	1423	pdq(ig,nlayer-3:,gazs_indx(iq)) = (4.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
	1424	endif
	1425	! C2H6 :
	1426	!-------
	1427	if(trim(nameOfTracer(gazs_indx(iq))) .eq. "C2H6") then
	1428	pdq(ig,nlayer-3:,gazs_indx(iq)) = (8.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
	1429	endif
	1430	! HCN :
	1431	!------
	1432	if(trim(nameOfTracer(gazs_indx(iq))) .eq. "HCN") then
	1433	pdq(ig,nlayer-3:,gazs_indx(iq)) = (2.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
	1434	endif
	1435	! AC6H6 :
	1436	!--------
	1437	if(trim(nameOfTracer(gazs_indx(iq))) .eq. "AC6H6") then
	1438	pdq(ig,nlayer-3:,gazs_indx(iq)) = (2.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
	1439	endif
	1440	enddo
[3318]	1441	enddo
[3497]	1442	endif
	1443
[1914]	1444	zq(:,:,:) = pq(:,:,:) + pdq(:,:,:)*ptimestep
[253]	1445
[1477]	1446	! Surface pressure.
[1914]	1447	ps(:) = pplev(:,1) + pdpsrf(:)*ptimestep
[253]	1448
	1449
[1477]	1450	! Surface and soil temperature information
[1542]	1451	call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
[1295]	1452	call planetwide_minval(tsurf(:),Ts2)
	1453	call planetwide_maxval(tsurf(:),Ts3)
[253]	1454	if(callsoil)then
[1542]	1455	TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
[1722]	1456	if (is_master) then
	1457	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
	1458	print*,Ts1,Ts2,Ts3,TsS
	1459	end if
[959]	1460	else
[1722]	1461	if (is_master) then
	1462	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
[1477]	1463	print*,Ts1,Ts2,Ts3
[1524]	1464	endif
[959]	1465	end if
[253]	1466
	1467
[1477]	1468	! Check the energy balance of the simulation during the run
[253]	1469	if(corrk)then
	1470
[1542]	1471	call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
	1472	call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
	1473	call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
	1474	call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
	1475	call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
[787]	1476	do ig=1,ngrid
[253]	1477	if(fluxtop_dn(ig).lt.0.0)then
	1478	print*,'fluxtop_dn has gone crazy'
	1479	print*,'fluxtop_dn=',fluxtop_dn(ig)
	1480	print*,'temp= ',pt(ig,:)
	1481	print*,'pplay= ',pplay(ig,:)
	1482	call abort
	1483	endif
	1484	end do
	1485
[787]	1486	if(ngrid.eq.1)then
[253]	1487	DYN=0.0
	1488	endif
[1524]	1489
	1490	if (is_master) then
[1477]	1491	print*,' ISR ASR OLR GND DYN [W m^-2]'
	1492	print*, ISR,ASR,OLR,GND,DYN
[1524]	1493	endif
[253]	1494
[1295]	1495	if(enertest .and. is_master)then
[651]	1496	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
	1497	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
	1498	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
[253]	1499	endif
	1500
[1295]	1501	if(meanOLR .and. is_master)then
[1216]	1502	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
[253]	1503	! to record global radiative balance
[588]	1504	open(92,file="rad_bal.out",form='formatted',position='append')
[651]	1505	write(92,*) zday,ISR,ASR,OLR
[253]	1506	close(92)
[588]	1507	open(93,file="tem_bal.out",form='formatted',position='append')
[1295]	1508	if(callsoil)then
[1524]	1509	write(93,*) zday,Ts1,Ts2,Ts3,TsS
	1510	else
	1511	write(93,*) zday,Ts1,Ts2,Ts3
	1512	endif
[253]	1513	close(93)
	1514	endif
	1515	endif
	1516
[1477]	1517	endif ! end of 'corrk'
[253]	1518
[651]	1519
[1477]	1520	! Diagnostic to test radiative-convective timescales in code.
[253]	1521	if(testradtimes)then
[588]	1522	open(38,file="tau_phys.out",form='formatted',position='append')
[253]	1523	ig=1
	1524	do l=1,nlayer
	1525	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
	1526	enddo
	1527	close(38)
[726]	1528	print*,'As testradtimes enabled,'
	1529	print*,'exiting physics on first call'
[253]	1530	call abort
	1531	endif
	1532
[1722]	1533	if (is_master) print,'--> Ls =',zls180./pi
[1477]	1534
	1535
	1536	!----------------------------------------------------------------------
[253]	1537	! Writing NetCDF file "RESTARTFI" at the end of the run
[1477]	1538	!----------------------------------------------------------------------
	1539
[253]	1540	! Note: 'restartfi' is stored just before dynamics are stored
	1541	! in 'restart'. Between now and the writting of 'restart',
	1542	! there will have been the itau=itau+1 instruction and
	1543	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	1544	! thus we store for time=time+dtvr
	1545
	1546
	1547
[1477]	1548	if(lastcall) then
	1549	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
[305]	1550
[2291]	1551	#ifndef MESOSCALE
[1477]	1552	if (ngrid.ne.1) then
	1553	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
	1554
	1555	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
	1556	ptimestep,ztime_fin, &
[1789]	1557	tsurf,tsoil,emis,q2,qsurf_hist,tankCH4)
[1477]	1558	endif
[2291]	1559	#endif
	1560
[1896]	1561	endif ! end of 'lastcall'
[253]	1562
[2291]	1563	#ifndef MESOSCALE
[1477]	1564	!-----------------------------------------------------------------------------------------------------
	1565	! OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
	1566	!
	1567	! Note 1 : output with period "ecritphy", set in "run.def"
	1568	!
	1569	! Note 2 : writediagfi can also be called from any other subroutine for any variable,
	1570	! but its preferable to keep all the calls in one place ...
	1571	!-----------------------------------------------------------------------------------------------------
[253]	1572
	1573
[1477]	1574	call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
	1575	call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
	1576	call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
	1577	call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
	1578	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1579	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
	1580	call writediagfi(ngrid,"temp","temperature","K",3,zt)
	1581	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
	1582	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1583	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
	1584	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
	1585	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
[3090]	1586
	1587	! Subsurface temperatures
	1588	!call writediagsoil(ngrid,"tempsoil","temperature soil","K",3,tsoil)
[1477]	1589
[3090]	1590	! Total energy balance diagnostics
[1477]	1591	if(callrad.and.(.not.newtonian))then
[3083]	1592	call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
[1477]	1593	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
	1594	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
	1595	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
	1596
[3083]	1597	call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
	1598	call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
	1599	call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
	1600	call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
	1601	call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
	1602	call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
[1477]	1603
[3083]	1604	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
[1477]	1605	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
	1606	endif ! end of 'callrad'
[1524]	1607
[1477]	1608	if(enertest) then
[1524]	1609	if (calldifv) then
	1610	call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
[1477]	1611	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
[3083]	1612	call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
	1613	call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
	1614	call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
[1524]	1615	endif
[1477]	1616
[1524]	1617	if (corrk) then
	1618	call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
	1619	call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
	1620	endif
[1477]	1621	endif ! end of 'enertest'
[253]	1622
[3083]	1623	! Diagnostics of optical thickness
	1624	! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
	1625	if (diagdtau) then
	1626	do nw=1,L_NSPECTV
[2133]	1627	write(str2,'(i2.2)') nw
[2138]	1628	call writediagfi(ngrid,'dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',1,int_dtauv(:,nlayer:1:-1,nw))
[3083]	1629	enddo
	1630	do nw=1,L_NSPECTI
[2133]	1631	write(str2,'(i2.2)') nw
[2138]	1632	call writediagfi(ngrid,'dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',1,int_dtaui(:,nlayer:1:-1,nw))
[3083]	1633	enddo
	1634	endif
[2133]	1635
[3083]	1636	! Temporary inclusions for winds diagnostics.
	1637	call writediagfi(ngrid,"zdudif","Turbdiff tend. zon. wind","m s-2",3,zdudif)
	1638	call writediagfi(ngrid,"zdudyn","Dyn. tend. zon. wind","m s-2",3,zdudyn)
[1477]	1639
[3083]	1640	! Temporary inclusions for heating diagnostics.
	1641	call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
	1642	call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
	1643	call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
	1644	call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
[1477]	1645
[3083]	1646	! For Debugging.
	1647	call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
[253]	1648
[1477]	1649	! Output tracers.
	1650	if (tracer) then
[1808]	1651
[1926]	1652	if (callmufi) then
	1653	! Microphysical tracers are expressed in unit/m3.
[1947]	1654	! convert X.kg-1 --> X.m-3 (whereas for optics was -> X.m-2)
	1655	i2e(:,:) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / gzlat(:,1:nlayer) /(zzlev(:,2:nlayer+1)-zzlev(:,1:nlayer))
[1926]	1656
[1897]	1657	#ifdef USE_QTEST
	1658	! Microphysical tracers passed through dyn+phys(except mufi)
[1947]	1659	call writediagfi(ngrid,"mu_m0as_dp","Dynphys only spherical mode 0th order moment",'m-3',3,zq(:,:,micro_indx(1))*i2e)
	1660	call writediagfi(ngrid,"mu_m3as_dp","Dynphys only spherical mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(2))*i2e)
	1661	call writediagfi(ngrid,"mu_m0af_dp","Dynphys only fractal mode 0th order moment",'m-3',3,zq(:,:,micro_indx(3))*i2e)
	1662	call writediagfi(ngrid,"mu_m3af_dp","Dynphys only fractal mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(4))*i2e)
[1897]	1663	! Microphysical tracers passed through mufi only
[1947]	1664	call writediagfi(ngrid,"mu_m0as_mo","Mufi only spherical mode 0th order moment",'m-3',3,tpq(:,:,micro_indx(1))*i2e)
	1665	call writediagfi(ngrid,"mu_m3as_mo","Mufi only spherical mode 3rd order moment",'m3/m3',3,tpq(:,:,micro_indx(2))*i2e)
	1666	call writediagfi(ngrid,"mu_m0af_mo","Mufi only fractal mode 0th order moment",'m-3',3,tpq(:,:,micro_indx(3))*i2e)
	1667	call writediagfi(ngrid,"mu_m3af_mo","Mufi only fractal mode 3rd order moment",'m3/m3',3,tpq(:,:,micro_indx(4))*i2e)
[1897]	1668	#else
[1947]	1669	call writediagfi(ngrid,"mu_m0as","Spherical mode 0th order moment",'m-3',3,zq(:,:,micro_indx(1))*i2e)
	1670	call writediagfi(ngrid,"mu_m3as","Spherical mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(2))*i2e)
	1671	call writediagfi(ngrid,"mu_m0af","Fractal mode 0th order moment",'m-3',3,zq(:,:,micro_indx(3))*i2e)
	1672	call writediagfi(ngrid,"mu_m3af","Fractal mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(4))*i2e)
[1897]	1673	#endif
[1926]	1674
	1675	! Microphysical diagnostics
	1676	call writediagfi(ngrid,"mmd_aer_prec","Total aerosols precipitations",'m',2,mmd_aer_prec)
	1677	call writediagfi(ngrid,"mmd_aer_s_flux","Spherical aerosols sedimentation flux",'kg.m-2.s-1',3,mmd_aer_s_flux)
	1678	call writediagfi(ngrid,"mmd_aer_f_flux","Fractal aerosols sedimentation flux",'kg.m-2.s-1',3,mmd_aer_f_flux)
	1679	call writediagfi(ngrid,"mmd_rc_sph","Spherical mode caracteristic radius",'m',3,mmd_rc_sph)
	1680	call writediagfi(ngrid,"mmd_rc_fra","Fractal mode caracteristic radius",'m',3,mmd_rc_fra)
[1795]	1681	endif ! end of 'callmufi'
[253]	1682
[1926]	1683	! Chemical tracers
	1684	if (callchim) then
	1685	do iq=1,nkim
	1686	call writediagfi(ngrid,cnames(iq),cnames(iq),'mol/mol',3,zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro))
[3083]	1687	call writediagfi(ngrid,'dqcond_'//cnames(iq),'dqcond_'//cnames(iq),'mol/mol/s',3,dyccond(:,:,iq+nmicro))
[1926]	1688	enddo
[1960]	1689	call writediagfi(ngrid,"evapCH4","Surface CH4 pseudo-evaporation rate",'mol/mol/s',2,dycevapCH4)
[1926]	1690	endif
[1808]	1691
[1477]	1692	endif ! end of 'tracer'
[253]	1693
[3083]	1694	#ifdef CPP_XIOS
	1695	!-----------------------------------------------------------------------------------------------------
[1622]	1696	! XIOS outputs
[3083]	1697	!-----------------------------------------------------------------------------------------------------
[1622]	1698	! Send fields to XIOS: (NB these fields must also be defined as
	1699	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
[3090]	1700
[3083]	1701	!--------------------------------------------------------
	1702	! General diagnostics :
	1703	!--------------------------------------------------------
[1915]	1704	CALL send_xios_field("ls",zls*180./pi)
	1705	CALL send_xios_field("lss",zlss*180./pi)
	1706	CALL send_xios_field("RA",right_ascen*180./pi)
	1707	CALL send_xios_field("Declin",declin*180./pi)
[1626]	1708
[3083]	1709	! Atmosphere (3D) :
	1710	CALL send_xios_field("temp",zt)
	1711	CALL send_xios_field("teta",zh)
	1712	CALL send_xios_field("p",pplay)
	1713	CALL send_xios_field("u",zu)
	1714	CALL send_xios_field("v",zv)
	1715	CALL send_xios_field("w",pw)
	1716
[1915]	1717	CALL send_xios_field("area",cell_area)
	1718	CALL send_xios_field("pphi",pphi)
[3083]	1719
	1720	! Surface (2D) :
[1622]	1721	CALL send_xios_field("ps",ps)
[1915]	1722	CALL send_xios_field("tsurf",tsurf)
[3083]	1723	CALL send_xios_field("pphis",phisfi)
[1622]	1724
[3083]	1725	! Total energy balance diagnostics (2D) :
	1726	IF (callrad.and.(.not.newtonian)) THEN
	1727	CALL send_xios_field("ISR_TOA",fluxtop_dn)
	1728	CALL send_xios_field("OLR_TOA",fluxtop_lw)
	1729	ENDIF
[1957]	1730
[3083]	1731	!--------------------------------------------------------
	1732	! Winds trends :
	1733	!--------------------------------------------------------
	1734	! Atmosphere (3D) :
	1735	! du_tot = zdudyn + pdu
	1736	CALL send_xios_field("dudyn",zdudyn)
	1737	CALL send_xios_field("pdu",pdu)
[1915]	1738
[3083]	1739	! pdu = zdudif + zduadj + zdundg
[1915]	1740	CALL send_xios_field("dudif",zdudif)
[3083]	1741	CALL send_xios_field("duadj",zduadj)
[3090]	1742	IF (nudging_u) THEN
	1743	CALL send_xios_field("dundg",zdundg)
	1744	ENDIF
[3083]	1745
	1746	! zhorizwind = sqrt(uu + vv)
[2002]	1747	CALL send_xios_field("horizwind",zhorizwind)
	1748
[3083]	1749
	1750	!--------------------------------------------------------
	1751	! Heating trends :
	1752	!--------------------------------------------------------
	1753	! Atmosphere (3D) :
	1754	! dt_tot = dtdyn + pdt
	1755	CALL send_xios_field("dtdyn",zdtdyn)
	1756	CALL send_xios_field("pdt",pdt)
	1757
	1758	! pdt = dtrad + zdtdif + dtadj + zdtlc
[1915]	1759	CALL send_xios_field("dtrad",dtrad)
[2097]	1760	CALL send_xios_field("dtdif",zdtdif)
[3083]	1761	CALL send_xios_field("dtadj",zdtadj(:,:))
[3090]	1762	IF (callclouds) THEN
	1763	CALL send_xios_field("dtlc",zdtlc)
	1764	ENDIF
[1926]	1765
[3083]	1766	! dtrad = zdtsw + zdtlw
	1767	CALL send_xios_field("dtsw",zdtsw)
	1768	CALL send_xios_field("dtlw",zdtlw)
	1769
	1770	! Surface (2D) :
	1771	IF(enertest) THEN
	1772	IF (calldifv) THEN
	1773	CALL send_xios_field("sensibFlux",sensibFlux)
	1774	ENDIF
	1775	ENDIF
	1776	CALL send_xios_field("fluxsurf_lw",fluxsurf_lw(:))
	1777	CALL send_xios_field("fluxsurfabs_sw",fluxsurfabs_sw(:))
	1778	CALL send_xios_field("emis",emis(:))
[1943]	1779
[3083]	1780	! dtsurf = dtsdif + dtsurfevap
	1781	CALL send_xios_field("dtsurf",zdtsurf(:))
	1782	CALL send_xios_field("dtsdif",zdtsdif(:))
	1783	CALL send_xios_field("dtsurfevap",zdtsurfevap(:))
	1784
[1947]	1785
[3083]	1786	!--------------------------------------------------------
	1787	! Optical diagnostics :
	1788	!--------------------------------------------------------
[3318]	1789	! Haze opacity :
[3497]	1790	CALL send_xios_field('ttauhv_08',zpopthv(:,:,8,2)) ! 08 --> 1.983 um
	1791	CALL send_xios_field('ttauhv_15',zpopthv(:,:,15,2)) ! 15 --> 1.000 um
[3318]	1792	CALL send_xios_field('ttauhv_20',zpopthv(:,:,20,2)) ! 20 --> 0.671 um
	1793	CALL send_xios_field('ttauhv_23',zpopthv(:,:,23,2)) ! 23 --> 0.346 um
	1794	CALL send_xios_field('ttauhi_02',zpopthi(:,:,2,2)) ! 02 --> 175.3 um
	1795	CALL send_xios_field('ttauhi_17',zpopthi(:,:,17,2)) ! 17 --> 11.00 um
	1796	CALL send_xios_field('ttauhi_23',zpopthi(:,:,23,2)) ! 23 --> 4.849 um
	1797	! Haze extinction :
[3497]	1798	CALL send_xios_field('kkhv_08',zpopthv(:,:,8,3))
	1799	CALL send_xios_field('kkhv_15',zpopthv(:,:,15,3))
[3318]	1800	CALL send_xios_field('kkhv_20',zpopthv(:,:,20,3))
	1801	CALL send_xios_field('kkhv_23',zpopthv(:,:,23,3))
	1802	CALL send_xios_field('kkhi_02',zpopthi(:,:,2,3))
	1803	CALL send_xios_field('kkhi_17',zpopthi(:,:,17,3))
	1804	CALL send_xios_field('kkhi_23',zpopthi(:,:,23,3))
	1805	! Haze single scattering albedo :
[3497]	1806	CALL send_xios_field('wwhv_08',zpopthv(:,:,8,4))
	1807	CALL send_xios_field('wwhv_15',zpopthv(:,:,15,4))
[3318]	1808	CALL send_xios_field('wwhv_20',zpopthv(:,:,20,4))
	1809	CALL send_xios_field('wwhv_23',zpopthv(:,:,23,4))
	1810	CALL send_xios_field('wwhi_02',zpopthi(:,:,2,4))
	1811	CALL send_xios_field('wwhi_17',zpopthi(:,:,17,4))
	1812	CALL send_xios_field('wwhi_23',zpopthi(:,:,23,4))
	1813	! Haze asymmetry parameter :
[3497]	1814	CALL send_xios_field('gghv_08',zpopthv(:,:,8,5))
	1815	CALL send_xios_field('gghv_15',zpopthv(:,:,15,5))
[3318]	1816	CALL send_xios_field('gghv_20',zpopthv(:,:,20,5))
	1817	CALL send_xios_field('gghv_23',zpopthv(:,:,23,5))
	1818	CALL send_xios_field('gghi_02',zpopthi(:,:,2,5))
	1819	CALL send_xios_field('gghi_17',zpopthi(:,:,17,5))
	1820	CALL send_xios_field('gghi_23',zpopthi(:,:,23,5))
	1821
	1822	! Diagnostics for haze and clouds :
	1823	IF (callclouds) THEN
	1824	! Opacity :
[3497]	1825	CALL send_xios_field('ttauv_08',zpopttv(:,:,8,2)) ! 08 --> 1.983 um
	1826	CALL send_xios_field('ttauv_15',zpopttv(:,:,15,2)) ! 15 --> 1.000 um
[3318]	1827	CALL send_xios_field('ttauv_20',zpopttv(:,:,20,2)) ! 20 --> 0.671 um
	1828	CALL send_xios_field('ttauv_23',zpopttv(:,:,23,2)) ! 23 --> 0.346 um
	1829	CALL send_xios_field('ttaui_02',zpoptti(:,:,2,2)) ! 02 --> 175.3 um
	1830	CALL send_xios_field('ttaui_17',zpoptti(:,:,17,2)) ! 17 --> 11.00 um
	1831	CALL send_xios_field('ttaui_23',zpoptti(:,:,23,2)) ! 23 --> 4.849 um
	1832	! Extinction :
[3497]	1833	CALL send_xios_field('kkv_08',zpopttv(:,:,8,3))
	1834	CALL send_xios_field('kkv_15',zpopttv(:,:,15,3))
[3318]	1835	CALL send_xios_field('kkv_20',zpopttv(:,:,20,3))
	1836	CALL send_xios_field('kkv_23',zpopttv(:,:,23,3))
	1837	CALL send_xios_field('kki_02',zpoptti(:,:,2,3))
	1838	CALL send_xios_field('kki_17',zpoptti(:,:,17,3))
	1839	CALL send_xios_field('kki_23',zpoptti(:,:,23,3))
	1840	! Single scattering albedo :
[3497]	1841	CALL send_xios_field('wwv_08',zpopttv(:,:,8,4))
	1842	CALL send_xios_field('wwv_15',zpopttv(:,:,15,4))
[3318]	1843	CALL send_xios_field('wwv_20',zpopttv(:,:,20,4))
	1844	CALL send_xios_field('wwv_23',zpopttv(:,:,23,4))
	1845	CALL send_xios_field('wwi_02',zpoptti(:,:,2,4))
	1846	CALL send_xios_field('wwi_17',zpoptti(:,:,17,4))
	1847	CALL send_xios_field('wwi_23',zpoptti(:,:,23,4))
	1848	! Asymmetry parameter :
[3497]	1849	CALL send_xios_field('ggv_08',zpopttv(:,:,8,5))
	1850	CALL send_xios_field('ggv_15',zpopttv(:,:,15,5))
[3318]	1851	CALL send_xios_field('ggv_20',zpopttv(:,:,20,5))
	1852	CALL send_xios_field('ggv_23',zpopttv(:,:,23,5))
	1853	CALL send_xios_field('ggi_02',zpoptti(:,:,2,5))
	1854	CALL send_xios_field('ggi_17',zpoptti(:,:,17,5))
	1855	CALL send_xios_field('ggi_23',zpoptti(:,:,23,5))
	1856	! DRAYAER, TAUGAS, DCONT :
	1857	CALL send_xios_field('drayaerv_20',zpopttv(:,:,20,6)) ! 20 --> 0.671um
	1858	CALL send_xios_field('taugasv_20',zpopttv(:,:,20,7))
	1859	CALL send_xios_field('dcontv_20',zpopttv(:,:,20,8))
	1860	CALL send_xios_field('drayaeri_17',zpoptti(:,:,17,6)) ! 17 --> 11.00um
	1861	CALL send_xios_field('taugasi_17',zpoptti(:,:,17,7))
	1862	CALL send_xios_field('dconti_17',zpoptti(:,:,17,8))
	1863	ENDIF
	1864
	1865	! Diagnostics for haze and clouds (4D) :
[3090]	1866	CALL send_xios_field('dtauhi',zpopthi(:,:,:,1))
	1867	CALL send_xios_field('tauhi',zpopthi(:,:,:,2))
	1868	CALL send_xios_field('khi',zpopthi(:,:,:,3))
	1869	CALL send_xios_field('whi',zpopthi(:,:,:,4))
	1870	CALL send_xios_field('ghi',zpopthi(:,:,:,5))
	1871	CALL send_xios_field('dtauhv',zpopthv(:,:,:,1))
	1872	CALL send_xios_field('tauhv',zpopthv(:,:,:,2))
	1873	CALL send_xios_field('khv',zpopthv(:,:,:,3))
	1874	CALL send_xios_field('whv',zpopthv(:,:,:,4))
	1875	CALL send_xios_field('ghv',zpopthv(:,:,:,5))
	1876	IF (callclouds) THEN
[3318]	1877	CALL send_xios_field('dtaui',zpoptti(:,:,:,1))
	1878	CALL send_xios_field('taui',zpoptti(:,:,:,2))
	1879	CALL send_xios_field('ki',zpoptti(:,:,:,3))
	1880	CALL send_xios_field('wi',zpoptti(:,:,:,4))
	1881	CALL send_xios_field('gi',zpoptti(:,:,:,5))
	1882	CALL send_xios_field('dtauv',zpopttv(:,:,:,1))
	1883	CALL send_xios_field('tauv',zpopttv(:,:,:,2))
	1884	CALL send_xios_field('kv',zpopttv(:,:,:,3))
	1885	CALL send_xios_field('wv',zpopttv(:,:,:,4))
	1886	CALL send_xios_field('gv',zpopttv(:,:,:,5))
	1887	ENDIF
[1926]	1888
[3083]	1889	!--------------------------------------------------------
	1890	! Microphysical tracers :
	1891	!--------------------------------------------------------
	1892	IF (callmufi) THEN
	1893	! Atmosphere (3D) :
	1894	! Moments M0 and M3 :
	1895	CALL send_xios_field("mu_m0as",zq(:,:,micro_indx(1))*i2e)
	1896	CALL send_xios_field("mu_m3as",zq(:,:,micro_indx(2))*i2e)
	1897	CALL send_xios_field("mu_m0af",zq(:,:,micro_indx(3))*i2e)
	1898	CALL send_xios_field("mu_m3af",zq(:,:,micro_indx(4))*i2e)
	1899	IF (callclouds) THEN
	1900	CALL send_xios_field("mu_m0n",zq(:,:,micro_indx(5))*i2e)
	1901	CALL send_xios_field("mu_m3n",zq(:,:,micro_indx(6))*i2e)
	1902	DO iq = 1, size(ices_indx)
	1903	CALL send_xios_field(TRIM(nameOfTracer(ices_indx(iq))),zq(:,:,ices_indx(iq))*i2e)
	1904	ENDDO
	1905	ENDIF
	1906
	1907	! Microphysical diagnostics :
	1908	CALL send_xios_field("rc_sph",mmd_rc_sph(:,:))
	1909	CALL send_xios_field("rc_fra",mmd_rc_fra(:,:))
	1910	CALL send_xios_field("vsed_aers",mmd_aer_s_w(:,:))
	1911	CALL send_xios_field("vsed_aerf",mmd_aer_f_w(:,:))
	1912	CALL send_xios_field("flux_aers",mmd_aer_s_flux(:,:))
	1913	CALL send_xios_field("flux_aerf",mmd_aer_f_flux(:,:))
	1914	IF (callclouds) THEN
	1915	CALL send_xios_field("rc_cld",mmd_rc_cld(:,:))
	1916	CALL send_xios_field("vsed_ccn",mmd_ccn_w(:,:))
	1917	CALL send_xios_field("flux_ccn",mmd_ccn_flux(:,:))
	1918	DO iq = 1, size(ices_indx)
[3090]	1919	CALL send_xios_field('flux_i'//TRIM(nameOfTracer(gazs_indx(iq))),mmd_ice_fluxes(:,:,iq))
[3083]	1920	CALL send_xios_field(TRIM(nameOfTracer(gazs_indx(iq)))//'_sat',mmd_gazs_sat(:,:,iq))
	1921	ENDDO
	1922	ENDIF
	1923
	1924	! Surface (2D) :
	1925	CALL send_xios_field("aer_prec",mmd_aer_prec(:))
[3090]	1926	IF (callclouds) THEN
[3340]	1927	CALL send_xios_field("ccn_prec",mmd_ccn_prec(:))
[3090]	1928	DO iq = 1, size(ices_indx)
	1929	CALL send_xios_field('i'//TRIM(nameOfTracer(gazs_indx(iq)))//'_prec',mmd_ice_prec(:,iq))
	1930	ENDDO
	1931	ENDIF
[3083]	1932	ENDIF ! of 'if callmufi'
	1933
	1934	!--------------------------------------------------------
	1935	! Chemical tracers :
	1936	!--------------------------------------------------------
	1937	IF (callchim) THEN
[3318]	1938	! Surface (2D) :
	1939	CALL send_xios_field("evapCH4",dycevapCH4(:)) ! Pseudo-evaporation flux (mol/mol/s)
	1940	CALL send_xios_field("tankCH4",tankCH4(:)) ! CH4 tank at the surface (m)
	1941
[3083]	1942	! Atmosphere (3D) :
	1943	! Chemical species :
	1944	DO iq = 1, nkim
[3090]	1945	! If no cloud : gzs_indx uninitialized
[3083]	1946	CALL send_xios_field(trim(cnames(iq)),zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro)) ! kg/kg -> mol/mol
	1947	ENDDO
	1948
	1949	! Condensation tendencies from microphysics (mol/mol/s) :
[3090]	1950	IF (callclouds) THEN
	1951	DO iq = 1, size(ices_indx)
	1952	CALL send_xios_field('dmuficond_'//trim(nameOfTracer(gazs_indx(iq))),dmuficond(:,:,iq)/rat_mmol(gazs_indx(iq))) ! kg/kg/s -> mol/mol/s
	1953	ENDDO
	1954	ENDIF
[3083]	1955
	1956	! Condensation tendencies (mol/mol/s) :
	1957	CALL send_xios_field("dqcond_CH4",dyccond(:,:,7+nmicro))
	1958	CALL send_xios_field("dqcond_C2H2",dyccond(:,:,10+nmicro))
	1959	CALL send_xios_field("dqcond_C2H4",dyccond(:,:,12+nmicro))
	1960	CALL send_xios_field("dqcond_C2H6",dyccond(:,:,14+nmicro))
	1961	CALL send_xios_field("dqcond_C3H6",dyccond(:,:,17+nmicro))
	1962	CALL send_xios_field("dqcond_C4H4",dyccond(:,:,21+nmicro))
	1963	CALL send_xios_field("dqcond_CH3CCH",dyccond(:,:,23+nmicro))
	1964	CALL send_xios_field("dqcond_C3H8",dyccond(:,:,24+nmicro))
	1965	CALL send_xios_field("dqcond_C4H2",dyccond(:,:,25+nmicro))
	1966	CALL send_xios_field("dqcond_C4H6",dyccond(:,:,26+nmicro))
	1967	CALL send_xios_field("dqcond_C4H10",dyccond(:,:,27+nmicro))
	1968	CALL send_xios_field("dqcond_AC6H6",dyccond(:,:,28+nmicro))
	1969	CALL send_xios_field("dqcond_HCN",dyccond(:,:,35+nmicro))
	1970	CALL send_xios_field("dqcond_CH3CN",dyccond(:,:,39+nmicro))
	1971	CALL send_xios_field("dqcond_HC3N",dyccond(:,:,41+nmicro))
	1972	CALL send_xios_field("dqcond_NCCN",dyccond(:,:,42+nmicro))
	1973	CALL send_xios_field("dqcond_C4N2",dyccond(:,:,43+nmicro))
	1974
	1975	! Upper atmosphere chemistry variables (3D) :
	1976	DO iq = 1, nkim
	1977	CALL send_xios_field(trim(cnames(iq))//"_up",ykim_up(iq,:,:)) ! mol/mol
	1978	ENDDO
	1979
	1980	! Total atmosphere chemistry variables (3D) :
	1981	! Append fields in ykim_tot for output on the total vertical grid (0->1300km)
	1982	DO iq = 1, nkim
	1983	! GCM levels
	1984	DO l = 1, nlayer
	1985	ykim_tot(iq,:,l) = zq(:,l,iq+nmicro)/rat_mmol(iq+nmicro)
	1986	ENDDO
	1987	! Upper levels
	1988	DO l = 1, nlaykim_up
	1989	ykim_tot(iq,:,nlayer+l) = ykim_up(iq,:,l)
	1990	ENDDO
	1991	CALL send_xios_field(trim(cnames(iq))//"_tot",ykim_tot(iq,:,:)) ! mol/mol
	1992	ENDDO
	1993	ENDIF ! of 'if callchim'
	1994
	1995
[1896]	1996	if (lastcall.and.is_omp_master) then
	1997	write(,) "physiq: call xios_context_finalize"
	1998	call xios_context_finalize
	1999	endif
[1622]	2000	#endif
[2291]	2001	#else
	2002	!MESOSCALE outputs
	2003	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
	2004	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
	2005	comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid)
	2006	comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid)
	2007	comm_FLUXTOP_LW(1:ngrid)=fluxtop_lw(1:ngrid)
	2008	comm_FLUXSURF_SW(1:ngrid)=fluxsurf_sw(1:ngrid)
	2009	comm_FLUXSURF_LW(1:ngrid)=fluxsurf_lw(1:ngrid)
	2010	comm_FLXGRD(1:ngrid)=fluxgrd(1:ngrid)
[2742]	2011	sensibFlux(1:ngrid) = zflubid(1:ngrid) - capcal(1:ngrid)*zdtsdif(1:ngrid)
[2291]	2012	#endif
[1622]	2013
[253]	2014	icount=icount+1
	2015
[751]	2016	end subroutine physiq
[1549]	2017
	2018	end module physiq_mod

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

Download in other formats:

Original Format