Context Navigation

physiq_mod.F90 @ 3300

Last change on this file since 3300 was 3299, checked in by emillour, 8 months ago

Generic PCM:

add some auxiliary functions in module tracer_h.F90: is_known_tracer(tname) returns ".true." if "tname" is a known tracer tracer_index(tname) returns the tracer index of tracer "tname"
add the possibility of "volcanic outgasing" as point sources. controled by a "callvolcano" flag (default .false.) details about source location and injection details (frequency, length) need be provided in a "volcano.def" ASCII file read at runtime. (see routine "read_volcano" in volcano.F90)

AB+EM

Property svn:executable set to *

File size: 122.3 KB

Rev	Line
[1549]	1	module physiq_mod
[3232]	2
[1549]	3	implicit none
[3232]	4
[1549]	5	contains
[3232]	6
[751]	7	subroutine physiq(ngrid,nlayer,nq, &
[253]	8	firstcall,lastcall, &
	9	pday,ptime,ptimestep, &
	10	pplev,pplay,pphi, &
	11	pu,pv,pt,pq, &
[1312]	12	flxw, &
[1576]	13	pdu,pdv,pdt,pdq,pdpsrf)
[2663]	14
[3100]	15	!!
[3232]	16	use write_field_phy, only: Writefield_phy
[3100]	17	!!
[3232]	18	use ioipsl_getin_p_mod, only: getin_p
[2736]	19	use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind, corrkdir, banddir
[2060]	20	use watercommon_h, only : RLVTT, Psat_water,epsi,su_watercycle, RV, T_h2o_ice_liq
[2728]	21	use generic_cloud_common_h, only : epsi_generic, Psat_generic
[2129]	22	use thermcell_mod, only: init_thermcell_mod
[3233]	23	use gases_h, only: gnom, gfrac, ngasmx
[2736]	24	use radcommon_h, only: sigma, glat, grav, BWNV, WNOI, DWNI, DWNV, WNOV
[2972]	25	use suaer_corrk_mod, only: suaer_corrk
[1308]	26	use radii_mod, only: h2o_reffrad, co2_reffrad
[2898]	27	use aerosol_mod, only: iniaerosol, iaero_co2, iaero_h2o
[1482]	28	use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, &
[2784]	29	dryness
[1216]	30	use comdiurn_h, only: coslat, sinlat, coslon, sinlon
[1327]	31	use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
[1216]	32	use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
[1543]	33	use geometry_mod, only: latitude, longitude, cell_area
[1542]	34	USE comgeomfi_h, only: totarea, totarea_planet
[1216]	35	USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
	36	alpha_lift, alpha_devil, qextrhor, &
	37	igcm_h2o_ice, igcm_h2o_vap, igcm_dustbin, &
[2728]	38	igcm_co2_ice, nesp, is_chim, is_condensable,constants_epsi_generic
[1525]	39	use time_phylmdz_mod, only: ecritphy, iphysiq, nday
[1669]	40	use phyetat0_mod, only: phyetat0
[2958]	41	use wstats_mod, only: callstats, wstats, mkstats
[1216]	42	use phyredem, only: physdem0, physdem1
[3100]	43	!! use slab_ice_h, only: capcalocean, capcalseaice,capcalsno, &
	44	!! noceanmx
	45	!! use ocean_slab_mod, only :ocean_slab_init, ocean_slab_ice, &
	46	!! capcalocean, capcalseaice,capcalsno, nslay, &
	47	!! ocean_slab_final
	48	!! use surf_heat_transp_mod,only :init_masquv
	49	use ocean_slab_mod, only :ocean_slab_init, ocean_slab_noice, ocean_slab_ice, &
	50	ocean_slab_frac, ocean_slab_get_vars, &
	51	capcalocean, capcalseaice,capcalsno, nslay, knon, &
	52	ocean_slab_final
[1295]	53	use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
	54	use mod_phys_lmdz_para, only : is_master
[1308]	55	use planete_mod, only: apoastr, periastr, year_day, peri_day, &
	56	obliquit, nres, z0
[1524]	57	use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
	58	use time_phylmdz_mod, only: daysec
[2446]	59	use callkeys_mod, only: albedo_spectral_mode, calladj, calldifv, &
[2603]	60	calllott_nonoro, callrad, callsoil, nosurf, &
[2446]	61	calltherm, CLFvarying, co2cond, corrk, diagdtau, &
	62	diurnal, enertest, fat1au, flatten, j2, &
[3204]	63	hydrology, intheat, iradia, kastprof, &
[2446]	64	lwrite, mass_redistrib, massplanet, meanOLR, &
	65	nearco2cond, newtonian, noseason_day, oblate, &
	66	ok_slab_ocean, photochem, rings_shadow, rmean, &
[2700]	67	season, sedimentation,generic_condensation, &
[3299]	68	specOLR, callthermos, callvolcano, &
[2446]	69	startphy_file, testradtimes, tlocked, &
	70	tracer, UseTurbDiff, water, watercond, &
[3233]	71	waterrain, generic_rain, global1d, szangle, &
[3277]	72	moistadjustment, moistadjustment_generic, varspec, varspec_data, nvarlayer
[2721]	73	use generic_tracer_index_mod, only: generic_tracer_index
[2299]	74	use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
[2663]	75	use check_fields_mod, only: check_physics_fields
[2595]	76	use conc_mod, only: rnew, cpnew, ini_conc_mod
[1836]	77	use phys_state_var_mod
[2032]	78	use callcorrk_mod, only: callcorrk
[3299]	79	use conduction_mod, only: conduction
	80	use volcano_mod, only: volcano
[3233]	81	use pindex_mod, only: pindex
[2427]	82	use vdifc_mod, only: vdifc
	83	use turbdiff_mod, only: turbdiff
[1836]	84	use turb_mod, only : q2,sensibFlux,turb_resolved
[2428]	85	use mass_redistribution_mod, only: mass_redistribution
[3232]	86	use condensation_generic_mod, only: condensation_generic
[2736]	87	use datafile_mod, only: datadir
[1836]	88	#ifndef MESOSCALE
[1622]	89	use vertical_layers_mod, only: presnivs, pseudoalt
[1682]	90	use mod_phys_lmdz_omp_data, ONLY: is_omp_master
[1836]	91	#else
[2865]	92	use comm_wrf, only : comm_HR_SW, comm_HR_LW, &
	93	comm_CLOUDFRAC,comm_TOTCLOUDFRAC, comm_RH, &
	94	comm_SURFRAIN,comm_REEVAP,comm_HR_DYN, &
	95	comm_DQICE,comm_DQVAP,comm_ALBEQ, &
	96	comm_FLUXTOP_DN,comm_FLUXABS_SW, &
	97	comm_FLUXTOP_LW,comm_FLUXSURF_SW, &
	98	comm_FLUXSURF_LW,comm_FLXGRD, &
[2871]	99	comm_DTRAIN,comm_DTLSC,comm_H2OICE_REFF, &
	100	comm_LATENT_HF
	101
[1836]	102	#endif
	103
[3232]	104	#ifdef CPP_XIOS
[1622]	105	use xios_output_mod, only: initialize_xios_output, &
	106	update_xios_timestep, &
	107	send_xios_field
[1682]	108	use wxios, only: wxios_context_init, xios_context_finalize
[1623]	109	#endif
[1836]	110
[253]	111	implicit none
	112
	113
	114	!==================================================================
[3232]	115	!
[253]	116	! Purpose
	117	! -------
	118	! Central subroutine for all the physics parameterisations in the
	119	! universal model. Originally adapted from the Mars LMDZ model.
	120	!
	121	! The model can be run without or with tracer transport
	122	! depending on the value of "tracer" in file "callphys.def".
	123	!
	124	!
	125	! It includes:
	126	!
[1477]	127	! I. Initialization :
	128	! I.1 Firstcall initializations.
	129	! I.2 Initialization for every call to physiq.
[253]	130	!
[1477]	131	! II. Compute radiative transfer tendencies (longwave and shortwave) :
	132	! II.a Option 1 : Call correlated-k radiative transfer scheme.
	133	! II.b Option 2 : Call Newtonian cooling scheme.
	134	! II.c Option 3 : Atmosphere has no radiative effect.
	135	!
	136	! III. Vertical diffusion (turbulent mixing) :
	137	!
[2060]	138	! IV. Convection :
	139	! IV.a Thermal plume model
	140	! IV.b Dry convective adjusment
[1477]	141	!
	142	! V. Condensation and sublimation of gases (currently just CO2).
	143	!
	144	! VI. Tracers
	145	! VI.1. Water and water ice.
[1801]	146	! VI.2 Photochemistry
	147	! VI.3. Aerosols and particles.
	148	! VI.4. Updates (pressure variations, surface budget).
	149	! VI.5. Slab Ocean.
	150	! VI.6. Surface Tracer Update.
[1477]	151	!
	152	! VII. Surface and sub-surface soil temperature.
	153	!
	154	! VIII. Perform diagnostics and write output files.
	155	!
	156	!
[253]	157	! arguments
	158	! ---------
	159	!
[1477]	160	! INPUT
[253]	161	! -----
[1477]	162	!
[253]	163	! ngrid Size of the horizontal grid.
	164	! nlayer Number of vertical layers.
[1477]	165	! nq Number of advected fields.
	166	!
	167	! firstcall True at the first call.
	168	! lastcall True at the last call.
	169	!
	170	! pday Number of days counted from the North. Spring equinoxe.
	171	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
	172	! ptimestep timestep (s).
	173	!
	174	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa).
	175	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
	176	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2.s-2).
	177	!
	178	! pu(ngrid,nlayer) u, zonal component of the wind (ms-1).
	179	! pv(ngrid,nlayer) v, meridional component of the wind (ms-1).
	180	!
	181	! pt(ngrid,nlayer) Temperature (K).
	182	!
	183	! pq(ngrid,nlayer,nq) Advected fields.
	184	!
[1216]	185	! pudyn(ngrid,nlayer) \
[253]	186	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
[1477]	187	! ptdyn(ngrid,nlayer) / corresponding variables.
[253]	188	! pqdyn(ngrid,nlayer,nq) /
[1312]	189	! flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
[253]	190	!
[1477]	191	! OUTPUT
[253]	192	! ------
	193	!
[1308]	194	! pdu(ngrid,nlayer) \
	195	! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding
	196	! pdt(ngrid,nlayer) / variables due to physical processes.
	197	! pdq(ngrid,nlayer) /
[253]	198	! pdpsrf(ngrid) /
	199	!
	200	!
	201	! Authors
	202	! -------
[1524]	203	! Frederic Hourdin 15/10/93
	204	! Francois Forget 1994
[3232]	205	! Christophe Hourdin 02/1997
[253]	206	! Subroutine completely rewritten by F. Forget (01/2000)
	207	! Water ice clouds: Franck Montmessin (update 06/2003)
	208	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
	209	! New correlated-k radiative scheme: R. Wordsworth (2009)
[3232]	210	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
[253]	211	! Improved water cycle: R. Wordsworth / B. Charnay (2010)
	212	! To F90: R. Wordsworth (2010)
[594]	213	! New turbulent diffusion scheme: J. Leconte (2012)
[716]	214	! Loops converted to F90 matrix format: J. Leconte (2012)
[787]	215	! No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
[1477]	216	! Purge of the code : M. Turbet (2015)
[1801]	217	! Photochemical core developped by F. Lefevre: B. Charnay (2017)
[253]	218	!==================================================================
	219
	220
	221	! 0. Declarations :
	222	! ------------------
	223
[1669]	224	include "netcdf.inc"
[253]	225
	226	! Arguments :
	227	! -----------
	228
[1477]	229	! INPUTS:
[253]	230	! -------
	231
[1477]	232	integer,intent(in) :: ngrid ! Number of atmospheric columns.
	233	integer,intent(in) :: nlayer ! Number of atmospheric layers.
	234	integer,intent(in) :: nq ! Number of tracers.
[3232]	235
[1477]	236	logical,intent(in) :: firstcall ! Signals first call to physics.
	237	logical,intent(in) :: lastcall ! Signals last call to physics.
[3232]	238
[1477]	239	real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0.
	240	real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
	241	real,intent(in) :: ptimestep ! Physics timestep (s).
	242	real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
	243	real,intent(in) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
	244	real,intent(in) :: pphi(ngrid,nlayer) ! Geopotential at mid-layer (m2s-2).
	245	real,intent(in) :: pu(ngrid,nlayer) ! Zonal wind component (m/s).
	246	real,intent(in) :: pv(ngrid,nlayer) ! Meridional wind component (m/s).
	247	real,intent(in) :: pt(ngrid,nlayer) ! Temperature (K).
	248	real,intent(in) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
	249	real,intent(in) :: flxw(ngrid,nlayer) ! Vertical mass flux (ks/s) at lower boundary of layer
[253]	250
[1477]	251	! OUTPUTS:
[253]	252	! --------
	253
[1477]	254	! Physical tendencies :
	255
	256	real,intent(out) :: pdu(ngrid,nlayer) ! Zonal wind tendencies (m/s/s).
	257	real,intent(out) :: pdv(ngrid,nlayer) ! Meridional wind tendencies (m/s/s).
	258	real,intent(out) :: pdt(ngrid,nlayer) ! Temperature tendencies (K/s).
	259	real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
	260	real,intent(out) :: pdpsrf(ngrid) ! Surface pressure tendency (Pa/s).
	261
[253]	262	! Local saved variables:
	263	! ----------------------
[1622]	264	integer,save :: day_ini ! Initial date of the run (sol since Ls=0).
	265	integer,save :: icount ! Counter of calls to physiq during the run.
	266	!$OMP THREADPRIVATE(day_ini,icount)
	267
[253]	268	! Local variables :
	269	! -----------------
	270
[3232]	271	! Aerosol (dust or ice) extinction optical depth at reference wavelength
[253]	272	! for the "naerkind" optically active aerosols:
	273
[2972]	274	real,save,allocatable :: aerosol(:,:,:) ! Aerosols
	275	!$OMP THREADPRIVATE(aerosol)
[1477]	276	real zh(ngrid,nlayer) ! Potential temperature (K).
	277	real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
[1877]	278	real omega(ngrid,nlayer) ! omega velocity (Pa/s, >0 when downward)
[1477]	279
[2721]	280	integer l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_vap, igcm_generic_ice
[2722]	281	logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
[3232]	282
[1477]	283	real zls ! Solar longitude (radians).
	284	real zlss ! Sub solar point longitude (radians).
	285	real zday ! Date (time since Ls=0, calculated in sols).
	286	real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers.
	287	real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries.
[253]	288
[2060]	289	! VARIABLES for the thermal plume model
[3232]	290
[2232]	291	real f(ngrid) ! Mass flux norm
	292	real fm(ngrid,nlayer+1) ! Mass flux
	293	real fm_bis(ngrid,nlayer) ! Recasted fm
	294	real entr(ngrid,nlayer) ! Entrainment
	295	real detr(ngrid,nlayer) ! Detrainment
[2105]	296	real dqevap(ngrid,nlayer,nq) ! water tracer mass mixing ratio variations due to evaporation
	297	real dtevap(ngrid,nlayer) ! temperature variation due to evaporation
	298	real zqtherm(ngrid,nlayer,nq) ! vapor mass mixing ratio after evaporation
	299	real zttherm(ngrid,nlayer) ! temperature after evaporation
[2232]	300	real fraca(ngrid,nlayer+1) ! Fraction of the surface that plumes occupies
	301	real zw2(ngrid,nlayer+1) ! Vertical speed
	302	real zw2_bis(ngrid,nlayer) ! Recasted zw2
[3232]	303
[1477]	304	! TENDENCIES due to various processes :
[253]	305
[1477]	306	! For Surface Temperature : (K/s)
	307	real zdtsurf(ngrid) ! Cumulated tendencies.
	308	real zdtsurfmr(ngrid) ! Mass_redistribution routine.
	309	real zdtsurfc(ngrid) ! Condense_co2 routine.
	310	real zdtsdif(ngrid) ! Turbdiff/vdifc routines.
	311	real zdtsurf_hyd(ngrid) ! Hydrol routine.
[3235]	312
	313	! For Thermosphere : (K/s)
	314	real zdtconduc(ngrid,nlayer)
[3232]	315
	316	! For Atmospheric Temperatures : (K/s)
[1477]	317	real dtlscale(ngrid,nlayer) ! Largescale routine.
[3232]	318	real dt_generic_condensation(ngrid,nlayer) ! condensation_generic routine.
[1477]	319	real zdtc(ngrid,nlayer) ! Condense_co2 routine.
[2060]	320	real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	321	real zdttherm(ngrid,nlayer) ! Calltherm routine.
[1477]	322	real zdtmr(ngrid,nlayer) ! Mass_redistribution routine.
	323	real zdtrain(ngrid,nlayer) ! Rain routine.
[2721]	324	real zdtrain_generic(ngrid,nlayer) ! Rain_generic routine.
[1477]	325	real dtmoist(ngrid,nlayer) ! Moistadj routine.
[3100]	326	real dt_ekman(ngrid,nslay), dt_hdiff(ngrid,nslay) ! Slab_ocean routine.
[1477]	327	real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
[2542]	328	real zdtchim(ngrid,nlayer) ! Calchim routine.
[3232]	329
[1477]	330	! For Surface Tracers : (kg/m2/s)
	331	real dqsurf(ngrid,nq) ! Cumulated tendencies.
[1484]	332	real zdqsurfc(ngrid) ! Condense_co2 routine.
[1477]	333	real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines.
	334	real zdqssed(ngrid,nq) ! Callsedim routine.
	335	real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine.
	336	real zdqsrain(ngrid), zdqssnow(ngrid) ! Rain routine.
[2721]	337	real zdqsrain_generic(ngrid,nq), zdqssnow_generic(ngrid,nq) ! Rain_generic routine.
[1477]	338	real dqs_hyd(ngrid,nq) ! Hydrol routine.
[1859]	339	real reevap_precip(ngrid) ! re-evaporation flux of precipitation (integrated over the atmospheric column)
[2721]	340	real reevap_precip_generic(ngrid,nq)
[3232]	341
[1477]	342	! For Tracers : (kg/kg_of_air/s)
	343	real zdqc(ngrid,nlayer,nq) ! Condense_co2 routine.
	344	real zdqadj(ngrid,nlayer,nq) ! Convadj routine.
	345	real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines.
	346	real zdqevap(ngrid,nlayer) ! Turbdiff routine.
[2129]	347	real zdqtherm(ngrid,nlayer,nq) ! Calltherm routine.
[1477]	348	real zdqsed(ngrid,nlayer,nq) ! Callsedim routine.
	349	real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine.
	350	real zdqrain(ngrid,nlayer,nq) ! Rain routine.
	351	real dqmoist(ngrid,nlayer,nq) ! Moistadj routine.
	352	real dqvaplscale(ngrid,nlayer) ! Largescale routine.
	353	real dqcldlscale(ngrid,nlayer) ! Largescale routine.
[2721]	354	real dqvaplscale_generic(ngrid,nlayer,nq) ! condensation_generic routine.
	355	real dqcldlscale_generic(ngrid,nlayer,nq) ! condensation_generic routine.
	356	real dq_rain_generic_vap(ngrid,nlayer,nq) ! rain_generic routine
	357	real dq_rain_generic_cld(ngrid,nlayer,nq) ! rain_generic routine
[2542]	358	REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine
	359	REAL,allocatable,save :: zdqschim(:,:) ! Calchim_asis routine
[2972]	360	!$OMP THREADPRIVATE(zdqchim,zdqschim)
[3299]	361	real zdqvolc(ngrid,nlayer,nq) ! injection by volcanoes (kg/kg_of_air/s)
[1801]	362
	363	REAL array_zero1(ngrid)
	364	REAL array_zero2(ngrid,nlayer)
[3232]	365
[1477]	366	! For Winds : (m/s/s)
[2060]	367	real zdvadj(ngrid,nlayer), zduadj(ngrid,nlayer) ! Convadj routine.
	368	real zdutherm(ngrid,nlayer), zdvtherm(ngrid,nlayer) ! Calltherm routine.
	369	real zdumr(ngrid,nlayer), zdvmr(ngrid,nlayer) ! Mass_redistribution routine.
	370	real zdvdif(ngrid,nlayer), zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	371	real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
	372	real zdhadj(ngrid,nlayer) ! Convadj routine.
[3232]	373
[1477]	374	! For Pressure and Mass :
	375	real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
	376	real zdmassmr_col(ngrid) ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
	377	real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s).
[253]	378
[3232]	379
	380
[1477]	381	! Local variables for LOCAL CALCULATIONS:
	382	! ---------------------------------------
[787]	383	real zflubid(ngrid)
[1308]	384	real zplanck(ngrid),zpopsk(ngrid,nlayer)
[253]	385	real ztim1,ztim2,ztim3, z1,z2
	386	real ztime_fin
[1308]	387	real zdh(ngrid,nlayer)
[1194]	388	real gmplanet
[1297]	389	real taux(ngrid),tauy(ngrid)
[1194]	390
[253]	391
[1477]	392	! local variables for DIAGNOSTICS : (diagfi & stat)
	393	! -------------------------------------------------
	394	real ps(ngrid) ! Surface Pressure.
	395	real zt(ngrid,nlayer) ! Atmospheric Temperature.
	396	real zu(ngrid,nlayer),zv(ngrid,nlayer) ! Zonal and Meridional Winds.
	397	real zq(ngrid,nlayer,nq) ! Atmospheric Tracers.
	398	real zdtadj(ngrid,nlayer) ! Convadj Diagnostic.
	399	real zdtdyn(ngrid,nlayer) ! Dynamical Heating (K/s).
[1637]	400	real zdudyn(ngrid,nlayer) ! Dynamical Zonal Wind tendency (m.s-2).
[253]	401
[1477]	402	real reff(ngrid,nlayer) ! Effective dust radius (used if doubleq=T).
	403	real vmr(ngrid,nlayer) ! volume mixing ratio
[253]	404	real time_phys
[3232]	405
[1477]	406	real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
[3232]	407
[1477]	408	real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
[253]	409
	410	! included by RW for H2O Manabe scheme
[1477]	411	real rneb_man(ngrid,nlayer) ! H2O cloud fraction (moistadj).
	412	real rneb_lsc(ngrid,nlayer) ! H2O cloud fraction (large scale).
[253]	413
	414
[594]	415	! to test energy conservation (RW+JL)
[1308]	416	real mass(ngrid,nlayer),massarea(ngrid,nlayer)
[651]	417	real dEtot, dEtots, AtmToSurf_TurbFlux
[959]	418	real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
[1315]	419	!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
[3232]	420
[594]	421	!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
[1477]	422
[3232]	423	real dtmoist_max,dtmoist_min
[1295]	424	real dItot, dItot_tmp, dVtot, dVtot_tmp
[253]	425
	426
[1477]	427	real h2otot ! Total Amount of water. For diagnostic.
	428	real icesrf,liqsrf,icecol,vapcol ! Total Amounts of water (ice,liq,vap). For diagnostic.
[1295]	429	real dWtot, dWtot_tmp, dWtots, dWtots_tmp
[1477]	430	logical,save :: watertest
[1315]	431	!$OMP THREADPRIVATE(watertest)
[253]	432
[1477]	433	real qsat(ngrid,nlayer) ! Water Vapor Volume Mixing Ratio at saturation (kg/kg_of_air).
	434	real RH(ngrid,nlayer) ! Relative humidity.
	435	real psat_tmp
[2724]	436
	437	real qsat_generic(ngrid,nlayer,nq) ! generic condensable tracers (GCS) specific concentration at saturation (kg/kg_of_air).
	438	real RH_generic(ngrid,nlayer,nq) ! generic condensable tracers (GCS) Relative humidity.
	439	real rneb_generic(ngrid,nlayer,nq) ! GCS cloud fraction (generic condensation).
	440	real psat_tmp_generic
[2726]	441	real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
	442	!$OMP THREADPRIVATE(metallicity)
	443
[2890]	444	real reffrad_generic_zeros_for_wrf(ngrid,nlayer) ! !!! this is temporary, it is only a list of zeros, it will be replaced when a generic aerosol will be implemented
	445
[1477]	446	logical clearsky ! For double radiative transfer call. By BC
[3232]	447
[1482]	448	! For Clear Sky Case.
[2485]	449	real fluxsurfabs_sw1(ngrid) ! For SW/LW flux.
[1482]	450	real tau_col1(ngrid) ! For aerosol optical depth diagnostic.
[2537]	451	real OLR_nu1(ngrid,L_NSPECTI) ! Clear sky TOA LW radiation in each IR band
	452	real OSR_nu1(ngrid,L_NSPECTV) ! Clear sky TOA SW radiation in each VI band
	453	real GSR_nu1(ngrid,L_NSPECTV) ! Clear sky Surface SW radiation in each VI band
[2133]	454	real int_dtaui1(ngrid,nlayer,L_NSPECTI),int_dtauv1(ngrid,nlayer,L_NSPECTV) ! For optical thickness diagnostics.
[253]	455	real tf, ntf
[3269]	456
	457	real net_fluxsurf_lw(ngrid) ! net longwave flux at the surface (for diagnostics only)
[253]	458
[1477]	459	real nconsMAX, vdifcncons(ngrid), cadjncons(ngrid) ! Vdfic water conservation test. By RW
[253]	460
[1477]	461	real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW
[3233]	462
	463	! For fixed variable molar mass
	464	real, dimension(:),allocatable,save :: p_var
	465	real, dimension(:),allocatable,save :: mu_var
	466	real, dimension(:,:),allocatable,save :: frac_var
	467	!$OMP THREADPRIVATE(p_var,mu_var,frac_var)
	468	real :: mu_vari(nlayer)
	469	real :: muvari(ngrid,nlayer)
	470	integer ios,k
	471	character*100 dt_file
[253]	472
[2972]	473	real,save,allocatable :: reffcol(:,:)
	474	!$OMP THREADPRIVATE(reffcol)
[253]	475
[2299]	476	! Non-oro GW tendencies
	477	REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
	478	REAL d_t_hin(ngrid,nlayer)
	479	! Diagnostics 2D of gw_nonoro
	480	REAL zustrhi(ngrid), zvstrhi(ngrid)
[996]	481
[1297]	482
	483	real :: tsurf2(ngrid)
[3100]	484	!! real :: flux_o(ngrid),flux_g(ngrid)
	485	real :: flux_g(ngrid)
[1297]	486	real :: flux_sens_lat(ngrid)
	487	real :: qsurfint(ngrid,nq)
[3232]	488	#ifdef MESOSCALE
[2019]	489	REAL :: lsf_dt(nlayer)
	490	REAL :: lsf_dq(nlayer)
	491	#endif
[1297]	492
[2663]	493	! flags to trigger extra sanity checks
	494	logical, save :: check_physics_inputs=.false.
	495	logical, save :: check_physics_outputs=.false.
[3232]	496	!$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)
[2663]	497
[2133]	498	! Misc
	499	character*2 :: str2
[3232]	500	character(len=10) :: tmp1
[2736]	501	character(len=10) :: tmp2
[1477]	502	!==================================================================================================
[253]	503
	504	! -----------------
[1477]	505	! I. INITIALISATION
	506	! -----------------
[253]	507
[1477]	508	! --------------------------------
	509	! I.1 First Call Initialisation.
	510	! --------------------------------
[253]	511	if (firstcall) then
[2867]	512	call getin_p("check_physics_inputs", check_physics_inputs)
	513	call getin_p("check_physics_outputs", check_physics_outputs)
[2663]	514
[2867]	515	! Allocate saved arrays (except for 1D model, where this has already
	516	! been done)
[2010]	517	#ifndef MESOSCALE
[2867]	518	if (ngrid>1) call phys_state_var_init(nq)
[2010]	519	#endif
[858]	520
[1477]	521	! Variables set to 0
[253]	522	! ~~~~~~~~~~~~~~~~~~
	523	dtrad(:,:) = 0.0
	524	fluxrad(:) = 0.0
	525	tau_col(:) = 0.0
	526	zdtsw(:,:) = 0.0
	527	zdtlw(:,:) = 0.0
[3293]	528	zdtconduc(:,:) = 0.0
[3233]	529
	530	! Initialize fixed variable mu
	531	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[726]	532
[3233]	533	if(varspec) then
	534	IF (.NOT.ALLOCATED(p_var)) ALLOCATE(p_var(nvarlayer))
	535	IF (.NOT.ALLOCATED(mu_var)) ALLOCATE(mu_var(nvarlayer))
	536	IF (.NOT.ALLOCATED(frac_var)) ALLOCATE(frac_var(nvarlayer,ngasmx))
	537	p_var = 0.0
	538	mu_var = 0.0
	539	frac_var = 0.0
	540	muvari = 0.0
	541	dt_file=TRIM(varspec_data)
	542	open(33,file=dt_file,form='formatted',status='old',iostat=ios)
	543	if (ios.ne.0) then ! file not found
	544	write(,) 'Error from varspec'
	545	write(,) 'Data file ',trim(varspec_data),' not found.'
	546	write(,) 'Check that the data is in your run repository.'
	547	call abort_physic !a verifier
	548	else
	549	do k=1,nvarlayer
	550	read(33,*) p_var(k), mu_var(k),frac_var(k,1:ngasmx)
	551	!The order of columns in frac_var must correspond to order of molecules gases.def
	552	!The format of your file must be:
	553	! pressure(k) molar_mass(k), molar_fraction_of_gas_1(k), molar_fraction_of_gas_2(k),..., molar_fraction_of_gas_ngasmx(k)
	554	enddo
	555	endif
	556	close(33)
[3236]	557	else
	558	IF (.NOT.ALLOCATED(p_var)) ALLOCATE(p_var(nlayer))
	559	IF (.NOT.ALLOCATED(mu_var)) ALLOCATE(mu_var(nlayer))
	560	IF (.NOT.ALLOCATED(frac_var)) ALLOCATE(frac_var(nlayer,ngasmx))
[3233]	561	endif
	562
[1477]	563	! Initialize tracer names, indexes and properties.
	564	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	565	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
[253]	566	if (tracer) then
[2785]	567	call initracer(ngrid,nq)
[1801]	568	if(photochem) then
	569	call ini_conc_mod(ngrid,nlayer)
[2542]	570	IF (.NOT.ALLOCATED(zdqchim)) ALLOCATE(zdqchim(ngrid,nlayer,nesp))
	571	IF (.NOT.ALLOCATED(zdqschim)) ALLOCATE(zdqschim(ngrid,nesp))
[1801]	572	endif
[1477]	573	endif
[2803]	574	! Initialize aerosol indexes.
	575	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2898]	576	call iniaerosol
[2972]	577	! allocate related local arrays
	578	! (need be allocated instead of automatic because of "naerkind")
	579	allocate(aerosol(ngrid,nlayer,naerkind))
	580	allocate(reffcol(ngrid,naerkind))
[253]	581
[1682]	582	#ifdef CPP_XIOS
[2867]	583	! Initialize XIOS context
	584	write(,) "physiq: call wxios_context_init"
	585	CALL wxios_context_init
[1682]	586	#endif
[726]	587
[1477]	588	! Read 'startfi.nc' file.
[253]	589	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3232]	590	#ifndef MESOSCALE
[1669]	591	call phyetat0(startphy_file, &
	592	ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
[1477]	593	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, &
	594	cloudfrac,totcloudfrac,hice, &
	595	rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
[3100]	596
[3232]	597	#else
[3100]	598
[2867]	599	day_ini = pday
[1836]	600	#endif
	601
	602	#ifndef MESOSCALE
[1669]	603	if (.not.startphy_file) then
	604	! additionnal "academic" initialization of physics
	605	if (is_master) write(,) "Physiq: initializing tsurf(:) to pt(:,1) !!"
	606	tsurf(:)=pt(:,1)
	607	if (is_master) write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
	608	do isoil=1,nsoilmx
	609	tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
	610	enddo
[3232]	611	if (is_master) write(,) "Physiq: initializing day_ini to pday !"
[1669]	612	day_ini=pday
	613	endif
[1836]	614	#endif
[253]	615	if (pday.ne.day_ini) then
[2299]	616	write(,) "ERROR in physiq.F90:"
	617	write(,) "bad synchronization between physics and dynamics"
	618	write(,) "dynamics day: ",pday
	619	write(,) "physics day: ",day_ini
	620	stop
[253]	621	endif
	622
	623	write (,) 'In physiq day_ini =', day_ini
	624
[3232]	625	! Initialize albedo calculation.
[1482]	626	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	627	albedo(:,:)=0.0
[2867]	628	albedo_bareground(:)=0.0
	629	albedo_snow_SPECTV(:)=0.0
	630	albedo_co2_ice_SPECTV(:)=0.0
[3232]	631	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
[3100]	632
[3232]	633	! Initialize orbital calculation.
[1482]	634	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[253]	635	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
	636
	637
	638	if(tlocked)then
	639	print*,'Planet is tidally locked at resonance n=',nres
	640	print*,'Make sure you have the right rotation rate!!!'
	641	endif
	642
[1477]	643	! Initialize oceanic variables.
	644	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1297]	645
	646	if (ok_slab_ocean)then
	647
[3100]	648	!! call ocean_slab_init(ngrid,ptimestep, tslab, &
	649	!! sea_ice, pctsrf_sic)
[1297]	650
[3100]	651	!! call ini_surf_heat_transp_mod()
[3232]	652
[1477]	653	knindex(:) = 0
[3100]	654	knon = 0
[1297]	655
[1477]	656	do ig=1,ngrid
[1297]	657	zmasq(ig)=1
[3100]	658	! knindex(ig) = ig
	659	if (nint(rnat(ig)).eq.0) then ! rnat = 0 (ocean), 1 (continent)
[1297]	660	zmasq(ig)=0
[3100]	661	knon=knon+1
	662	knindex(knon)=ig
[1297]	663	endif
[1477]	664	enddo
[1297]	665
[3100]	666	call ocean_slab_init(ptimestep, pctsrf_sic, tslab, tsea_ice, sea_ice, zmasq)
[1297]	667
[3100]	668	!! CALL init_masquv(ngrid,zmasq)
	669
[1477]	670	endif ! end of 'ok_slab_ocean'.
[1297]	671
	672
[1477]	673	! Initialize soil.
	674	! ~~~~~~~~~~~~~~~~
[253]	675	if (callsoil) then
[3232]	676
[787]	677	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
[1477]	678	ptimestep,tsurf,tsoil,capcal,fluxgrd)
[1297]	679
	680	if (ok_slab_ocean) then
	681	do ig=1,ngrid
	682	if (nint(rnat(ig)).eq.2) then
[1477]	683	capcal(ig)=capcalocean
	684	if (pctsrf_sic(ig).gt.0.5) then
	685	capcal(ig)=capcalseaice
	686	if (qsurf(ig,igcm_h2o_ice).gt.0.) then
	687	capcal(ig)=capcalsno
	688	endif
	689	endif
[1297]	690	endif
	691	enddo
[1477]	692	endif ! end of 'ok_slab_ocean'.
[1297]	693
[1477]	694	else ! else of 'callsoil'.
[3232]	695
[253]	696	print*,'WARNING! Thermal conduction in the soil turned off'
[918]	697	capcal(:)=1.e6
[952]	698	fluxgrd(:)=intheat
	699	print*,'Flux from ground = ',intheat,' W m^-2'
[3232]	700
[1477]	701	endif ! end of 'callsoil'.
[3232]	702
[253]	703	icount=1
	704
[1477]	705	! Here is defined the type of the surface : Continent or Ocean.
	706	! BC2014 : This is now already done in newstart.
	707	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1297]	708	if (.not.ok_slab_ocean) then
[3232]	709
[1297]	710	rnat(:)=1.
[2603]	711	if (.not.nosurf) then ! inertiedat only defined if there is a surface
	712	do ig=1,ngrid
[1477]	713	if(inertiedat(ig,1).gt.1.E4)then
	714	rnat(ig)=0
	715	endif
[2603]	716	enddo
[253]	717
[2603]	718	print*,'WARNING! Surface type currently decided by surface inertia'
	719	print*,'This should be improved e.g. in newstart.F'
	720	endif
[1477]	721	endif ! end of 'ok_slab_ocean'.
[253]	722
[1477]	723
	724	! Initialize surface history variable.
[253]	725	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[728]	726	qsurf_hist(:,:)=qsurf(:,:)
[253]	727
[1637]	728	! Initialize variable for dynamical heating and zonal wind tendency diagnostic
	729	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[253]	730	ztprevious(:,:)=pt(:,:)
[1637]	731	zuprevious(:,:)=pu(:,:)
[253]	732
	733	! Set temperature just above condensation temperature (for Early Mars)
	734	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	735	if(nearco2cond) then
	736	write(,)' WARNING! Starting at Tcond+1K'
	737	do l=1, nlayer
	738	do ig=1,ngrid
	739	pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4 &
	740	-pt(ig,l)) / ptimestep
	741	enddo
	742	enddo
	743	endif
	744
[3232]	745	if(meanOLR)then
	746	call system('rm -f rad_bal.out') ! to record global radiative balance.
	747	call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
[1477]	748	call system('rm -f h2o_bal.out') ! to record global hydrological balance.
[253]	749	endif
	750
	751
[3232]	752	watertest=.false.
[1477]	753	if(water)then ! Initialize variables for water cycle
[3232]	754
[365]	755	if(enertest)then
	756	watertest = .true.
	757	endif
	758
[253]	759	endif
[2726]	760
	761	! Set metallicity for GCS
	762	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	763	metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
	764	call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic
[3232]	765
[2060]	766	! Set some parameters for the thermal plume model
	767	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	768	if (calltherm) then
[2176]	769	CALL init_thermcell_mod(g, rcp, r, pi, T_h2o_ice_liq, RV)
[2060]	770	endif
[3100]	771
[253]	772	call su_watercycle ! even if we don't have a water cycle, we might
	773	! need epsi for the wvp definitions in callcorrk.F
[2060]	774	! or RETV, RLvCp for the thermal plume model
[3100]	775
[1836]	776	#ifndef MESOSCALE
[1477]	777	if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
[1542]	778	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
[2299]	779	ptimestep,pday+nday,time_phys,cell_area, &
	780	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
[1216]	781	endif
[3100]	782
[3232]	783	#endif
	784	if (corrk) then
	785	! We initialise the spectral grid here instead of
	786	! at firstcall of callcorrk so we can output XspecIR, XspecVI
[2736]	787	! when using Dynamico
	788	print*, "physiq_mod: Correlated-k data base folder:",trim(datadir)
	789	call getin_p("corrkdir",corrkdir)
	790	print*,"corrkdir = ", corrkdir
[2802]	791	write (tmp1, '(i4)') L_NSPECTI
	792	write (tmp2, '(i4)') L_NSPECTV
[2736]	793	banddir=trim(trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2)))
	794	banddir=trim(trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir)))
	795	call setspi !Basic infrared properties.
	796	call setspv ! Basic visible properties.
	797	call sugas_corrk ! Set up gaseous absorption properties.
	798	call suaer_corrk ! Set up aerosol optical properties.
	799	endif
[3100]	800
[1622]	801	! XIOS outputs
	802	#ifdef CPP_XIOS
	803
	804	write(,) "physiq: call initialize_xios_output"
	805	call initialize_xios_output(pday,ptime,ptimestep,daysec, &
[2736]	806	year_day,presnivs,pseudoalt,WNOI,WNOV)
[1622]	807	#endif
[3232]	808
[1682]	809	write(,) "physiq: end of firstcall"
[1477]	810	endif ! end of 'firstcall'
[253]	811
[3100]	812
[2663]	813	if (check_physics_inputs) then
	814	!check the validity of input fields coming from the dynamics
	815	call check_physics_fields("begin physiq:", pt, pu, pv, pplev, pq)
	816	endif
	817
[3100]	818
[1477]	819	! ------------------------------------------------------
	820	! I.2 Initializations done at every physical timestep:
	821	! ------------------------------------------------------
	822
[3232]	823	#ifdef CPP_XIOS
[1622]	824	! update XIOS time/calendar
	825	call update_xios_timestep
[3232]	826	#endif
[1622]	827
[1477]	828	! Initialize various variables
	829	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3232]	830
[253]	831	if ( .not.nearco2cond ) then
[1308]	832	pdt(1:ngrid,1:nlayer) = 0.0
[3232]	833	endif
[1477]	834	zdtsurf(1:ngrid) = 0.0
[1308]	835	pdq(1:ngrid,1:nlayer,1:nq) = 0.0
[1477]	836	dqsurf(1:ngrid,1:nq)= 0.0
	837	pdu(1:ngrid,1:nlayer) = 0.0
	838	pdv(1:ngrid,1:nlayer) = 0.0
[787]	839	pdpsrf(1:ngrid) = 0.0
[3232]	840	zflubid(1:ngrid) = 0.0
[1297]	841	flux_sens_lat(1:ngrid) = 0.0
	842	taux(1:ngrid) = 0.0
	843	tauy(1:ngrid) = 0.0
[3294]	844	net_fluxsurf_lw(1:ngrid) = 0.0
[253]	845
[1477]	846	zday=pday+ptime ! Compute time, in sols (and fraction thereof).
[1297]	847
[1477]	848	! Compute Stellar Longitude (Ls), and orbital parameters.
	849	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[253]	850	if (season) then
	851	call stellarlong(zday,zls)
	852	else
[2300]	853	call stellarlong(noseason_day,zls)
[253]	854	end if
	855
[1329]	856	call orbite(zls,dist_star,declin,right_ascen)
[3232]	857
[1329]	858	if (tlocked) then
	859	zlss=Mod(-(2.pi(zday/year_day)nres - right_ascen),2.pi)
	860	elseif (diurnal) then
[1524]	861	zlss=-2.pi(zday-.5)
[1329]	862	else if(diurnal .eqv. .false.) then
	863	zlss=9999.
[3232]	864	endif
[1194]	865
	866
[1477]	867	! Compute variations of g with latitude (oblate case).
	868	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3232]	869	if (oblate .eqv. .false.) then
	870	glat(:) = g
	871	else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
[1477]	872	print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
	873	call abort
	874	else
	875	gmplanet = MassPlanetgrav1e24
	876	do ig=1,ngrid
[3232]	877	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]	878	end do
	879	endif
[1297]	880
[1477]	881	! Compute geopotential between layers.
	882	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1308]	883	zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
[3232]	884	do l=1,nlayer
[1477]	885	zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
[1194]	886	enddo
[728]	887
[787]	888	zzlev(1:ngrid,1)=0.
[728]	889
[253]	890	do l=2,nlayer
	891	do ig=1,ngrid
	892	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
	893	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
	894	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	895	enddo
[3232]	896	enddo
[253]	897
[2987]	898	!Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
[3232]	899
[2987]	900	zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
	901
[1477]	902	! Compute potential temperature
	903	! Note : Potential temperature calculation may not be the same in physiq and dynamic...
	904	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3232]	905	do l=1,nlayer
[787]	906	do ig=1,ngrid
[253]	907	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
[597]	908	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
[1194]	909	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
[1542]	910	massarea(ig,l)=mass(ig,l)*cell_area(ig)
[253]	911	enddo
	912	enddo
	913
[1312]	914	! Compute vertical velocity (m/s) from vertical mass flux
[1346]	915	! w = F / (rhoarea) and rho = P/(rT)
[1477]	916	! But first linearly interpolate mass flux to mid-layers
	917	do l=1,nlayer-1
	918	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
	919	enddo
	920	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
	921	do l=1,nlayer
	922	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) / &
[1542]	923	(pplay(1:ngrid,l)*cell_area(1:ngrid))
[1477]	924	enddo
[1877]	925	! omega in Pa/s
	926	do l=1,nlayer-1
	927	omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
	928	enddo
	929	omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
	930	do l=1,nlayer
	931	omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid)
	932	enddo
[1194]	933
[1801]	934	! ----------------------------------------------------------------
	935	! Compute mean mass, cp, and R
	936	! --------------------------------
[2058]	937	#ifndef MESOSCALE
[1801]	938	if(photochem) then
	939	call concentrations(ngrid,nlayer,nq,pplay,pt,pdt,pq,pdq,ptimestep)
[3232]	940	endif
[2058]	941	#endif
[1801]	942
[1477]	943	!---------------------------------
	944	! II. Compute radiative tendencies
	945	!---------------------------------
[253]	946
[3081]	947	! Compute local stellar zenith angles
	948	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	949	if (tlocked) then
	950	! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb
	951	ztim1=SIN(declin)
	952	ztim2=COS(declin)*COS(zlss)
	953	ztim3=COS(declin)*SIN(zlss)
[253]	954
[3081]	955	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
	956	ztim1,ztim2,ztim3,mu0,fract, flatten)
[253]	957
[3081]	958	elseif (diurnal) then
	959	ztim1=SIN(declin)
	960	ztim2=COS(declin)COS(2.pi*(zday-.5))
	961	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
[253]	962
[3081]	963	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
	964	ztim1,ztim2,ztim3,mu0,fract, flatten)
	965	else if(diurnal .eqv. .false.) then
[253]	966
[3081]	967	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad,flatten)
	968	! WARNING: this function appears not to work in 1D
[253]	969
[3081]	970	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
	971	mu0 = cos(pi*szangle/180.0)
	972	endif
[253]	973
[3232]	974	endif
[2470]	975
[3081]	976	if (callrad) then
	977	if( mod(icount-1,iradia).eq.0.or.lastcall) then
[3232]	978
	979	! Eclipse incoming sunlight (e.g. Saturn ring shadowing).
[1429]	980	if(rings_shadow) then
	981	call call_rings(ngrid, ptime, pday, diurnal)
[3232]	982	endif
[1133]	983
[1329]	984
[1477]	985	if (corrk) then
[3232]	986
[1477]	987	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	988	! II.a Call correlated-k radiative transfer scheme
	989	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	990	if(kastprof)then
	991	print*,'kastprof should not = true here'
	992	call abort
	993	endif
[1524]	994	if(water) then
[3233]	995	! take into account water effect on mean molecular weight
[3232]	996	muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap))
	997	muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap))
[2728]	998	elseif(generic_condensation) then
[3233]	999	! take into account generic condensable specie (GCS) effect on mean molecular weight
[2728]	1000	do iq=1,nq
	1001
	1002	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
[3232]	1003
[2728]	1004	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
	1005
	1006	epsi_generic=constants_epsi_generic(iq)
	1007
[3232]	1008	muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap))
	1009	muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap))
[2728]	1010
	1011	endif
[3232]	1012	end do ! do iq=1,nq loop on tracers
[3233]	1013	elseif(varspec) then
	1014	!take into account fixed variable mean molecular weight
	1015	do ig=1,ngrid
	1016	call pindex(p_var,mu_var,pplay(ig,:),nvarlayer,nlayer,mu_vari)
	1017	muvari(ig,:) = mu_vari
	1018	enddo
[3232]	1019
[3233]	1020	muvar(1:ngrid,1:nlayer) = muvari(1:ngrid,1:nlayer)
	1021	muvar(1:ngrid,nlayer+1) = muvari(1:ngrid,nlayer)
	1022
[1524]	1023	else
[1308]	1024	muvar(1:ngrid,1:nlayer+1)=mugaz
[3232]	1025	endif
[1297]	1026
[3232]	1027
[1477]	1028	if(ok_slab_ocean) then
[3232]	1029	tsurf2(:)=tsurf(:)
[1477]	1030	do ig=1,ngrid
	1031	if (nint(rnat(ig))==0) then
	1032	tsurf(ig)=((1.-pctsrf_sic(ig))tslab(ig,1)4+pctsrf_sic(ig)tsea_ice(ig)4)0.25
	1033	endif
	1034	enddo
	1035	endif !(ok_slab_ocean)
[3232]	1036
[1477]	1037	! standard callcorrk
	1038	clearsky=.false.
[1482]	1039	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
	1040	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
	1041	tsurf,fract,dist_star,aerosol,muvar, &
	1042	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, &
	1043	fluxsurfabs_sw,fluxtop_lw, &
[2537]	1044	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,GSR_nu, &
[2133]	1045	int_dtaui,int_dtauv, &
[1482]	1046	tau_col,cloudfrac,totcloudfrac, &
[3233]	1047	clearsky,p_var,frac_var,firstcall,lastcall)
[1297]	1048
[2831]	1049	!GG (feb2021): Option to "artificially" decrease the raditive time scale in
	1050	!the deep atmosphere press > 0.1 bar. Suggested by MT.
[3232]	1051	!! COEFF_RAD to be "tuned" to facilitate convergence of tendency
	1052
[2831]	1053	!coeff_rad=0. ! 0 values, it doesn't accelerate the convergence
	1054	!coeff_rad=0.5
[3232]	1055	!coeff_rad=1.
[2831]	1056	!do l=1, nlayer
	1057	! do ig=1,ngrid
	1058	! if(pplay(ig,l).ge.1.d4) then
	1059	! zdtsw(ig,l)=zdtsw(ig,l)(pplay(ig,l)/1.d4)*coeff_rad
	1060	! zdtlw(ig,l)=zdtlw(ig,l)(pplay(ig,l)/1.d4)*coeff_rad
	1061	! endif
	1062	! enddo
	1063	!enddo
	1064
[3232]	1065	! Option to call scheme once more for clear regions
[1477]	1066	if(CLFvarying)then
[253]	1067
[1477]	1068	! ---> PROBLEMS WITH ALLOCATED ARRAYS : temporary solution in callcorrk: do not deallocate if CLFvarying ...
	1069	clearsky=.true.
[1482]	1070	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
[2537]	1071	albedo,albedo_equivalent1,emis,mu0,pplev,pplay,pt, &
[1482]	1072	tsurf,fract,dist_star,aerosol,muvar, &
	1073	zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1, &
	1074	fluxsurfabs_sw1,fluxtop_lw1, &
[2537]	1075	fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,GSR_nu1, &
[2133]	1076	int_dtaui1,int_dtauv1, &
[1482]	1077	tau_col1,cloudfrac,totcloudfrac, &
[3233]	1078	clearsky,p_var,frac_var,firstcall,lastcall)
[3232]	1079	clearsky = .false. ! just in case.
[253]	1080
[1477]	1081	! Sum the fluxes and heating rates from cloudy/clear cases
	1082	do ig=1,ngrid
	1083	tf=totcloudfrac(ig)
[3232]	1084	ntf=1.-tf
[1482]	1085	fluxsurf_lw(ig) = ntffluxsurf_lw1(ig) + tffluxsurf_lw(ig)
	1086	fluxsurf_sw(ig) = ntffluxsurf_sw1(ig) + tffluxsurf_sw(ig)
	1087	albedo_equivalent(ig) = ntfalbedo_equivalent1(ig) + tfalbedo_equivalent(ig)
	1088	fluxsurfabs_sw(ig) = ntffluxsurfabs_sw1(ig) + tffluxsurfabs_sw(ig)
	1089	fluxtop_lw(ig) = ntffluxtop_lw1(ig) + tffluxtop_lw(ig)
	1090	fluxabs_sw(ig) = ntffluxabs_sw1(ig) + tffluxabs_sw(ig)
	1091	tau_col(ig) = ntftau_col1(ig) + tftau_col(ig)
[3232]	1092
[1477]	1093	zdtlw(ig,1:nlayer) = ntfzdtlw1(ig,1:nlayer) + tfzdtlw(ig,1:nlayer)
	1094	zdtsw(ig,1:nlayer) = ntfzdtsw1(ig,1:nlayer) + tfzdtsw(ig,1:nlayer)
[253]	1095
[3232]	1096	OSR_nu(ig,1:L_NSPECTV) = ntfOSR_nu1(ig,1:L_NSPECTV) + tfOSR_nu(ig,1:L_NSPECTV)
	1097	GSR_nu(ig,1:L_NSPECTV) = ntfGSR_nu1(ig,1:L_NSPECTV) + tfGSR_nu(ig,1:L_NSPECTV)
	1098	OLR_nu(ig,1:L_NSPECTI) = ntfOLR_nu1(ig,1:L_NSPECTI) + tfOLR_nu(ig,1:L_NSPECTI)
[2446]	1099	if (diagdtau) then
[3232]	1100	int_dtauv(ig,:,1:L_NSPECTV) = ntfint_dtauv1(ig,:,1:L_NSPECTV) + tfint_dtauv(ig,:,1:L_NSPECTV)
	1101	int_dtaui(ig,:,1:L_NSPECTI) = ntfint_dtaui1(ig,:,1:L_NSPECTI) + tfint_dtaui(ig,:,1:L_NSPECTI)
[2446]	1102	endif
[3232]	1103	enddo
[253]	1104
[1477]	1105	endif ! end of CLFvarying.
[253]	1106
[1477]	1107	if(ok_slab_ocean) then
	1108	tsurf(:)=tsurf2(:)
	1109	endif
[1297]	1110
	1111
[1482]	1112	! Radiative flux from the sky absorbed by the surface (W.m-2).
[1477]	1113	GSR=0.0
[1482]	1114	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurfabs_sw(1:ngrid)
[3269]	1115	net_fluxsurf_lw(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
[253]	1116
[1477]	1117	!if(noradsurf)then ! no lower surface; SW flux just disappears
[1542]	1118	! GSR = SUM(fluxsurf_sw(1:ngrid)*cell_area(1:ngrid))/totarea
[1477]	1119	! fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
	1120	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
	1121	!endif
[253]	1122
[1477]	1123	! Net atmospheric radiative heating rate (K.s-1)
	1124	dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
[3232]	1125
	1126	! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that !
[1498]	1127	if (firstcall .and. albedo_spectral_mode) then
	1128	call spectral_albedo_calc(albedo_snow_SPECTV)
	1129	endif
[253]	1130
[1477]	1131	elseif(newtonian)then
[3232]	1132
[1482]	1133	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1134	! II.b Call Newtonian cooling scheme
	1135	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1477]	1136	call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
[253]	1137
[1477]	1138	zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep
	1139	! e.g. surface becomes proxy for 1st atmospheric layer ?
[253]	1140
[1477]	1141	else
[3232]	1142
[1477]	1143	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3232]	1144	! II.c Atmosphere has no radiative effect
[1477]	1145	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1146	fluxtop_dn(1:ngrid) = fract(1:ngrid)mu0(1:ngrid)Fat1AU/dist_star**2
	1147	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
	1148	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
	1149	endif
	1150	fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
[1482]	1151	print*,'------------WARNING---WARNING------------' ! by MT2015.
	1152	print*,'You are in corrk=false mode, '
[3232]	1153	print*,'and the surface albedo is taken equal to the first visible spectral value'
	1154
[1498]	1155	fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
	1156	fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid)
[1477]	1157	fluxtop_lw(1:ngrid) = emis(1:ngrid)sigmatsurf(1:ngrid)**4
[253]	1158
[1477]	1159	dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating
[253]	1160
[1477]	1161	endif ! end of corrk
[253]	1162
[1477]	1163	endif ! of if(mod(icount-1,iradia).eq.0)
[253]	1164
[3232]	1165
[1477]	1166	! Transformation of the radiative tendencies
	1167	! ------------------------------------------
	1168	zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
	1169	zplanck(1:ngrid)=emis(1:ngrid)sigmazplanck(1:ngrid)*zplanck(1:ngrid)
[3269]	1170	net_fluxsurf_lw(1:ngrid)=net_fluxsurf_lw(1:ngrid)-zplanck(1:ngrid)
[1477]	1171	fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
	1172	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
[3232]	1173
[1477]	1174	! Test of energy conservation
	1175	!----------------------------
[253]	1176	if(enertest)then
[1524]	1177	call planetwide_sumval(cppmassarea(:,:)zdtsw(:,:)/totarea_planet,dEtotSW)
	1178	call planetwide_sumval(cppmassarea(:,:)zdtlw(:,:)/totarea_planet,dEtotLW)
[1542]	1179	!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
	1180	call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
	1181	call planetwide_sumval((fluxsurf_lw(:)emis(:)-zplanck(:))cell_area(:)/totarea_planet,dEtotsLW)
[1524]	1182	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
	1183	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
	1184	if (is_master) then
[1477]	1185	print*,'---------------------------------------------------------------'
	1186	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
	1187	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
	1188	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
	1189	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
	1190	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
	1191	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
[1524]	1192	endif
[1477]	1193	endif ! end of 'enertest'
[253]	1194
	1195	endif ! of if (callrad)
	1196
	1197
[1477]	1198	! --------------------------------------------
	1199	! III. Vertical diffusion (turbulent mixing) :
	1200	! --------------------------------------------
	1201
[253]	1202	if (calldifv) then
[3232]	1203
[787]	1204	zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
[253]	1205
[1477]	1206	! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
[1524]	1207	if (UseTurbDiff) then
[3232]	1208
[1477]	1209	call turbdiff(ngrid,nlayer,nq,rnat, &
[2427]	1210	ptimestep,capcal, &
[1477]	1211	pplay,pplev,zzlay,zzlev,z0, &
	1212	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
	1213	pdt,pdq,zflubid, &
	1214	zdudif,zdvdif,zdtdif,zdtsdif, &
[1524]	1215	sensibFlux,q2,zdqdif,zdqevap,zdqsdif, &
[2427]	1216	taux,tauy)
[594]	1217
[1524]	1218	else
[3232]	1219
[1477]	1220	zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
[3232]	1221
[1477]	1222	call vdifc(ngrid,nlayer,nq,rnat,zpopsk, &
	1223	ptimestep,capcal,lwrite, &
	1224	pplay,pplev,zzlay,zzlev,z0, &
	1225	pu,pv,zh,pq,tsurf,emis,qsurf, &
	1226	zdh,pdq,zflubid, &
	1227	zdudif,zdvdif,zdhdif,zdtsdif, &
[2427]	1228	sensibFlux,q2,zdqdif,zdqsdif)
[253]	1229
[1477]	1230	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
[1524]	1231	zdqevap(1:ngrid,1:nlayer)=0.
[594]	1232
[1477]	1233	end if !end of 'UseTurbDiff'
[594]	1234
[2121]	1235	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
	1236
[1836]	1237	!!! this is always done, except for turbulence-resolving simulations
	1238	if (.not. turb_resolved) then
	1239	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
	1240	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
	1241	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
	1242	endif
[1477]	1243
[1297]	1244	if(ok_slab_ocean)then
	1245	flux_sens_lat(1:ngrid)=(zdtsdif(1:ngrid)*capcal(1:ngrid)-fluxrad(1:ngrid))
	1246	endif
	1247
[3232]	1248	if (tracer) then
[1308]	1249	pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
[787]	1250	dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
[253]	1251	end if ! of if (tracer)
	1252
[1477]	1253
	1254	! test energy conservation
[253]	1255	!-------------------------
	1256	if(enertest)then
[3232]	1257
[1524]	1258	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
[253]	1259	do ig = 1, ngrid
[1524]	1260	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
	1261	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
[253]	1262	enddo
[3232]	1263
[1542]	1264	call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
[1524]	1265	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
[1542]	1266	call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
	1267	call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
[3232]	1268
[1524]	1269	if (is_master) then
[3232]	1270
[1477]	1271	if (UseTurbDiff) then
[1524]	1272	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
	1273	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
[1477]	1274	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
[1524]	1275	else
	1276	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
	1277	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
	1278	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
	1279	end if
	1280	endif ! end of 'is_master'
[3232]	1281
[1477]	1282	! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
	1283	endif ! end of 'enertest'
[253]	1284
[1477]	1285
	1286	! Test water conservation.
[253]	1287	if(watertest.and.water)then
[3232]	1288
[1524]	1289	call planetwide_sumval(massarea(:,:)zdqdif(:,:,igcm_h2o_vap)ptimestep/totarea_planet,dWtot_tmp)
[1542]	1290	call planetwide_sumval(zdqsdif(:,igcm_h2o_vap)cell_area(:)ptimestep/totarea_planet,dWtots_tmp)
[253]	1291	do ig = 1, ngrid
[1524]	1292	vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
	1293	enddo
	1294	call planetwide_sumval(massarea(:,:)zdqdif(:,:,igcm_h2o_ice)ptimestep/totarea_planet,dWtot)
[1542]	1295	call planetwide_sumval(zdqsdif(:,igcm_h2o_ice)cell_area(:)ptimestep/totarea_planet,dWtots)
[1524]	1296	dWtot = dWtot + dWtot_tmp
	1297	dWtots = dWtots + dWtots_tmp
[651]	1298	do ig = 1, ngrid
[1524]	1299	vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
[3232]	1300	enddo
[1524]	1301	call planetwide_maxval(vdifcncons(:),nconsMAX)
[253]	1302
[1524]	1303	if (is_master) then
[1477]	1304	print*,'---------------------------------------------------------------'
	1305	print*,'In difv atmospheric water change =',dWtot,' kg m-2'
	1306	print*,'In difv surface water change =',dWtots,' kg m-2'
	1307	print*,'In difv non-cons factor =',dWtot+dWtots,' kg m-2'
	1308	print*,'In difv MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
[1524]	1309	endif
[253]	1310
[1477]	1311	endif ! end of 'watertest'
[253]	1312	!-------------------------
	1313
[1477]	1314	else ! calldifv
[253]	1315
	1316	if(.not.newtonian)then
	1317
[787]	1318	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
[253]	1319
	1320	endif
	1321
[1477]	1322	endif ! end of 'calldifv'
[253]	1323
	1324
[2060]	1325	!-------------------
	1326	! IV. Convection :
	1327	!-------------------
[3232]	1328
[2060]	1329	! ~~~~~~~~~~~~~~~~~~~~~~~~~~
	1330	! IV.a Thermal plume model :
	1331	! ~~~~~~~~~~~~~~~~~~~~~~~~~~
[3232]	1332
[2060]	1333	IF (calltherm) THEN
[3232]	1334
[2127]	1335	! AB: We need to evaporate ice before calling thermcell_main.
[2105]	1336	IF (water) THEN
	1337	CALL evap(ngrid,nlayer,nq,ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,zqtherm,zttherm)
	1338	ELSE
[2114]	1339	zttherm(:,:) = pt(:,:) + pdt(:,:) * ptimestep
	1340	zqtherm(:,:,:) = pq(:,:,:) + pdq(:,:,:) * ptimestep
[2105]	1341	ENDIF
[3232]	1342
[2127]	1343	CALL thermcell_main(ngrid, nlayer, nq, ptimestep, firstcall, &
	1344	pplay, pplev, pphi, zpopsk, &
	1345	pu, pv, zttherm, zqtherm, &
	1346	zdutherm, zdvtherm, zdttherm, zdqtherm, &
[2232]	1347	fm, entr, detr, zw2, fraca)
[3232]	1348
[2060]	1349	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdutherm(1:ngrid,1:nlayer)
	1350	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvtherm(1:ngrid,1:nlayer)
	1351	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdttherm(1:ngrid,1:nlayer)
[3232]	1352
[2127]	1353	IF (tracer) THEN
	1354	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqtherm(1:ngrid,1:nlayer,1:nq)
[2060]	1355	ENDIF
[3232]	1356
[2060]	1357	ENDIF ! end of 'calltherm'
[3232]	1358
[2060]	1359	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1360	! IV.b Dry convective adjustment :
	1361	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[253]	1362
	1363	if(calladj) then
	1364
[1308]	1365	zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
	1366	zduadj(1:ngrid,1:nlayer)=0.0
	1367	zdvadj(1:ngrid,1:nlayer)=0.0
	1368	zdhadj(1:ngrid,1:nlayer)=0.0
[253]	1369
	1370
[1477]	1371	call convadj(ngrid,nlayer,nq,ptimestep, &
	1372	pplay,pplev,zpopsk, &
	1373	pu,pv,zh,pq, &
	1374	pdu,pdv,zdh,pdq, &
	1375	zduadj,zdvadj,zdhadj, &
[2232]	1376	zdqadj)
[253]	1377
[1308]	1378	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
	1379	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer)
	1380	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
	1381	zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
[1283]	1382
[3232]	1383	if(tracer) then
	1384	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
[253]	1385	end if
	1386
[1477]	1387	! Test energy conservation
[253]	1388	if(enertest)then
[1524]	1389	call planetwide_sumval(cppmassarea(:,:)zdtadj(:,:)/totarea_planet,dEtot)
[1295]	1390	if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
[253]	1391	endif
	1392
[1477]	1393	! Test water conservation
[253]	1394	if(watertest)then
[1524]	1395	call planetwide_sumval(massarea(:,:)zdqadj(:,:,igcm_h2o_vap)ptimestep/totarea_planet,dWtot_tmp)
[253]	1396	do ig = 1, ngrid
[1524]	1397	cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
	1398	enddo
	1399	call planetwide_sumval(massarea(:,:)zdqadj(:,:,igcm_h2o_ice)ptimestep/totarea_planet,dWtot)
	1400	dWtot = dWtot + dWtot_tmp
[651]	1401	do ig = 1, ngrid
[1524]	1402	cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
[3232]	1403	enddo
[1524]	1404	call planetwide_maxval(cadjncons(:),nconsMAX)
[253]	1405
[1295]	1406	if (is_master) then
[1524]	1407	print*,'In convadj atmospheric water change =',dWtot,' kg m-2'
[1477]	1408	print*,'In convadj MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
[1524]	1409	endif
[3232]	1410
[1477]	1411	endif ! end of 'watertest'
[3232]	1412
[1477]	1413	endif ! end of 'calladj'
[2299]	1414	!----------------------------------------------
	1415	! Non orographic Gravity Waves:
	1416	!---------------------------------------------
	1417	IF (calllott_nonoro) THEN
	1418
[2595]	1419	CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep, &
	1420	cpnew, rnew, &
	1421	pplay, &
[2299]	1422	zmax_th, &! max altitude reached by thermals (m)
	1423	pt, pu, pv, &
	1424	pdt, pdu, pdv, &
	1425	zustrhi,zvstrhi, &
	1426	d_t_hin, d_u_hin, d_v_hin)
	1427
	1428	! Update tendencies
	1429	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) &
	1430	+ d_t_hin(1:ngrid,1:nlayer)
	1431	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) &
	1432	+ d_u_hin(1:ngrid,1:nlayer)
	1433	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) &
	1434	+ d_v_hin(1:ngrid,1:nlayer)
	1435	print*,'du_nonoro: max ', maxval(d_u_hin), 'min ', minval(d_u_hin)
	1436	print*,'dv_nonoro: max ', maxval(d_v_hin), 'min ', minval(d_v_hin)
	1437
	1438	ENDIF ! of IF (calllott_nonoro)
	1439
	1440
[3232]	1441
[1477]	1442	!-----------------------------------------------
	1443	! V. Carbon dioxide condensation-sublimation :
	1444	!-----------------------------------------------
[253]	1445
	1446	if (co2cond) then
[3232]	1447
[253]	1448	if (.not.tracer) then
	1449	print*,'We need a CO2 ice tracer to condense CO2'
	1450	call abort
	1451	endif
[1477]	1452	call condense_co2(ngrid,nlayer,nq,ptimestep, &
	1453	capcal,pplay,pplev,tsurf,pt, &
[1485]	1454	pdt,zdtsurf,pq,pdq, &
	1455	qsurf,zdqsurfc,albedo,emis, &
[1482]	1456	albedo_bareground,albedo_co2_ice_SPECTV, &
[1485]	1457	zdtc,zdtsurfc,pdpsrf,zdqc)
[253]	1458
[1484]	1459	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
	1460	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
[728]	1461
[1484]	1462	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq)
	1463	dqsurf(1:ngrid,igcm_co2_ice) = dqsurf(1:ngrid,igcm_co2_ice) + zdqsurfc(1:ngrid)
[253]	1464
[3100]	1465
[253]	1466	! test energy conservation
	1467	if(enertest)then
[1524]	1468	call planetwide_sumval(cppmassarea(:,:)zdtc(:,:)/totarea_planet,dEtot)
[1542]	1469	call planetwide_sumval(capcal(:)zdtsurfc(:)cell_area(:)/totarea_planet,dEtots)
[1524]	1470	if (is_master) then
	1471	print*,'In co2cloud atmospheric energy change =',dEtot,' W m-2'
[1477]	1472	print*,'In co2cloud surface energy change =',dEtots,' W m-2'
[1524]	1473	endif
[253]	1474	endif
	1475
[1477]	1476	endif ! end of 'co2cond'
[253]	1477
	1478
[1477]	1479	!---------------------------------------------
[3232]	1480	! VI. Specific parameterizations for tracers
[1477]	1481	!---------------------------------------------
[253]	1482
[1477]	1483	if (tracer) then
[3232]	1484
[3299]	1485	! -----------------------------------------------------
	1486	! VI.0. Volcanoes injecting tracers in the atmosphere
	1487	! -----------------------------------------------------
	1488	if (callvolcano) then
	1489	call volcano(ngrid,nlayer,pplev,pu,pv,pt,zzlev,zdqvolc,nq,massarea,&
	1490	zday,ptimestep)
	1491	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + &
	1492	zdqvolc(1:ngrid,1:nlayer,1:nq)
	1493	endif
	1494
[1477]	1495	! ---------------------
	1496	! VI.1. Water and ice
	1497	! ---------------------
[253]	1498	if (water) then
[3232]	1499
[1477]	1500	! Water ice condensation in the atmosphere
[728]	1501	if(watercond.and.(RLVTT.gt.1.e-8))then
[3232]	1502
[2871]	1503	if ((.not.calltherm).and.moistadjustment) then
[2105]	1504	dqmoist(1:ngrid,1:nlayer,1:nq)=0.
	1505	dtmoist(1:ngrid,1:nlayer)=0.
[3232]	1506
[2105]	1507	! Moist Convective Adjustment.
	1508	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3270]	1509
[2105]	1510	call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
[3232]	1511
[2105]	1512	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) &
	1513	+ dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)
	1514	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) &
	1515	+ dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice)
	1516	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer)
[3232]	1517
[2105]	1518	! Test energy conservation.
	1519	if(enertest)then
	1520	call planetwide_sumval(cppmassarea(:,:)dtmoist(:,:)/totarea_planet,dEtot)
	1521	call planetwide_maxval(dtmoist(:,:),dtmoist_max)
	1522	call planetwide_minval(dtmoist(:,:),dtmoist_min)
	1523	madjdEz(:,:)=cppmass(:,:)dtmoist(:,:)
[3232]	1524
[2105]	1525	do ig=1,ngrid
	1526	madjdE(ig) = cppSUM(mass(:,:)dtmoist(:,:))
	1527	enddo
[3232]	1528
[2105]	1529	if (is_master) then
	1530	print*,'In moistadj atmospheric energy change =',dEtot,' W m-2'
	1531	print,'In moistadj MAX atmospheric energy change =',dtmoist_maxptimestep,'K/step'
	1532	print,'In moistadj MIN atmospheric energy change =',dtmoist_minptimestep,'K/step'
	1533	endif
[3232]	1534
[2105]	1535	call planetwide_sumval(massarea(:,:)dqmoist(:,:,igcm_h2o_vap)ptimestep/totarea_planet+ &
	1536	massarea(:,:)dqmoist(:,:,igcm_h2o_ice)ptimestep/totarea_planet,dWtot)
	1537	if (is_master) print*,'In moistadj atmospheric water change =',dWtot,' kg m-2'
[3232]	1538
[2105]	1539	endif ! end of 'enertest'
[2875]	1540	else
	1541	! rneb_man, dqmoist are output of moistadj, but used later on
	1542	! so set them to zero if moistadj is not called
	1543	rneb_man(:,:)=0
	1544	dqmoist(:,:,:)=0
	1545	endif ! end of '(.not.calltherm).and.moistadjustment'
[3232]	1546
[1477]	1547	! Large scale condensation/evaporation.
	1548	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1308]	1549	call largescale(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
[253]	1550
[1308]	1551	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtlscale(1:ngrid,1:nlayer)
	1552	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)+dqvaplscale(1:ngrid,1:nlayer)
	1553	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice)+dqcldlscale(1:ngrid,1:nlayer)
[253]	1554
[1477]	1555	! Test energy conservation.
[253]	1556	if(enertest)then
[1016]	1557	lscaledEz(:,:) = cppmass(:,:)dtlscale(:,:)
[787]	1558	do ig=1,ngrid
[728]	1559	lscaledE(ig) = cppSUM(mass(:,:)dtlscale(:,:))
[253]	1560	enddo
[1524]	1561	call planetwide_sumval(cppmassarea(:,:)dtlscale(:,:)/totarea_planet,dEtot)
[1477]	1562
[1524]	1563	if (is_master) print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
[728]	1564
[1477]	1565	! Test water conservation.
[1524]	1566	call planetwide_sumval(massarea(:,:)dqvaplscale(:,:)ptimestep/totarea_planet+ &
	1567	SUM(massarea(:,:)dqcldlscale(:,:))ptimestep/totarea_planet,dWtot)
[3232]	1568
[1524]	1569	if (is_master) print*,'In largescale atmospheric water change =',dWtot,' kg m-2'
[1477]	1570	endif ! end of 'enertest'
[253]	1571
[1477]	1572	! Compute cloud fraction.
[253]	1573	do l = 1, nlayer
[787]	1574	do ig=1,ngrid
[3232]	1575	cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
[253]	1576	enddo
	1577	enddo
	1578
[1477]	1579	endif ! end of 'watercond'
[3232]	1580
[1477]	1581	! Water ice / liquid precipitation.
	1582	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3232]	1583	zdqrain(1:ngrid,1:nlayer,1:nq) = 0.0 !JL18 need to do that everytimestep if mass redis is on.
[1989]	1584
[728]	1585	if(waterrain)then
[253]	1586
[787]	1587	zdqsrain(1:ngrid) = 0.0
	1588	zdqssnow(1:ngrid) = 0.0
[253]	1589
[2871]	1590	call rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,pphi,pt,pdt,pq,pdq, &
[1859]	1591	zdtrain,zdqrain,zdqsrain,zdqssnow,reevap_precip,cloudfrac)
[253]	1592
[1308]	1593	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) &
[1524]	1594	+ zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)
[1308]	1595	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) &
[1524]	1596	+ zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice)
[1308]	1597	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain(1:ngrid,1:nlayer)
[3232]	1598
[1477]	1599	dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid)
[3232]	1600	dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
	1601
[1477]	1602	! Test energy conservation.
[651]	1603	if(enertest)then
[3232]	1604
[1524]	1605	call planetwide_sumval(cppmassarea(:,:)zdtrain(:,:)/totarea_planet,dEtot)
	1606	if (is_master) print*,'In rain atmospheric T energy change =',dEtot,' W m-2'
	1607	call planetwide_sumval(massarea(:,:)zdqrain(:,:,igcm_h2o_ice)/totarea_planetRLVTT/cpp,dItot_tmp)
[1542]	1608	call planetwide_sumval(cell_area(:)zdqssnow(:)/totarea_planetRLVTT/cpp,dItot)
[1524]	1609	dItot = dItot + dItot_tmp
	1610	call planetwide_sumval(massarea(:,:)zdqrain(:,:,igcm_h2o_vap)ptimestep/totarea_planet,dVtot_tmp)
[1542]	1611	call planetwide_sumval(cell_area(:)zdqsrain(:)/totarea_planetRLVTT/cpp,dVtot)
[1524]	1612	dVtot = dVtot + dVtot_tmp
	1613	dEtot = dItot + dVtot
[3232]	1614
[1524]	1615	if (is_master) then
[1477]	1616	print*,'In rain dItot =',dItot,' W m-2'
	1617	print*,'In rain dVtot =',dVtot,' W m-2'
	1618	print*,'In rain atmospheric L energy change =',dEtot,' W m-2'
[1524]	1619	endif
[253]	1620
[1477]	1621	! Test water conservation
[1524]	1622	call planetwide_sumval(massarea(:,:)zdqrain(:,:,igcm_h2o_vap)ptimestep/totarea_planet+ &
	1623	massarea(:,:)zdqrain(:,:,igcm_h2o_ice)ptimestep/totarea_planet,dWtot)
[1542]	1624	call planetwide_sumval((zdqsrain(:)+zdqssnow(:))cell_area(:)ptimestep/totarea_planet,dWtots)
[3232]	1625
[1524]	1626	if (is_master) then
	1627	print*,'In rain atmospheric water change =',dWtot,' kg m-2'
	1628	print*,'In rain surface water change =',dWtots,' kg m-2'
	1629	print*,'In rain non-cons factor =',dWtot+dWtots,' kg m-2'
	1630	endif
[3232]	1631
[1477]	1632	endif ! end of 'enertest'
[253]	1633
[1477]	1634	end if ! enf of 'waterrain'
[3232]	1635
[1477]	1636	end if ! end of 'water'
[253]	1637
[1801]	1638	! -------------------------
	1639	! VI.2. Photochemistry
[1477]	1640	! -------------------------
[1801]	1641
[2058]	1642	#ifndef MESOSCALE
[1801]	1643	IF (photochem) then
	1644
	1645	DO ig=1,ngrid
	1646	array_zero1(ig)=0.0
	1647	DO l=1,nlayer
	1648	array_zero2(ig,l)=0.
	1649	ENDDO
	1650	ENDDO
	1651
	1652	call calchim_asis(ngrid,nlayer,nq, &
	1653	ptimestep,pplay,pplev,pt,pdt,dist_star,mu0, &
	1654	fract,zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim, &
	1655	array_zero1,array_zero1, &
[2542]	1656	pu,pdu,pv,pdv,array_zero2,array_zero2,icount,zdtchim)
[1801]	1657
	1658	! increment values of tracers:
[2542]	1659	iesp = 0
	1660	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	1661	! tracers is zero anyways
	1662	! September 2020: flag is_chim to increment only on chemical species
	1663	IF (is_chim(iq)==1) THEN
	1664	iesp = iesp + 1
	1665	DO l=1,nlayer
	1666	DO ig=1,ngrid
	1667	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iesp)
	1668	ENDDO
[1801]	1669	ENDDO
[2542]	1670	ENDIF
[1801]	1671	ENDDO ! of DO iq=1,nq
	1672
	1673
	1674	! increment surface values of tracers:
	1675	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	1676	! tracers is zero anyways
	1677	DO ig=1,ngrid
	1678	! dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
	1679	ENDDO
	1680	ENDDO ! of DO iq=1,nq
	1681
[2542]	1682	! increment values of temperature:
	1683	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtchim(1:ngrid,1:nlayer)
	1684
[1801]	1685	END IF ! of IF (photochem)
[2058]	1686	#endif
[1801]	1687
	1688
[3277]	1689	! ------------------------------------------
	1690	! VI.3.a Generic Condensable Species (GCS)
	1691	! ------------------------------------------
[2701]	1692
[3232]	1693	if (generic_condensation) then
[3277]	1694	if ((.not.water).and.(.not.calltherm).and.moistadjustment_generic) then
	1695	write(,) "moist adjustment for GCS"
	1696	dqmoist(1:ngrid,1:nlayer,1:nq)=0.
	1697	dtmoist(1:ngrid,1:nlayer)=0.
	1698
	1699	! Moist Convective Adjustment.
	1700	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1701	call moistadj_generic(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
	1702
	1703	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) &
	1704	+ dqmoist(1:ngrid,1:nlayer,1:nq)
	1705	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer)
	1706
	1707	! Test energy conservation.
	1708	if(enertest)then
	1709	call planetwide_sumval(cppmassarea(:,:)dtmoist(:,:)/totarea_planet,dEtot)
	1710	call planetwide_maxval(dtmoist(:,:),dtmoist_max)
	1711	call planetwide_minval(dtmoist(:,:),dtmoist_min)
	1712	madjdEz(:,:)=cppmass(:,:)dtmoist(:,:)
	1713
	1714	do ig=1,ngrid
	1715	madjdE(ig) = cppSUM(mass(:,:)dtmoist(:,:))
	1716	enddo
	1717
	1718	if (is_master) then
[3278]	1719	print*,'In moistadj_generic atmospheric energy change =',dEtot,' W m-2'
	1720	print,'In moistadj_generic MAX atmospheric energy change =',dtmoist_maxptimestep,'K/step'
	1721	print,'In moistadj_generic MIN atmospheric energy change =',dtmoist_minptimestep,'K/step'
[3277]	1722	endif
	1723
	1724	! igcm_generic_vap & igcm_generic_ice are not declared below
[3278]	1725	call planetwide_sumval(massarea(:,:)dqmoist(:,:,1)ptimestep/totarea_planet+ &
	1726	massarea(:,:)dqmoist(:,:,2)ptimestep/totarea_planet,dWtot)
	1727	if (is_master) print*,'In moistadj_generic atmospheric GCS change =',dWtot,' kg m-2'
[3277]	1728
	1729	endif ! end of 'enertest'
	1730	else
	1731	! rneb_man, dqmoist are output of moistadj, but used later on
	1732	! so set them to zero if moistadj is not called
	1733	rneb_man(:,:)=0
	1734	dqmoist(:,:,:)=0
	1735	endif ! end of '(.not.calltherm).and.moistadjustment_generic'
	1736
[2701]	1737	call condensation_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay, &
[2890]	1738	pt,pq,pdt,pdq,dt_generic_condensation, &
[2724]	1739	dqvaplscale_generic,dqcldlscale_generic,rneb_generic)
[2890]	1740	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dt_generic_condensation(1:ngrid,1:nlayer)
[2700]	1741	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+dqvaplscale_generic(1:ngrid,1:nlayer,1:nq)
	1742	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+dqcldlscale_generic(1:ngrid,1:nlayer,1:nq)
[2714]	1743
	1744	if(enertest)then
	1745	do ig=1,ngrid
[2890]	1746	genericconddE(ig) = cppSUM(mass(:,:)dt_generic_condensation(:,:))
[2714]	1747	enddo
	1748
[2890]	1749	call planetwide_sumval(cppmassarea(:,:)dt_generic_condensation(:,:)/totarea_planet,dEtot)
[2714]	1750
	1751	if (is_master) print*,'In generic condensation atmospheric energy change =',dEtot,' W m-2'
	1752	end if
[2724]	1753
[3232]	1754	if (.not. water) then
[2726]	1755	! Compute GCS (Generic Condensable Specie) cloud fraction. For now we can not have both water cloud fraction and GCS cloud fraction
	1756	! Water is the priority
	1757	! If you have set water and generic_condensation, then cloudfrac will be water cloudfrac
	1758	!
	1759	! If you have set generic_condensation (and not water) and you have set several GCS
	1760	! then cloudfrac will be only the cloudfrac of the last generic tracer
	1761	! (Because it is rewritten every tracer in the loop)
	1762	!
	1763	! Maybe one should create a cloudfrac_generic(ngrid,nlayer,nq) with 3 dimensions, the last one for tracers
[2724]	1764
	1765	! Let's loop on tracers
[2802]	1766	cloudfrac(:,:)=0.0
[2724]	1767	do iq=1,nq
	1768	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
	1769	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
	1770	do l = 1, nlayer
	1771	do ig=1,ngrid
[2802]	1772	cloudfrac(ig,l)=rneb_generic(ig,l,iq)
[2724]	1773	enddo
	1774	enddo
	1775	endif
[3232]	1776	end do ! do iq=1,nq loop on tracers
[2724]	1777	endif ! .not. water
	1778
[2700]	1779	endif !generic_condensation
[253]	1780
[2721]	1781	!Generic Rain
	1782
	1783	if (generic_rain) then
	1784
	1785	zdqsrain_generic(1:ngrid,1:nq) = 0.0
	1786	zdqssnow_generic(1:ngrid,1:nq) = 0.0
	1787
[2891]	1788	call rain_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay,pphi,pt,pdt,pq,pdq, &
[2721]	1789	zdtrain_generic,dq_rain_generic_vap,dq_rain_generic_cld, &
	1790	zdqsrain_generic,zdqssnow_generic,reevap_precip_generic,cloudfrac)
	1791
	1792	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) &
	1793	+ dq_rain_generic_vap(1:ngrid,1:nlayer,1:nq)
	1794	! only the parts with indexes of generic vapor tracers are filled in dq_rain_generic_vap, other parts are 0.
	1795	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) &
	1796	+ dq_rain_generic_cld(1:ngrid,1:nlayer,1:nq)
	1797	! only the parts with indexes of generic ice(cloud) tracers are filled in dq_rain_generic_cld, other parts are 0.
	1798
	1799	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain_generic(1:ngrid,1:nlayer)
[3232]	1800
[2721]	1801	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq)+zdqsrain_generic(1:ngrid,1:nq)
	1802
	1803	! Test energy conservation.
	1804	if(enertest)then
[3232]	1805
[2721]	1806	call planetwide_sumval(cppmassarea(:,:)zdtrain_generic(:,:)/totarea_planet,dEtot)
	1807	if (is_master) print*,'In rain_generic atmospheric T energy change =',dEtot,' W m-2'
	1808
	1809	! Test conservationfor each generic condensable specie (GCS) tracer
	1810
	1811	! Let's loop on tracers
	1812
	1813	do iq=1,nq
	1814
[2722]	1815	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
[3232]	1816
[2722]	1817	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
[2721]	1818
	1819	call planetwide_sumval(massarea(:,:)dq_rain_generic_cld(:,:,igcm_generic_ice)/totarea_planetRLVTT/cpp,dItot_tmp)
	1820	call planetwide_sumval(cell_area(:)zdqssnow_generic(:,igcm_generic_ice)/totarea_planetRLVTT/cpp,dItot)
	1821	dItot = dItot + dItot_tmp
[3278]	1822	call planetwide_sumval(massarea(:,:)dq_rain_generic_vap(:,:,igcm_generic_vap)/totarea_planetRLVTT/cpp,dVtot_tmp)
[2721]	1823	call planetwide_sumval(cell_area(:)zdqsrain_generic(:,igcm_generic_ice)/totarea_planetRLVTT/cpp,dVtot)
	1824	dVtot = dVtot + dVtot_tmp
	1825	dEtot = dItot + dVtot
[3232]	1826
[2721]	1827	if (is_master) then
	1828	print*,'In rain_generic dItot =',dItot,' W m-2'
	1829	print*,'In rain_generic dVtot =',dVtot,' W m-2'
	1830	print*,'In rain_generic atmospheric L energy change =',dEtot,' W m-2'
	1831	endif
	1832
	1833	call planetwide_sumval(massarea(:,:)dq_rain_generic_vap(:,:,igcm_generic_vap)ptimestep/totarea_planet+ &
	1834	massarea(:,:)dq_rain_generic_cld(:,:,igcm_generic_ice)ptimestep/totarea_planet,dWtot)
	1835	call planetwide_sumval((zdqsrain_generic(:,igcm_generic_ice)+zdqssnow_generic(:,igcm_generic_ice))cell_area(:)ptimestep/totarea_planet,dWtots)
[3232]	1836
[2721]	1837	if (is_master) then
	1838	print*,'In rain_generic atmospheric generic tracer change =',dWtot,' kg m-2'
	1839	print*,'In rain_generic surface generic tracer change =',dWtots,' kg m-2'
	1840	print*,'In rain_generic non-cons factor =',dWtot+dWtots,' kg m-2'
	1841	endif
	1842
	1843	endif
	1844
[3232]	1845	end do ! do iq=1,nq loop on tracers
	1846
[2721]	1847	endif ! end of 'enertest'
[3232]	1848
[2721]	1849	endif !generic_rain
	1850
[3277]	1851	! -------------------------
	1852	! VI.3.b. Sedimentation.
	1853	! -------------------------
[1477]	1854	if (sedimentation) then
[2802]	1855
[1477]	1856	zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
	1857	zdqssed(1:ngrid,1:nq) = 0.0
[253]	1858
[1477]	1859	if(watertest)then
[3232]	1860
[1477]	1861	iq=igcm_h2o_ice
[1524]	1862	call planetwide_sumval(massarea(:,:)pq(:,:,iq)ptimestep/totarea_planet,dWtot)
	1863	call planetwide_sumval(massarea(:,:)pdq(:,:,iq)ptimestep/totarea_planet,dWtots)
	1864	if (is_master) then
	1865	print*,'Before sedim pq =',dWtot,' kg m-2'
[1477]	1866	print*,'Before sedim pdq =',dWtots,' kg m-2'
[1524]	1867	endif
[1477]	1868	endif
[3232]	1869
[1477]	1870	call callsedim(ngrid,nlayer,ptimestep, &
	1871	pplev,zzlev,pt,pdt,pq,pdq, &
	1872	zdqsed,zdqssed,nq)
[253]	1873
[1477]	1874	if(watertest)then
	1875	iq=igcm_h2o_ice
[1524]	1876	call planetwide_sumval(massarea(:,:)pq(:,:,iq)ptimestep/totarea_planet,dWtot)
	1877	call planetwide_sumval(massarea(:,:)pdq(:,:,iq)ptimestep/totarea_planet,dWtots)
	1878	if (is_master) then
	1879	print*,'After sedim pq =',dWtot,' kg m-2'
	1880	print*,'After sedim pdq =',dWtots,' kg m-2'
	1881	endif
[1477]	1882	endif
[253]	1883
[1477]	1884	! Whether it falls as rain or snow depends only on the surface temperature
	1885	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
	1886	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
[253]	1887
[1477]	1888	! Test water conservation
	1889	if(watertest)then
[1524]	1890	call planetwide_sumval(massarea(:,:)(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice))ptimestep/totarea_planet,dWtot)
[1542]	1891	call planetwide_sumval((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))cell_area(:)ptimestep/totarea_planet,dWtots)
[1524]	1892	if (is_master) then
	1893	print*,'In sedim atmospheric ice change =',dWtot,' kg m-2'
	1894	print*,'In sedim surface ice change =',dWtots,' kg m-2'
	1895	print*,'In sedim non-cons factor =',dWtot+dWtots,' kg m-2'
	1896	endif
[1477]	1897	endif
[253]	1898
[1477]	1899	end if ! end of 'sedimentation'
[253]	1900
	1901
[1477]	1902	! ---------------
[1801]	1903	! VI.4. Updates
[1477]	1904	! ---------------
[253]	1905
[1477]	1906	! Updating Atmospheric Mass and Tracers budgets.
[728]	1907	if(mass_redistrib) then
	1908
[1477]	1909	zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * &
[1524]	1910	( zdqevap(1:ngrid,1:nlayer) &
	1911	+ zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap) &
	1912	+ dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap) &
	1913	+ dqvaplscale(1:ngrid,1:nlayer) )
[863]	1914
	1915	do ig = 1, ngrid
[1524]	1916	zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
[863]	1917	enddo
[3232]	1918
[1524]	1919	call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
	1920	call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
	1921	call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
[728]	1922
[1524]	1923	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
[1477]	1924	rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
[1524]	1925	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
	1926	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
[3232]	1927
[1308]	1928	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
[1477]	1929	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
	1930	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
	1931	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
	1932	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
[1524]	1933	pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
[1477]	1934	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
[3232]	1935
[1524]	1936	endif
[728]	1937
[3100]	1938
[1477]	1939	! ------------------
[1801]	1940	! VI.5. Slab Ocean
[1477]	1941	! ------------------
[728]	1942
[1477]	1943	if (ok_slab_ocean)then
[728]	1944
[1477]	1945	do ig=1,ngrid
	1946	qsurfint(:,igcm_h2o_ice)=qsurf(:,igcm_h2o_ice)
	1947	enddo
[1297]	1948
[3268]	1949	!! These two commented functions are from BC2014 (only kept for record)
	1950	!!
[3100]	1951	!! call ocean_slab_ice(ptimestep, &
	1952	!! ngrid, knindex, tsea_ice, fluxrad, &
	1953	!! zdqssnow, qsurf(:,igcm_h2o_ice), &
	1954	!! - zdqsdif(:,igcm_h2o_vap), &
	1955	!! flux_sens_lat,tsea_ice, pctsrf_sic, &
	1956	!! taux,tauy,icount)
	1957	!! call ocean_slab_get_vars(ngrid,tslab, &
	1958	!! sea_ice, flux_o, &
	1959	!! flux_g, dt_hdiff, &
	1960	!! dt_ekman)
[3232]	1961
[3100]	1962	call ocean_slab_noice(icount, ptimestep, knon, knindex, &
	1963	zdqssnow, tsea_ice, &
	1964	fluxrad, qsurf(:,igcm_h2o_ice), flux_sens_lat, &
	1965	tsea_ice, taux, tauy, zmasq)
	1966
	1967	call ocean_slab_ice(icount, ptimestep, knon, knindex, &
	1968	zdqssnow, tsea_ice, &
	1969	fluxrad, qsurf(:,igcm_h2o_ice), flux_sens_lat, &
	1970	tsea_ice, -zdqsdif(:,igcm_h2o_vap), taux, tauy, zmasq)
	1971
	1972	call ocean_slab_frac(pctsrf_sic, zmasq)
	1973
	1974	call ocean_slab_get_vars(ngrid, tslab, sea_ice, flux_g, &
	1975	dt_hdiff, dt_ekman)
[3268]	1976
	1977	!! sea_ice defines the sea ice thickness in kg/m2
	1978	!! pctsrf_sic defines the percentage of the oceanic grid that is covered by sea ice (between 0 and almost 1)
	1979	!! qsurf(:,igcm_h2o_ice) defines the amount of snow that accumulates on top of the sea ice (in kg/m2) ; note that the snow fraction depends on the amount of snow.
	1980	!! therefore, beware that qsurf(:,igcm_h2o_ice) has its own meaning for oceanic grid points (when ok_slab_ocean activated)
[3232]	1981
[3268]	1982
[1297]	1983	do ig=1,ngrid
	1984	if (nint(rnat(ig)).eq.1)then
[1477]	1985	tslab(ig,1) = 0.
	1986	tslab(ig,2) = 0.
	1987	tsea_ice(ig) = 0.
[3100]	1988	! tice(ig) = 0.
[1477]	1989	sea_ice(ig) = 0.
	1990	pctsrf_sic(ig) = 0.
	1991	qsurf(ig,igcm_h2o_ice) = qsurfint(ig,igcm_h2o_ice)
[1297]	1992	endif
	1993	enddo
	1994
[1477]	1995	endif ! end of 'ok_slab_ocean'
[1297]	1996
[3100]	1997
[1477]	1998	! -----------------------------
[1801]	1999	! VI.6. Surface Tracer Update
[1477]	2000	! -----------------------------
[1297]	2001
[253]	2002	if(hydrology)then
[3232]	2003
[1482]	2004	call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, &
	2005	capcal,albedo,albedo_bareground, &
[1524]	2006	albedo_snow_SPECTV,albedo_co2_ice_SPECTV, &
[1482]	2007	mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic, &
[1477]	2008	sea_ice)
[3232]	2009
[1484]	2010	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurf_hyd(1:ngrid)
	2011	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + dqs_hyd(1:ngrid,1:nq)
[3100]	2012
[3232]	2013
[1484]	2014	qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
[253]	2015
[1477]	2016	! Test energy conservation
[253]	2017	if(enertest)then
[1542]	2018	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf_hyd(:)/totarea_planet,dEtots)
[1524]	2019	if (is_master) print*,'In hydrol surface energy change =',dEtots,' W m-2'
[253]	2020	endif
	2021
	2022	! test water conservation
	2023	if(watertest)then
[1542]	2024	call planetwide_sumval(dqs_hyd(:,igcm_h2o_ice)cell_area(:)ptimestep/totarea_planet,dWtots)
[1524]	2025	if (is_master) print*,'In hydrol surface ice change =',dWtots,' kg m-2'
[1542]	2026	call planetwide_sumval(dqs_hyd(:,igcm_h2o_vap)cell_area(:)ptimestep/totarea_planet,dWtots)
[1524]	2027	if (is_master) then
[1477]	2028	print*,'In hydrol surface water change =',dWtots,' kg m-2'
	2029	print*,'---------------------------------------------------------------'
[1524]	2030	endif
[253]	2031	endif
	2032
[1477]	2033	else ! of if (hydrology)
[253]	2034
[1484]	2035	qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
[253]	2036
[1477]	2037	end if ! of if (hydrology)
[253]	2038
[3232]	2039	! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
[1477]	2040	! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
[622]	2041	qsurf_hist(:,:) = qsurf(:,:)
[253]	2042
[1477]	2043	endif! end of if 'tracer'
[253]	2044
	2045
[1477]	2046	!------------------------------------------------
[3232]	2047	! VII. Surface and sub-surface soil temperature
[1477]	2048	!------------------------------------------------
[253]	2049
[1477]	2050
	2051	! Increment surface temperature
[1297]	2052	if(ok_slab_ocean)then
	2053	do ig=1,ngrid
	2054	if (nint(rnat(ig)).eq.0)then
	2055	zdtsurf(ig)= (tslab(ig,1) &
	2056	+ pctsrf_sic(ig)*(tsea_ice(ig)-tslab(ig,1))-tsurf(ig))/ptimestep
	2057	endif
	2058	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
	2059	enddo
[3232]	2060
[1297]	2061	else
[3232]	2062	tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
[1477]	2063	endif ! end of 'ok_slab_ocean'
[1297]	2064
[1477]	2065
	2066	! Compute soil temperatures and subsurface heat flux.
[253]	2067	if (callsoil) then
[787]	2068	call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, &
[3232]	2069	ptimestep,tsurf,tsoil,capcal,fluxgrd)
[253]	2070	endif
	2071
[1297]	2072
	2073	if (ok_slab_ocean) then
[3232]	2074
[1477]	2075	do ig=1,ngrid
[3232]	2076
[1477]	2077	fluxgrdocean(ig)=fluxgrd(ig)
	2078	if (nint(rnat(ig)).eq.0) then
[1297]	2079	capcal(ig)=capcalocean
	2080	fluxgrd(ig)=0.
[3100]	2081	!! fluxgrdocean(ig)=pctsrf_sic(ig)flux_g(ig)+(1-pctsrf_sic(ig))(dt_hdiff(ig,1)+dt_ekman(ig,1))
	2082	! Dividing by cell area to have flux in W/m2
	2083	fluxgrdocean(ig)=flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1))/cell_area(ig)
[1477]	2084	do iq=1,nsoilmx
	2085	tsoil(ig,iq)=tsurf(ig)
	2086	enddo
	2087	if (pctsrf_sic(ig).gt.0.5) then
	2088	capcal(ig)=capcalseaice
	2089	if (qsurf(ig,igcm_h2o_ice).gt.0.) then
	2090	capcal(ig)=capcalsno
	2091	endif
[3232]	2092	endif
[1477]	2093	endif
[3232]	2094
[1477]	2095	enddo
[3232]	2096
[1477]	2097	endif !end of 'ok_slab_ocean'
[1297]	2098
[1477]	2099
	2100	! Test energy conservation
[253]	2101	if(enertest)then
[3232]	2102	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf(:)/totarea_planet,dEtots)
[1524]	2103	if (is_master) print*,'Surface energy change =',dEtots,' W m-2'
[253]	2104	endif
	2105
	2106
[1477]	2107	!---------------------------------------------------
	2108	! VIII. Perform diagnostics and write output files
	2109	!---------------------------------------------------
	2110
	2111	! Note : For output only: the actual model integration is performed in the dynamics.
	2112
	2113
[3232]	2114
[1477]	2115	! Temperature, zonal and meridional winds.
[1308]	2116	zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
	2117	zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
	2118	zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
[3232]	2119
[2176]	2120	! Recast thermal plume vertical velocity array for outputs
	2121	IF (calltherm) THEN
	2122	DO ig=1,ngrid
	2123	DO l=1,nlayer
	2124	zw2_bis(ig,l) = zw2(ig,l)
[2232]	2125	fm_bis(ig,l) = fm(ig,l)
[2176]	2126	ENDDO
	2127	ENDDO
	2128	ENDIF
[253]	2129
[1477]	2130	! Diagnostic.
[1637]	2131	zdtdyn(1:ngrid,1:nlayer) = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep
[1308]	2132	ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
[253]	2133
[1637]	2134	zdudyn(1:ngrid,1:nlayer) = (pu(1:ngrid,1:nlayer)-zuprevious(1:ngrid,1:nlayer)) / ptimestep
	2135	zuprevious(1:ngrid,1:nlayer) = zu(1:ngrid,1:nlayer)
	2136
[253]	2137	if(firstcall)then
[1308]	2138	zdtdyn(1:ngrid,1:nlayer)=0.0
[1637]	2139	zdudyn(1:ngrid,1:nlayer)=0.0
[253]	2140	endif
	2141
[1477]	2142	! Dynamical heating diagnostic.
[253]	2143	do ig=1,ngrid
[1637]	2144	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp
[253]	2145	enddo
	2146
[1477]	2147	! Tracers.
[1308]	2148	zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
[253]	2149
[1477]	2150	! Surface pressure.
[787]	2151	ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
[253]	2152
	2153
	2154
[1477]	2155	! Surface and soil temperature information
[1542]	2156	call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
[1295]	2157	call planetwide_minval(tsurf(:),Ts2)
	2158	call planetwide_maxval(tsurf(:),Ts3)
[253]	2159	if(callsoil)then
[1542]	2160	TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
[1699]	2161	if (is_master) then
	2162	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
	2163	print*,Ts1,Ts2,Ts3,TsS
	2164	end if
[959]	2165	else
[1699]	2166	if (is_master) then
	2167	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
[1477]	2168	print*,Ts1,Ts2,Ts3
[1524]	2169	endif
[959]	2170	end if
[253]	2171
	2172
[1477]	2173	! Check the energy balance of the simulation during the run
[253]	2174	if(corrk)then
	2175
[1542]	2176	call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
	2177	call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
	2178	call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
	2179	call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
	2180	call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
[787]	2181	do ig=1,ngrid
[253]	2182	if(fluxtop_dn(ig).lt.0.0)then
	2183	print*,'fluxtop_dn has gone crazy'
	2184	print*,'fluxtop_dn=',fluxtop_dn(ig)
	2185	print*,'tau_col=',tau_col(ig)
	2186	print*,'aerosol=',aerosol(ig,:,:)
	2187	print*,'temp= ',pt(ig,:)
	2188	print*,'pplay= ',pplay(ig,:)
	2189	call abort
	2190	endif
	2191	end do
[3232]	2192
[787]	2193	if(ngrid.eq.1)then
[253]	2194	DYN=0.0
	2195	endif
[3232]	2196
[1524]	2197	if (is_master) then
[1477]	2198	print*,' ISR ASR OLR GND DYN [W m^-2]'
	2199	print*, ISR,ASR,OLR,GND,DYN
[1524]	2200	endif
[253]	2201
[1295]	2202	if(enertest .and. is_master)then
[651]	2203	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
	2204	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
	2205	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
[253]	2206	endif
	2207
[1295]	2208	if(meanOLR .and. is_master)then
[1216]	2209	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
[253]	2210	! to record global radiative balance
[588]	2211	open(92,file="rad_bal.out",form='formatted',position='append')
[651]	2212	write(92,*) zday,ISR,ASR,OLR
[253]	2213	close(92)
[588]	2214	open(93,file="tem_bal.out",form='formatted',position='append')
[1295]	2215	if(callsoil)then
[1524]	2216	write(93,*) zday,Ts1,Ts2,Ts3,TsS
	2217	else
	2218	write(93,*) zday,Ts1,Ts2,Ts3
	2219	endif
[253]	2220	close(93)
	2221	endif
	2222	endif
	2223
[1477]	2224	endif ! end of 'corrk'
[253]	2225
[651]	2226
[1477]	2227	! Diagnostic to test radiative-convective timescales in code.
[253]	2228	if(testradtimes)then
[588]	2229	open(38,file="tau_phys.out",form='formatted',position='append')
[253]	2230	ig=1
	2231	do l=1,nlayer
	2232	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
	2233	enddo
	2234	close(38)
[726]	2235	print*,'As testradtimes enabled,'
	2236	print*,'exiting physics on first call'
[253]	2237	call abort
	2238	endif
	2239
[1477]	2240
	2241	! Compute column amounts (kg m-2) if tracers are enabled.
[253]	2242	if(tracer)then
[787]	2243	qcol(1:ngrid,1:nq)=0.0
[253]	2244	do iq=1,nq
[1477]	2245	do ig=1,ngrid
	2246	qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
	2247	enddo
[253]	2248	enddo
	2249
[1477]	2250	! Generalised for arbitrary aerosols now. By LK
[787]	2251	reffcol(1:ngrid,1:naerkind)=0.0
[728]	2252	if(co2cond.and.(iaero_co2.ne.0))then
[1308]	2253	call co2_reffrad(ngrid,nlayer,nq,zq,reffrad(1,1,iaero_co2))
[787]	2254	do ig=1,ngrid
[1308]	2255	reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayer,igcm_co2_ice)reffrad(ig,1:nlayer,iaero_co2)mass(ig,1:nlayer))
[253]	2256	enddo
[728]	2257	endif
	2258	if(water.and.(iaero_h2o.ne.0))then
[1308]	2259	call h2o_reffrad(ngrid,nlayer,zq(1,1,igcm_h2o_ice),zt, &
[858]	2260	reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
[787]	2261	do ig=1,ngrid
[1308]	2262	reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayer,igcm_h2o_ice)reffrad(ig,1:nlayer,iaero_h2o)mass(ig,1:nlayer))
[728]	2263	enddo
	2264	endif
[253]	2265
[1477]	2266	endif ! end of 'tracer'
[253]	2267
	2268
[1477]	2269	! Test for water conservation.
[253]	2270	if(water)then
	2271
[1542]	2272	call planetwide_sumval(cell_area(:)*qsurf_hist(:,igcm_h2o_ice)/totarea_planet,icesrf)
	2273	call planetwide_sumval(cell_area(:)*qsurf_hist(:,igcm_h2o_vap)/totarea_planet,liqsrf)
	2274	call planetwide_sumval(cell_area(:)*qcol(:,igcm_h2o_ice)/totarea_planet,icecol)
	2275	call planetwide_sumval(cell_area(:)*qcol(:,igcm_h2o_vap)/totarea_planet,vapcol)
[253]	2276
[651]	2277	h2otot = icesrf + liqsrf + icecol + vapcol
[3232]	2278
[1524]	2279	if (is_master) then
[1477]	2280	print*,' Total water amount [kg m^-2]: ',h2otot
	2281	print*,' Surface ice Surface liq. Atmos. con. Atmos. vap. [kg m^-2] '
	2282	print*, icesrf,liqsrf,icecol,vapcol
[1524]	2283	endif
[253]	2284
[1295]	2285	if(meanOLR .and. is_master)then
[1216]	2286	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
[253]	2287	! to record global water balance
[588]	2288	open(98,file="h2o_bal.out",form='formatted',position='append')
[651]	2289	write(98,*) zday,icesrf,liqsrf,icecol,vapcol
[253]	2290	close(98)
	2291	endif
	2292	endif
	2293
	2294	endif
	2295
	2296
[1477]	2297	! Calculate RH (Relative Humidity) for diagnostic.
[253]	2298	if(water)then
[3232]	2299
[253]	2300	do l = 1, nlayer
[787]	2301	do ig=1,ngrid
[728]	2302	call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l))
[253]	2303	RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
	2304	enddo
	2305	enddo
	2306
[1477]	2307	! Compute maximum possible H2O column amount (100% saturation).
[253]	2308	do ig=1,ngrid
[1477]	2309	H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
[253]	2310	enddo
	2311
[1477]	2312	endif ! end of 'water'
[253]	2313
[2724]	2314	! Calculate RH_generic (Generic Relative Humidity) for diagnostic.
	2315	if(generic_condensation)then
[2802]	2316	RH_generic(:,:,:)=0.0
[2724]	2317	do iq=1,nq
[996]	2318
[2724]	2319	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
[3232]	2320
[2724]	2321	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
[3232]	2322
[2724]	2323	do l = 1, nlayer
	2324	do ig=1,ngrid
	2325	call Psat_generic(zt(ig,l),pplay(ig,l),metallicity,psat_tmp_generic,qsat_generic(ig,l,iq))
	2326	RH_generic(ig,l,iq) = zq(ig,l,igcm_generic_vap) / qsat_generic(ig,l,iq)
	2327	enddo
	2328	enddo
	2329
	2330	end if
[3232]	2331
[2724]	2332	end do ! iq=1,nq
	2333
	2334	endif ! end of 'generic_condensation'
	2335
	2336
[1699]	2337	if (is_master) print,'--> Ls =',zls180./pi
[3232]	2338
	2339
[3235]	2340	! -------------------------
	2341	! IX. Thermosphere
	2342	! -------------------------
	2343	if (callthermos) then
	2344	call conduction(ngrid,nlayer,nq,ptimestep,pplay,pplev, &
	2345	pt,tsurf,zzlev,zzlay,muvar,pq,firstcall,zdtconduc)
	2346	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtconduc(1:ngrid,1:nlayer)
	2347	endif ! of if (callthermos)
	2348
	2349
[1477]	2350	!----------------------------------------------------------------------
[253]	2351	! Writing NetCDF file "RESTARTFI" at the end of the run
[1477]	2352	!----------------------------------------------------------------------
	2353
[253]	2354	! Note: 'restartfi' is stored just before dynamics are stored
	2355	! in 'restart'. Between now and the writting of 'restart',
	2356	! there will have been the itau=itau+1 instruction and
	2357	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	2358	! thus we store for time=time+dtvr
	2359
	2360
	2361
[1477]	2362	if(lastcall) then
	2363	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
[305]	2364
[1836]	2365	#ifndef MESOSCALE
[3232]	2366
[1477]	2367	if (ngrid.ne.1) then
	2368	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
[3232]	2369
[1477]	2370	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
	2371	ptimestep,ztime_fin, &
	2372	tsurf,tsoil,emis,q2,qsurf_hist, &
	2373	cloudfrac,totcloudfrac,hice, &
	2374	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
	2375	endif
[1836]	2376	#endif
[1477]	2377	if(ok_slab_ocean) then
	2378	call ocean_slab_final!(tslab, seaice)
	2379	end if
[1297]	2380
[1682]	2381	endif ! end of 'lastcall'
[253]	2382
[861]	2383
[2958]	2384	! -----------------------------------------------------------------
	2385	! WSTATS: Saving statistics
	2386	! -----------------------------------------------------------------
	2387	! ("stats" stores and accumulates key variables in file "stats.nc"
	2388	! which can later be used to make the statistic files of the run:
	2389	! if flag "callstats" from callphys.def is .true.)
[253]	2390
[3232]	2391
[1477]	2392	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
	2393	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	2394	call wstats(ngrid,"fluxsurf_lw", &
	2395	"Thermal IR radiative flux to surface","W.m-2",2, &
	2396	fluxsurf_lw)
	2397	call wstats(ngrid,"fluxtop_lw", &
	2398	"Thermal IR radiative flux to space","W.m-2",2, &
	2399	fluxtop_lw)
[3232]	2400
[253]	2401	! call wstats(ngrid,"fluxsurf_sw", &
	2402	! "Solar radiative flux to surface","W.m-2",2, &
[3232]	2403	! fluxsurf_sw_tot)
[253]	2404	! call wstats(ngrid,"fluxtop_sw", &
	2405	! "Solar radiative flux to space","W.m-2",2, &
	2406	! fluxtop_sw_tot)
[526]	2407
[253]	2408
[1477]	2409	call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
	2410	call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
	2411	call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
[1482]	2412	!call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
	2413	!call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
[1477]	2414	call wstats(ngrid,"p","Pressure","Pa",3,pplay)
	2415	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
	2416	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
	2417	call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
	2418	call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
	2419	call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
[526]	2420
[1477]	2421	if (tracer) then
	2422	do iq=1,nq
	2423	call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
[3232]	2424	call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
[1477]	2425	'kg m^-2',2,qsurf(1,iq) )
[3232]	2426	call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
[526]	2427	'kg m^-2',2,qcol(1,iq) )
[3232]	2428
	2429	! call wstats(ngrid,trim(noms(iq))//'_reff', &
	2430	! trim(noms(iq))//'_reff', &
[726]	2431	! 'm',3,reffrad(1,1,iq))
[1477]	2432
	2433	end do
[3232]	2434
[253]	2435	if (water) then
[1308]	2436	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
[1477]	2437	call wstats(ngrid,"vmr_h2ovapor", &
[3232]	2438	"H2O vapour volume mixing ratio","mol/mol", &
[1477]	2439	3,vmr)
	2440	endif
[253]	2441
[3232]	2442	endif ! end of 'tracer'
[253]	2443
[1477]	2444	if(watercond.and.CLFvarying)then
	2445	call wstats(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
	2446	call wstats(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
	2447	call wstats(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
	2448	call wstats(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
	2449	call wstats(ngrid,"RH","relative humidity"," ",3,RH)
	2450	endif
[1297]	2451
[1477]	2452	if (ok_slab_ocean) then
[1297]	2453	call wstats(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
	2454	call wstats(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
	2455	call wstats(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1))
	2456	call wstats(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2))
	2457	call wstats(ngrid,"tslab1","tslab1","K",2,tslab(:,1))
	2458	call wstats(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
	2459	call wstats(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
[3100]	2460	call wstats(ngrid,"tsea_ice","top layer temp, snow/ice","K",2,tsea_ice)
	2461	! call wstats(ngrid,"tice","ice temp if snow on top","K",2,tice)
[1297]	2462	call wstats(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
	2463	call wstats(ngrid,"rnat","nature of the surface","",2,rnat)
[1477]	2464	endif
[1297]	2465
[2958]	2466	if(lastcall.and.callstats) then
[1477]	2467	write (,) "Writing stats..."
	2468	call mkstats(ierr)
	2469	endif
[253]	2470
[3232]	2471
[1836]	2472	#ifndef MESOSCALE
[3232]	2473
[1477]	2474	!-----------------------------------------------------------------------------------------------------
	2475	! OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
	2476	!
	2477	! Note 1 : output with period "ecritphy", set in "run.def"
	2478	!
	2479	! Note 2 : writediagfi can also be called from any other subroutine for any variable,
	2480	! but its preferable to keep all the calls in one place ...
	2481	!-----------------------------------------------------------------------------------------------------
[253]	2482
[1477]	2483	call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
	2484	call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
	2485	call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
	2486	call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
	2487	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	2488	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
	2489	call writediagfi(ngrid,"temp","temperature","K",3,zt)
	2490	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
	2491	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
	2492	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
	2493	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
	2494	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
	2495
[965]	2496	! Subsurface temperatures
[969]	2497	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	2498	! call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
[965]	2499
[1477]	2500	! Oceanic layers
	2501	if(ok_slab_ocean) then
	2502	call writediagfi(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
[3100]	2503	call writediagfi(ngrid,"tsea_ice","top layer temp, snow/ice","K",2,tsea_ice)
	2504	! call writediagfi(ngrid,"tice","ice temp if snow on top","K",2,tice)
[1477]	2505	call writediagfi(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
	2506	call writediagfi(ngrid,"tslab1","tslab1","K",2,tslab(:,1))
	2507	call writediagfi(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
	2508	call writediagfi(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
	2509	call writediagfi(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
	2510	call writediagfi(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1))
	2511	call writediagfi(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2))
	2512	call writediagfi(ngrid,"rnat","nature of the surface","",2,rnat)
	2513	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
	2514	call writediagfi(ngrid,"latentFlux","latent heat flux","w.m^-2",2,zdqsdif(:,igcm_h2o_vap)*RLVTT)
	2515	endif
[3232]	2516
[2060]	2517	! Thermal plume model
	2518	if (calltherm) then
[2232]	2519	call writediagfi(ngrid,'entr','Entrainment','kg m$^{-2}$ s$^{-1}$', 3, entr)
	2520	call writediagfi(ngrid,'detr','Detrainment','kg m$^{-2}$ s$^{-1}$', 3, detr)
	2521	call writediagfi(ngrid,'fm','Mass flux','kg m$^{-2}$ s$^{-1}$', 3, fm_bis)
[2176]	2522	call writediagfi(ngrid,'w_plm','Squared vertical velocity','m s$^{-1}$', 3, zw2_bis)
[2060]	2523	call writediagfi(ngrid,'fraca','Updraft fraction','', 3, fraca)
	2524	endif
[2299]	2525
	2526	! GW non-oro outputs
	2527	if (calllott_nonoro) then
[2595]	2528	call WRITEDIAGFI(ngrid,"dugwno","GW non-oro dU","m/s2", 3,d_u_hin)
	2529	call WRITEDIAGFI(ngrid,"dvgwno","GW non-oro dV","m/s2", 3,d_v_hin)
[2299]	2530	endif
[3232]	2531
[1477]	2532	! Total energy balance diagnostics
	2533	if(callrad.and.(.not.newtonian))then
[3232]	2534
[1482]	2535	!call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
	2536	!call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
[1477]	2537	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
	2538	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
	2539	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
	2540	call writediagfi(ngrid,"shad","rings"," ", 2, fract)
[3232]	2541
[1016]	2542	! call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
	2543	! call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
	2544	! call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
	2545	! call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
	2546	! call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
	2547	! call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
[3269]	2548	! call writediagfi(ngrid,"netfluxsurflw","lw net surface flux","W m-2",2,net_fluxsurf_lw)
[1477]	2549
	2550	if(ok_slab_ocean) then
	2551	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrdocean)
	2552	else
	2553	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
	2554	endif
[3232]	2555
[1477]	2556	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
[3232]	2557
[1477]	2558	endif ! end of 'callrad'
[3232]	2559
[1477]	2560	if(enertest) then
[3232]	2561
[1524]	2562	if (calldifv) then
[3232]	2563
[1524]	2564	call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
[1477]	2565	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
[3232]	2566
[1524]	2567	! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
	2568	! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
	2569	! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
[3232]	2570
[1524]	2571	endif
[3232]	2572
[1524]	2573	if (corrk) then
	2574	call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
	2575	call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
	2576	endif
[3232]	2577
[1477]	2578	if(watercond) then
	2579
[3232]	2580	call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
[2871]	2581	if ((.not.calltherm).and.moistadjustment) then
	2582	call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
	2583	endif
[1524]	2584	call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
[3232]	2585
[3270]	2586	! call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz)
	2587	! call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)
[1524]	2588	! call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
[253]	2589
[1477]	2590	endif
[2714]	2591
	2592	if (generic_condensation) then
	2593
	2594	call writediagfi(ngrid,"genericconddE","heat from generic condensation","W m-2",2,genericconddE)
[2890]	2595	call writediagfi(ngrid,"dt_generic_condensation","heating from generic condensation","K s-1",3,dt_generic_condensation)
[3232]	2596
[2714]	2597	endif
[3232]	2598
[1477]	2599	endif ! end of 'enertest'
[253]	2600
[2133]	2601	! Diagnostics of optical thickness
[2138]	2602	! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
[3232]	2603	if (diagdtau) then
[2133]	2604	do nw=1,L_NSPECTV
	2605	write(str2,'(i2.2)') nw
[2138]	2606	call writediagfi(ngrid,'dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',1,int_dtauv(:,nlayer:1:-1,nw))
[2133]	2607	enddo
	2608	do nw=1,L_NSPECTI
	2609	write(str2,'(i2.2)') nw
[2138]	2610	call writediagfi(ngrid,'dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',1,int_dtaui(:,nlayer:1:-1,nw))
[2133]	2611	enddo
	2612	endif
[1477]	2613
[2133]	2614
[1477]	2615	! Temporary inclusions for heating diagnostics.
[2831]	2616	call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
	2617	call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
[2299]	2618	call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
[2831]	2619	call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
[3232]	2620
[3284]	2621	if (tracer .or. water) call writediagfi(ngrid,"dtmoistadj","moist adj heating","K s-1",3,dtmoist)
	2622	if (calladj) call writediagfi(ngrid,"dtdryadj","dry adj heating","K s-1",3,zdtadj)
	2623
[1477]	2624	! For Debugging.
[368]	2625	!call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
[253]	2626	!call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
	2627
[3232]	2628
[1477]	2629	! Output aerosols.
	2630	if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
[1524]	2631	call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
[1477]	2632	if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
[1524]	2633	call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
[1477]	2634	if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
[1524]	2635	call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
[1477]	2636	if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
[1524]	2637	call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
[253]	2638
[1477]	2639	! Output tracers.
	2640	if (tracer) then
[3232]	2641
[1477]	2642	do iq=1,nq
	2643	call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
[3232]	2644	call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
[1477]	2645	'kg m^-2',2,qsurf_hist(1,iq) )
[3232]	2646	call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
[1477]	2647	'kg m^-2',2,qcol(1,iq) )
[3232]	2648	! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
	2649	! 'kg m^-2',2,qsurf(1,iq) )
[253]	2650
[1477]	2651	if(watercond.or.CLFvarying)then
	2652	call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
	2653	call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
	2654	call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
	2655	call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
[1524]	2656	call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
[1477]	2657	endif
[253]	2658
[1477]	2659	if(waterrain)then
	2660	call writediagfi(ngrid,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
	2661	call writediagfi(ngrid,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
[1859]	2662	call writediagfi(ngrid,"reevap","reevaporation of precipitation","kg m-2 s-1",2,reevap_precip)
[1477]	2663	endif
[253]	2664
[2724]	2665	if(generic_condensation)then
	2666	call writediagfi(ngrid,"rneb_generic","GCS cloud fraction (generic condensation)"," ",3,rneb_generic)
	2667	call writediagfi(ngrid,"CLF","GCS cloud fraction"," ",3,cloudfrac)
	2668	call writediagfi(ngrid,"RH_generic","GCS relative humidity"," ",3,RH_generic)
	2669	endif
	2670
[2721]	2671	if(generic_rain)then
	2672	call writediagfi(ngrid,"rain","generic rainfall","kg m-2 s-1",2,zdqsrain_generic)
	2673	call writediagfi(ngrid,"snow","generic snowfall","kg m-2 s-1",2,zdqssnow_generic)
	2674	call writediagfi(ngrid,"reevap","generic reevaporation of precipitation","kg m-2 s-1",2,reevap_precip_generic)
	2675	endif
	2676
[1477]	2677	if((hydrology).and.(.not.ok_slab_ocean))then
	2678	call writediagfi(ngrid,"hice","oceanic ice height","m",2,hice)
	2679	endif
[253]	2680
[1477]	2681	call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
[253]	2682
[1477]	2683	enddo ! end of 'nq' loop
[3232]	2684
[1477]	2685	endif ! end of 'tracer'
[253]	2686
[1477]	2687
	2688	! Output spectrum.
[3232]	2689	if(specOLR.and.corrk)then
[728]	2690	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
	2691	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
[2537]	2692	call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W/m^2/cm^-1",3,GSR_nu)
[526]	2693	endif
[253]	2694
[1836]	2695	#else
[2865]	2696	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
	2697	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
[2871]	2698	comm_ALBEQ(1:ngrid)=albedo_equivalent(1:ngrid)
	2699	if (.not.calldifv) comm_LATENT_HF(:)=0.0
[2865]	2700	if ((tracer).and.(water)) then
	2701	comm_CLOUDFRAC(1:ngrid,1:nlayer)=cloudfrac(1:ngrid,1:nlayer)
	2702	comm_TOTCLOUDFRAC(1:ngrid)=totcloudfrac(1:ngrid)
	2703	comm_SURFRAIN(1:ngrid)=zdqsrain(1:ngrid)
	2704	comm_DQVAP(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_h2o_vap)
[2871]	2705	comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_h2o_ice)
[2865]	2706	comm_H2OICE_REFF(1:ngrid,1:nlayer)=reffrad(1:ngrid,1:nlayer,iaero_h2o)
	2707	comm_REEVAP(1:ngrid)=reevap_precip(1:ngrid)
	2708	comm_DTRAIN(1:ngrid,1:nlayer)=zdtrain(1:ngrid,1:nlayer)
	2709	comm_DTLSC(1:ngrid,1:nlayer)=dtlscale(1:ngrid,1:nlayer)
	2710	comm_RH(1:ngrid,1:nlayer)=RH(1:ngrid,1:nlayer)
[2890]	2711
	2712	else if ((tracer).and.(generic_condensation).and.(.not. water)) then
	2713
	2714	! If you have set generic_condensation (and not water) and you have set several GCS
	2715	! then the outputs given to WRF will be only the ones for the last generic tracer
	2716	! (Because it is rewritten every tracer in the loop)
	2717	! WRF can take only one moist tracer
	2718
	2719	do iq=1,nq
	2720	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
[3232]	2721
[2890]	2722	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
	2723
	2724	reffrad_generic_zeros_for_wrf(:,:) = 1.
	2725
	2726	comm_CLOUDFRAC(1:ngrid,1:nlayer) = cloudfrac(1:ngrid,1:nlayer)
	2727	comm_TOTCLOUDFRAC(1:ngrid) = totcloudfrac(1:ngrid) !??????
	2728	comm_SURFRAIN(1:ngrid) = zdqsrain_generic(1:ngrid,iq)
	2729	comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_vap)
	2730	comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_generic_ice)
	2731	comm_H2OICE_REFF(1:ngrid,1:nlayer) = reffrad_generic_zeros_for_wrf(1:ngrid,1:nlayer) ! for now zeros !
	2732	!comm_H2OICE_REFF(1:ngrid,1:nlayer) = 0*zdtrain_generic(1:ngrid,1:nlayer) ! for now zeros !
	2733	comm_REEVAP(1:ngrid) = reevap_precip_generic(1:ngrid,iq)
	2734	comm_DTRAIN(1:ngrid,1:nlayer) = zdtrain_generic(1:ngrid,1:nlayer)
	2735	comm_DTLSC(1:ngrid,1:nlayer) = dt_generic_condensation(1:ngrid,1:nlayer)
	2736	comm_RH(1:ngrid,1:nlayer) = RH_generic(1:ngrid,1:nlayer,iq)
	2737
[3232]	2738	endif
	2739	end do ! do iq=1,nq loop on tracers
[2890]	2740
[2871]	2741	else
	2742	comm_CLOUDFRAC(1:ngrid,1:nlayer)=0.
	2743	comm_TOTCLOUDFRAC(1:ngrid)=0.
	2744	comm_SURFRAIN(1:ngrid)=0.
	2745	comm_DQVAP(1:ngrid,1:nlayer)=0.
	2746	comm_DQICE(1:ngrid,1:nlayer)=0.
	2747	comm_H2OICE_REFF(1:ngrid,1:nlayer)=0.
	2748	comm_REEVAP(1:ngrid)=0.
	2749	comm_DTRAIN(1:ngrid,1:nlayer)=0.
	2750	comm_DTLSC(1:ngrid,1:nlayer)=0.
	2751	comm_RH(1:ngrid,1:nlayer)=0.
[2890]	2752
	2753	endif ! if water, if generic_condensation, else
	2754
[2865]	2755	comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid)
	2756	comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid)
	2757	comm_FLUXTOP_LW(1:ngrid)=fluxtop_lw(1:ngrid)
	2758	comm_FLUXSURF_SW(1:ngrid)=fluxsurf_sw(1:ngrid)
	2759	comm_FLUXSURF_LW(1:ngrid)=fluxsurf_lw(1:ngrid)
	2760	comm_FLXGRD(1:ngrid)=fluxgrd(1:ngrid)
	2761	sensibFlux(1:ngrid) = zflubid(1:ngrid) - capcal(1:ngrid)*zdtsdif(1:ngrid) !!! ????
	2762	comm_HR_DYN(1:ngrid,1:nlayer) = zdtdyn(1:ngrid,1:nlayer)
[2019]	2763
[2867]	2764	! if (turb_resolved) then
	2765	! open(17,file='lsf.txt',form='formatted',status='old')
	2766	! rewind(17)
	2767	! DO l=1,nlayer
	2768	! read(17,*) lsf_dt(l),lsf_dq(l)
	2769	! ENDDO
	2770	! close(17)
	2771	! do ig=1,ngrid
	2772	! if ((tracer).and.(water)) then
	2773	! pdq(ig,:,igcm_h2o_vap) = pdq(ig,:,igcm_h2o_vap) + lsf_dq(:)
	2774	! endif
	2775	! pdt(ig,:) = pdt(ig,:) + lsf_dt(:)
	2776	! comm_HR_DYN(ig,:) = lsf_dt(:)
	2777	! enddo
	2778	! endif
[1836]	2779	#endif
	2780
[1622]	2781	! XIOS outputs
[3232]	2782	#ifdef CPP_XIOS
[1622]	2783	! Send fields to XIOS: (NB these fields must also be defined as
	2784	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
[1626]	2785	CALL send_xios_field("ls",zls)
[3232]	2786
[1622]	2787	CALL send_xios_field("ps",ps)
	2788	CALL send_xios_field("area",cell_area)
[2839]	2789	CALL send_xios_field("p",pplay)
[1622]	2790	CALL send_xios_field("temperature",zt)
	2791	CALL send_xios_field("u",zu)
	2792	CALL send_xios_field("v",zv)
[1877]	2793	CALL send_xios_field("omega",omega)
[3232]	2794
[2114]	2795	IF (calltherm) THEN
[2176]	2796	CALL send_xios_field('w_plm',zw2_bis)
[2232]	2797	CALL send_xios_field('entr',entr)
	2798	CALL send_xios_field('detr',detr)
	2799	! CALL send_xios_field('fm',fm_bis)
	2800	! CALL send_xios_field('fraca',fraca)
[2114]	2801	ENDIF
[3232]	2802
[2114]	2803	IF (water) THEN
	2804	CALL send_xios_field('h2o_vap',zq(:,:,igcm_h2o_vap))
	2805	CALL send_xios_field('h2o_ice',zq(:,:,igcm_h2o_ice))
[2985]	2806
	2807	CALL send_xios_field('h2o_layer1',zq(:,1,igcm_h2o_vap))
	2808	CALL send_xios_field('co2_layer1',zq(:,1,igcm_co2_ice))
	2809	CALL send_xios_field('tsurf',tsurf)
	2810	CALL send_xios_field('co2ice',qsurf(1:ngrid,igcm_co2_ice))
	2811	CALL send_xios_field('h2o_ice_s',qsurf(1:ngrid,igcm_h2o_ice))
[2114]	2812	ENDIF
[2891]	2813
[1877]	2814	CALL send_xios_field("ISR",fluxtop_dn)
	2815	CALL send_xios_field("OLR",fluxtop_lw)
[2733]	2816	CALL send_xios_field("ASR",fluxabs_sw)
[2734]	2817
[3232]	2818	if (specOLR .and. corrk) then
[2734]	2819	call send_xios_field("OLR3D",OLR_nu)
	2820	call send_xios_field("IR_Bandwidth",DWNI)
	2821	call send_xios_field("VI_Bandwidth",DWNV)
	2822	call send_xios_field("OSR3D",OSR_nu)
	2823	call send_xios_field("GSR3D",GSR_nu)
	2824	endif
	2825
[1682]	2826	if (lastcall.and.is_omp_master) then
	2827	write(,) "physiq: call xios_context_finalize"
	2828	call xios_context_finalize
	2829	endif
[1622]	2830	#endif
[3232]	2831
[2663]	2832	if (check_physics_outputs) then
	2833	! Check the validity of updated fields at the end of the physics step
	2834	call check_physics_fields("end of physiq:", zt, zu, zv, pplev, zq)
	2835	endif
	2836
[253]	2837	icount=icount+1
[3232]	2838
[2060]	2839	end subroutine physiq
[3232]	2840
[2060]	2841	end module physiq_mod

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

Download in other formats:

Original Format