Context Navigation

physiq.F90 @ 853

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

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

Download in other formats:

Original Format