Context Navigation

physiq.F90 @ 863

Last change on this file since 863 was 863, checked in by jleconte, 13 years ago

23/01/2013 == JL

Correction in largescale. a rneb factor was forgotten
Added some spectra in ave_stelpec
Corrected reevaporation in rain. Now conserve water better

Property svn:executable set to *

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

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

Download in other formats:

Original Format