Context Navigation

physiq_mod.F90 @ 1980

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

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

Download in other formats:

Original Format