Context Navigation

physiq.F90 @ 1448

Last change on this file since 1448 was 1429, checked in by msylvestre, 10 years ago
Some computations relative to the rings' shadow have been moved from physiq.f90 to call_rings.f90
Property svn:executable set to ``*
File size: 91.9 KB

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

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

Download in other formats:

Original Format