Context Navigation

physiq_mod.F90 @ 3026

Last change on this file since 3026 was 2987, checked in by emillour, 18 months ago

Fixed the hardcoded value for the altitude of topmost level

In physiq_mod, the variable zzlev(:,nlayer+1) was hardcoded to 1e7. This is not enough when simulating the atmosphere of inflated hot Jupiter. Thus, to compute this value, we now use the thickness of the level below the top one to compute the altitude of topmost level. This variable is used only in the turbulent diffusion (we use 1/thickness) and the sedimentation (use of thickness).

Property svn:executable set to *

File size: 118.6 KB

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

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

Download in other formats:

Original Format