Context Navigation

physiq_mod.F90 @ 1782

Last change on this file since 1782 was 1699, checked in by jleconte, 8 years ago

21/04/2017 == JL

Add some arch files for a cluster in Bordeaux
added some is_master before printouts in callcorrk and physiq_mod
corrected a bug in bilinear big

Property svn:executable set to *

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

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

Download in other formats:

Original Format