Context Navigation

physiq_mod.F90 @ 2095

Last change on this file since 2095 was 2069, checked in by aboissinot, 7 years ago
f0 is now allocated only if calltherm=true. Useless thermal plume model flag "iflag_thermals_alim" is removed.
Property svn:executable set to ``*
File size: 99.2 KB

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

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

Download in other formats:

Original Format