Context Navigation

physiq.F90 @ 988

Last change on this file since 988 was 969, checked in by emillour, 12 years ago
Generic/Universal? GCM: Further cleanup in outputs: enable output of a scalar in writediagfi in parallel mode and remove obsolete io* routines. EM
Property svn:executable set to ``*
File size: 76.8 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)
[867]	242	real zdmassmr_col(ngrid)
[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
[959]	328	real,save :: 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'
[918]	521	capcal(:)=1.e6
[952]	522	fluxgrd(:)=intheat
	523	print*,'Flux from ground = ',intheat,' W m^-2'
[253]	524	endif
	525	icount=1
	526
	527	! decide whether to update ice at end of run
	528	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	529	ice_update=.false.
	530	if(sourceevol)then
[955]	531	open(128,file='num_run',form='formatted', &
	532	status="old",iostat=ierr)
	533	if (ierr.ne.0) then
	534	write(,) "physiq: Error! No num_run file!"
	535	write(,) " (which is needed for sourceevol option)"
	536	stop
	537	endif
[253]	538	read(128,*) num_run
	539	close(128)
	540
[365]	541	if(num_run.ne.0.and.mod(num_run,2).eq.0)then
	542	!if(num_run.ne.0.and.mod(num_run,3).eq.0)then
[253]	543	print*,'Updating ice at end of this year!'
	544	ice_update=.true.
	545	ice_min(:)=1.e4
	546	endif
	547	endif
	548
	549	! define surface as continent or ocean
	550	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[728]	551	rnat(:)=1
[787]	552	do ig=1,ngrid
[253]	553	! if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
	554	if(inertiedat(ig,1).gt.1.E4)then
	555	rnat(ig)=0
	556	endif
	557	enddo
	558
	559	print*,'WARNING! Surface type currently decided by surface inertia'
	560	print*,'This should be improved e.g. in newstart.F'
	561
	562
	563	! initialise surface history variable
	564	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[728]	565	qsurf_hist(:,:)=qsurf(:,:)
[253]	566
	567	! initialise variable for dynamical heating diagnostic
	568	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	569	ztprevious(:,:)=pt(:,:)
	570
	571	! Set temperature just above condensation temperature (for Early Mars)
	572	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	573	if(nearco2cond) then
	574	write(,)' WARNING! Starting at Tcond+1K'
	575	do l=1, nlayer
	576	do ig=1,ngrid
	577	pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4 &
	578	-pt(ig,l)) / ptimestep
	579	enddo
	580	enddo
	581	endif
	582
	583	if(meanOLR)then
	584	! to record global radiative balance
	585	call system('rm -f rad_bal.out')
	586	! to record global mean/max/min temperatures
	587	call system('rm -f tem_bal.out')
	588	! to record global hydrological balance
	589	call system('rm -f h2o_bal.out')
	590	endif
	591
	592	watertest=.false.
	593	if(water)then
	594	! initialise variables for water cycle
	595
[365]	596	if(enertest)then
	597	watertest = .true.
	598	endif
	599
[728]	600	if(ice_update)then
	601	ice_initial(:)=qsurf(:,igcm_h2o_ice)
	602	endif
[253]	603
	604	endif
	605	call su_watercycle ! even if we don't have a water cycle, we might
	606	! need epsi for the wvp definitions in callcorrk.F
	607
	608	endif ! (end of "if firstcall")
	609
	610	! ---------------------------------------------------
	611	! 1.2 Initializations done at every physical timestep:
	612	! ---------------------------------------------------
	613
	614	! Initialize various variables
	615	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	616
[787]	617	pdu(1:ngrid,1:nlayermx) = 0.0
	618	pdv(1:ngrid,1:nlayermx) = 0.0
[253]	619	if ( .not.nearco2cond ) then
[787]	620	pdt(1:ngrid,1:nlayermx) = 0.0
[728]	621	endif
[787]	622	pdq(1:ngrid,1:nlayermx,1:nq) = 0.0
	623	pdpsrf(1:ngrid) = 0.0
	624	zflubid(1:ngrid) = 0.0
	625	zdtsurf(1:ngrid) = 0.0
	626	dqsurf(1:ngrid,1:nq)= 0.0
[253]	627
	628	zday=pday+ptime ! compute time, in sols (and fraction thereof)
	629
	630	! Compute Stellar Longitude (Ls)
	631	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	632	if (season) then
	633	call stellarlong(zday,zls)
	634	else
	635	call stellarlong(float(day_ini),zls)
	636	end if
	637
	638	! Compute geopotential between layers
	639	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	640
[787]	641	zzlay(1:ngrid,1:nlayermx)=pphi(1:ngrid,1:nlayermx)/g
[728]	642
[787]	643	zzlev(1:ngrid,1)=0.
	644	zzlev(1:ngrid,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
[728]	645
[253]	646	do l=2,nlayer
	647	do ig=1,ngrid
	648	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
	649	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
	650	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	651	enddo
	652	enddo
	653	! Potential temperature calculation may not be the same in physiq and dynamic...
	654
	655	! Compute potential temperature
	656	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	657
[597]	658	do l=1,nlayer
[787]	659	do ig=1,ngrid
[253]	660	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
[597]	661	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
[651]	662	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/g
	663	massarea(ig,l)=mass(ig,l)*area(ig)
[253]	664	enddo
	665	enddo
	666
	667	!-----------------------------------------------------------------------
	668	! 2. Compute radiative tendencies
	669	!-----------------------------------------------------------------------
	670
	671	if (callrad) then
[526]	672	if( mod(icount-1,iradia).eq.0.or.lastcall) then
[253]	673
	674	! Compute local stellar zenith angles
	675	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	676	call orbite(zls,dist_star,declin)
	677
	678	if (tlocked) then
	679	ztim1=SIN(declin)
[773]	680	! ztim2=COS(declin)COS(2.pi(zday/year_day) - zlsnres)
	681	! ztim3=-COS(declin)SIN(2.pi(zday/year_day) - zlsnres)
	682	! JL12 corrects tidally resonant cases. nres=omega_rot/omega_orb
[775]	683	ztim2=COS(declin)COS(2.pi(zday/year_day)nres - zls)
	684	ztim3=-COS(declin)SIN(2.pi(zday/year_day)nres - zls)
[253]	685
[728]	686	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
[253]	687	ztim1,ztim2,ztim3,mu0,fract)
	688
	689	elseif (diurnal) then
	690	ztim1=SIN(declin)
	691	ztim2=COS(declin)COS(2.pi*(zday-.5))
	692	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
	693
	694	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
	695	ztim1,ztim2,ztim3,mu0,fract)
	696
	697	else
	698
	699	call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
	700	! WARNING: this function appears not to work in 1D
	701
	702	endif
	703
	704	if (corrk) then
	705
	706	! a) Call correlated-k radiative transfer scheme
	707	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	708
	709	if(kastprof)then
[305]	710	print*,'kastprof should not = true here'
	711	call abort
[253]	712	endif
[538]	713	muvar(:,:)=0.0 ! only used for climate evolution studies (kcm1d) for now
	714
[526]	715	! standard callcorrk
	716	clearsky=.false.
[538]	717	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
[526]	718	albedo,emis,mu0,pplev,pplay,pt, &
[586]	719	tsurf,fract,dist_star,aerosol,muvar, &
[526]	720	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, &
[538]	721	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, &
[858]	722	tau_col,cloudfrac,totcloudfrac, &
[526]	723	clearsky,firstcall,lastcall)
[253]	724
[526]	725	! Option to call scheme once more for clear regions
[253]	726	if(CLFvarying)then
	727
[716]	728	! ---> PROBLEMS WITH ALLOCATED ARRAYS
	729	! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
[253]	730	clearsky=.true.
[538]	731	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
[253]	732	albedo,emis,mu0,pplev,pplay,pt, &
[586]	733	tsurf,fract,dist_star,aerosol,muvar, &
[253]	734	zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
[526]	735	fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1, &
[858]	736	tau_col1,cloudfrac,totcloudfrac, &
[538]	737	clearsky,firstcall,lastcall)
[716]	738	clearsky = .false. ! just in case.
[253]	739
	740	! Sum the fluxes and heating rates from cloudy/clear cases
	741	do ig=1,ngrid
	742	tf=totcloudfrac(ig)
	743	ntf=1.-tf
	744
[526]	745	fluxsurf_lw(ig) = ntffluxsurf_lw1(ig) + tffluxsurf_lw(ig)
	746	fluxsurf_sw(ig) = ntffluxsurf_sw1(ig) + tffluxsurf_sw(ig)
	747	fluxtop_lw(ig) = ntffluxtop_lw1(ig) + tffluxtop_lw(ig)
	748	fluxabs_sw(ig) = ntffluxabs_sw1(ig) + tffluxabs_sw(ig)
	749	tau_col(ig) = ntftau_col1(ig) + tftau_col(ig)
[253]	750
[728]	751	zdtlw(ig,1:nlayermx) = ntfzdtlw1(ig,1:nlayermx) + tfzdtlw(ig,1:nlayermx)
	752	zdtsw(ig,1:nlayermx) = ntfzdtsw1(ig,1:nlayermx) + tfzdtsw(ig,1:nlayermx)
[253]	753
[728]	754	OSR_nu(ig,1:L_NSPECTV) = ntfOSR_nu1(ig,1:L_NSPECTV) + tfOSR_nu(ig,1:L_NSPECTV)
[743]	755	OLR_nu(ig,1:L_NSPECTI) = ntfOLR_nu1(ig,1:L_NSPECTI) + tfOLR_nu(ig,1:L_NSPECTI)
[253]	756
[526]	757	enddo
[253]	758
[526]	759	endif !CLFvarying
[253]	760
	761	! Radiative flux from the sky absorbed by the surface (W.m-2)
	762	GSR=0.0
[787]	763	fluxrad_sky(1:ngrid)=emis(1:ngrid)fluxsurf_lw(1:ngrid)+fluxsurf_sw(1:ngrid)(1.-albedo(1:ngrid))
[253]	764
[952]	765	!if(noradsurf)then ! no lower surface; SW flux just disappears
	766	! GSR = SUM(fluxsurf_sw(1:ngrid)*area(1:ngrid))/totarea
	767	! fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
	768	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
	769	!endif
[253]	770
	771	! Net atmospheric radiative heating rate (K.s-1)
[787]	772	dtrad(1:ngrid,1:nlayermx)=zdtsw(1:ngrid,1:nlayermx)+zdtlw(1:ngrid,1:nlayermx)
[253]	773
	774	elseif(newtonian)then
	775
	776	! b) Call Newtonian cooling scheme
	777	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[787]	778	call newtrelax(ngrid,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
[253]	779
[787]	780	zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep
[253]	781	! e.g. surface becomes proxy for 1st atmospheric layer ?
	782
	783	else
	784
	785	! c) Atmosphere has no radiative effect
	786	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[787]	787	fluxtop_dn(1:ngrid) = fract(1:ngrid)mu0(1:ngrid)Fat1AU/dist_star**2
[728]	788	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
	789	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
	790	endif
[787]	791	fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
	792	fluxrad_sky(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid))
	793	fluxtop_lw(1:ngrid) = emis(1:ngrid)sigmatsurf(1:ngrid)**4
[728]	794	! radiation skips the atmosphere entirely
[253]	795
	796
[787]	797	dtrad(1:ngrid,1:nlayermx)=0.0
[728]	798	! hence no atmospheric radiative heating
	799
[253]	800	endif ! if corrk
	801
	802	endif ! of if(mod(icount-1,iradia).eq.0)
[787]	803
[253]	804
	805	! Transformation of the radiative tendencies
	806	! ------------------------------------------
	807
[787]	808	zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
	809	zplanck(1:ngrid)=emis(1:ngrid)sigmazplanck(1:ngrid)*zplanck(1:ngrid)
	810	fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
	811	pdt(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)+dtrad(1:ngrid,1:nlayermx)
[253]	812
	813	!-------------------------
	814	! test energy conservation
	815	if(enertest)then
[651]	816	dEtotSW = cppSUM(massarea(:,:)zdtsw(:,:))/totarea
	817	dEtotLW = cppSUM(massarea(:,:)zdtlw(:,:))/totarea
[959]	818	dEtotsSW = SUM(fluxsurf_sw(:)(1.-albedo(:))area(:))/totarea !JL13 carefull, albedo can have changed since the last time we called corrk
[726]	819	dEtotsLW = SUM((fluxsurf_lw(:)emis(:)-zplanck(:))area(:))/totarea
[651]	820	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
	821	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
[253]	822	print*,'---------------------------------------------------------------'
[594]	823	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
	824	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
	825	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
	826	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
	827	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
	828	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
[253]	829	endif
	830	!-------------------------
	831
	832	endif ! of if (callrad)
	833
	834	!-----------------------------------------------------------------------
	835	! 4. Vertical diffusion (turbulent mixing):
	836	! -----------------------------------------
	837
	838	if (calldifv) then
[526]	839
[787]	840	zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
[253]	841
[787]	842	zdum1(1:ngrid,1:nlayermx)=0.0
	843	zdum2(1:ngrid,1:nlayermx)=0.0
[253]	844
[594]	845
	846	!JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff
	847	if (UseTurbDiff) then
	848
	849	call turbdiff(ngrid,nlayer,nq,rnat, &
[253]	850	ptimestep,capcal,lwrite, &
	851	pplay,pplev,zzlay,zzlev,z0, &
[594]	852	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
	853	zdum1,zdum2,pdt,pdq,zflubid, &
	854	zdudif,zdvdif,zdtdif,zdtsdif, &
[728]	855	sensibFlux,q2,zdqdif,zdqevap,zdqsdif,lastcall)
[594]	856
	857	else
	858
[787]	859	zdh(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)/zpopsk(1:ngrid,1:nlayermx)
[594]	860
	861	call vdifc(ngrid,nlayer,nq,rnat,zpopsk, &
	862	ptimestep,capcal,lwrite, &
	863	pplay,pplev,zzlay,zzlev,z0, &
[253]	864	pu,pv,zh,pq,tsurf,emis,qsurf, &
	865	zdum1,zdum2,zdh,pdq,zflubid, &
[594]	866	zdudif,zdvdif,zdhdif,zdtsdif, &
	867	sensibFlux,q2,zdqdif,zdqsdif,lastcall)
[253]	868
[787]	869	zdtdif(1:ngrid,1:nlayermx)=zdhdif(1:ngrid,1:nlayermx)*zpopsk(1:ngrid,1:nlayermx) ! for diagnostic only
	870	zdqevap(1:ngrid,1:nlayermx)=0.
[594]	871
	872	end if !turbdiff
	873
[787]	874	pdv(1:ngrid,1:nlayermx)=pdv(1:ngrid,1:nlayermx)+zdvdif(1:ngrid,1:nlayermx)
	875	pdu(1:ngrid,1:nlayermx)=pdu(1:ngrid,1:nlayermx)+zdudif(1:ngrid,1:nlayermx)
	876	pdt(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)+zdtdif(1:ngrid,1:nlayermx)
	877	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
[253]	878	if (tracer) then
[787]	879	pdq(1:ngrid,1:nlayermx,1:nq)=pdq(1:ngrid,1:nlayermx,1:nq)+ zdqdif(1:ngrid,1:nlayermx,1:nq)
	880	dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
[253]	881	end if ! of if (tracer)
	882
	883	!-------------------------
	884	! test energy conservation
	885	if(enertest)then
[651]	886	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
[253]	887	do ig = 1, ngrid
[651]	888	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
[622]	889	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
[253]	890	enddo
[651]	891	dEtot = SUM(dEdiff(:)*area(:))/totarea
	892	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
	893	dEtots = SUM(dEdiffs(:)*area(:))/totarea
	894	AtmToSurf_TurbFlux=SUM(sensibFlux(:)*area(:))/totarea
[597]	895	if (UseTurbDiff) then
	896	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
	897	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
	898	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
	899	else
	900	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
	901	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
	902	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
	903	end if
[594]	904	! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
	905	! but not given back elsewhere
[253]	906	endif
	907	!-------------------------
	908
	909	!-------------------------
	910	! test water conservation
	911	if(watertest.and.water)then
[651]	912	dWtot = SUM(massarea(:,:)zdqdif(:,:,igcm_h2o_vap))ptimestep/totarea
	913	dWtots = SUM(zdqsdif(:,igcm_h2o_vap)area(:))ptimestep/totarea
[253]	914	do ig = 1, ngrid
[651]	915	vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
	916	Enddo
	917	dWtot = dWtot + SUM(massarea(:,:)zdqdif(:,:,igcm_h2o_ice))ptimestep/totarea
	918	dWtots = dWtots + SUM(zdqsdif(:,igcm_h2o_ice)area(:))ptimestep/totarea
	919	do ig = 1, ngrid
	920	vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
	921	Enddo
	922	nconsMAX=MAXVAL(vdifcncons(:))
[253]	923
	924	print*,'---------------------------------------------------------------'
	925	print*,'In difv atmospheric water change =',dWtot,' kg m-2'
	926	print*,'In difv surface water change =',dWtots,' kg m-2'
	927	print*,'In difv non-cons factor =',dWtot+dWtots,' kg m-2'
	928	print*,'In difv MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
	929
	930	endif
	931	!-------------------------
	932
	933	else
	934
	935	if(.not.newtonian)then
	936
[787]	937	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
[253]	938
	939	endif
	940
	941	endif ! of if (calldifv)
	942
	943
	944	!-----------------------------------------------------------------------
	945	! 5. Dry convective adjustment:
	946	! -----------------------------
	947
	948	if(calladj) then
	949
[787]	950	zdh(1:ngrid,1:nlayermx) = pdt(1:ngrid,1:nlayermx)/zpopsk(1:ngrid,1:nlayermx)
	951	zduadj(1:ngrid,1:nlayermx)=0.0
	952	zdvadj(1:ngrid,1:nlayermx)=0.0
	953	zdhadj(1:ngrid,1:nlayermx)=0.0
[253]	954
	955
[586]	956	call convadj(ngrid,nlayer,nq,ptimestep, &
	957	pplay,pplev,zpopsk, &
	958	pu,pv,zh,pq, &
	959	pdu,pdv,zdh,pdq, &
	960	zduadj,zdvadj,zdhadj, &
	961	zdqadj)
[253]	962
[787]	963	pdu(1:ngrid,1:nlayermx) = pdu(1:ngrid,1:nlayermx) + zduadj(1:ngrid,1:nlayermx)
	964	pdv(1:ngrid,1:nlayermx) = pdv(1:ngrid,1:nlayermx) + zdvadj(1:ngrid,1:nlayermx)
	965	pdt(1:ngrid,1:nlayermx) = pdt(1:ngrid,1:nlayermx) + zdhadj(1:ngrid,1:nlayermx)*zpopsk(1:ngrid,1:nlayermx)
	966	zdtadj(1:ngrid,1:nlayermx) = zdhadj(1:ngrid,1:nlayermx)*zpopsk(1:ngrid,1:nlayermx) ! for diagnostic only
[728]	967
[253]	968	if(tracer) then
[787]	969	pdq(1:ngrid,1:nlayermx,1:nq) = pdq(1:ngrid,1:nlayermx,1:nq) + zdqadj(1:ngrid,1:nlayermx,1:nq)
[253]	970	end if
	971
	972	!-------------------------
	973	! test energy conservation
	974	if(enertest)then
[651]	975	dEtot=cppSUM(massarea(:,:)zdtadj(:,:))/totarea
[594]	976	print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
[253]	977	endif
	978	!-------------------------
	979
	980	!-------------------------
	981	! test water conservation
	982	if(watertest)then
[651]	983	dWtot = SUM(massarea(:,:)zdqadj(:,:,igcm_h2o_vap))ptimestep/totarea
[253]	984	do ig = 1, ngrid
[651]	985	cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
	986	Enddo
	987	dWtot = dWtot + SUM(massarea(:,:)zdqadj(:,:,igcm_h2o_ice))ptimestep/totarea
	988	do ig = 1, ngrid
	989	cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
	990	Enddo
	991	nconsMAX=MAXVAL(cadjncons(:))
[253]	992
	993	print*,'In convadj atmospheric water change =',dWtot,' kg m-2'
[651]	994	print*,'In convadj MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
[253]	995	endif
	996	!-------------------------
[787]	997
[253]	998	endif ! of if(calladj)
	999
	1000	!-----------------------------------------------------------------------
	1001	! 6. Carbon dioxide condensation-sublimation:
	1002	! -------------------------------------------
	1003
	1004	if (co2cond) then
	1005	if (.not.tracer) then
	1006	print*,'We need a CO2 ice tracer to condense CO2'
	1007	call abort
	1008	endif
[305]	1009	call condense_cloud(ngrid,nlayer,nq,ptimestep, &
[253]	1010	capcal,pplay,pplev,tsurf,pt, &
	1011	pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, &
	1012	qsurf(1,igcm_co2_ice),albedo,emis, &
	1013	zdtc,zdtsurfc,pdpsrf,zduc,zdvc, &
[858]	1014	zdqc)
[253]	1015
	1016
[787]	1017	pdt(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)+zdtc(1:ngrid,1:nlayermx)
	1018	pdv(1:ngrid,1:nlayermx)=pdv(1:ngrid,1:nlayermx)+zdvc(1:ngrid,1:nlayermx)
	1019	pdu(1:ngrid,1:nlayermx)=pdu(1:ngrid,1:nlayermx)+zduc(1:ngrid,1:nlayermx)
	1020	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
[728]	1021
[787]	1022	pdq(1:ngrid,1:nlayermx,1:nq)=pdq(1:ngrid,1:nlayermx,1:nq)+ zdqc(1:ngrid,1:nlayermx,1:nq)
[253]	1023	! Note: we do not add surface co2ice tendency
	1024	! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
	1025
	1026	!-------------------------
	1027	! test energy conservation
	1028	if(enertest)then
[651]	1029	dEtot = cppSUM(massarea(:,:)zdtc(:,:))/totarea
	1030	dEtots = SUM(capcal(:)zdtsurfc(:)area(:))/totarea
[253]	1031	print*,'In co2cloud atmospheric energy change =',dEtot,' W m-2'
	1032	print*,'In co2cloud surface energy change =',dEtots,' W m-2'
	1033	endif
	1034	!-------------------------
	1035
	1036	endif ! of if (co2cond)
	1037
	1038
	1039	!-----------------------------------------------------------------------
	1040	! 7. Specific parameterizations for tracers
	1041	! -----------------------------------------
	1042
	1043	if (tracer) then
	1044
	1045	! 7a. Water and ice
	1046	! ---------------
	1047	if (water) then
	1048
	1049	! ----------------------------------------
	1050	! Water ice condensation in the atmosphere
	1051	! ----------------------------------------
[728]	1052	if(watercond.and.(RLVTT.gt.1.e-8))then
[253]	1053
[728]	1054	! ----------------
	1055	! Moist convection
	1056	! ----------------
[773]	1057
[787]	1058	dqmoist(1:ngrid,1:nlayermx,1:nq)=0.
	1059	dtmoist(1:ngrid,1:nlayermx)=0.
[253]	1060
[787]	1061	call moistadj(ngrid,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
[253]	1062
[787]	1063	pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) = pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) &
	1064	+dqmoist(1:ngrid,1:nlayermx,igcm_h2o_vap)
	1065	pdq(1:ngrid,1:nlayermx,igcm_h2o_ice) =pdq(1:ngrid,1:nlayermx,igcm_h2o_ice) &
	1066	+dqmoist(1:ngrid,1:nlayermx,igcm_h2o_ice)
	1067	pdt(1:ngrid,1:nlayermx) = pdt(1:ngrid,1:nlayermx)+dtmoist(1:ngrid,1:nlayermx)
[728]	1068
[253]	1069	!-------------------------
	1070	! test energy conservation
	1071	if(enertest)then
[728]	1072	dEtot=cppSUM(massarea(:,:)dtmoist(:,:))/totarea
[787]	1073	do ig=1,ngrid
[728]	1074	madjdE(ig) = cppSUM(mass(:,:)dtmoist(:,:))
[253]	1075	enddo
	1076	print*,'In moistadj atmospheric energy change =',dEtot,' W m-2'
[773]	1077	print,'In moistadj MAX atmospheric energy change =',MAXVAL(dtmoist(:,:))ptimestep,'K/step'
	1078	print,'In moistadj MIN atmospheric energy change =',MINVAL(dtmoist(:,:))ptimestep,'K/step'
[651]	1079
	1080	! test energy conservation
	1081	dWtot = SUM(massarea(:,:)dqmoist(:,:,igcm_h2o_vap))ptimestep/totarea
	1082	dWtot = dWtot + SUM(massarea(:,:)dqmoist(:,:,igcm_h2o_ice))ptimestep/totarea
[253]	1083	print*,'In moistadj atmospheric water change =',dWtot,' kg m-2'
	1084	endif
	1085	!-------------------------
	1086
	1087
[728]	1088	! --------------------------------
	1089	! Large scale condensation/evaporation
	1090	! --------------------------------
[253]	1091
[787]	1092	call largescale(ngrid,nq,ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
[253]	1093
[787]	1094	pdt(1:ngrid,1:nlayermx) = pdt(1:ngrid,1:nlayermx)+dtlscale(1:ngrid,1:nlayermx)
	1095	pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) = pdq(1:ngrid,1:nlayermx,igcm_h2o_vap)+dqvaplscale(1:ngrid,1:nlayermx)
	1096	pdq(1:ngrid,1:nlayermx,igcm_h2o_ice) = pdq(1:ngrid,1:nlayermx,igcm_h2o_ice)+dqcldlscale(1:ngrid,1:nlayermx)
[253]	1097
	1098	!-------------------------
	1099	! test energy conservation
	1100	if(enertest)then
[787]	1101	do ig=1,ngrid
[728]	1102	lscaledE(ig) = cppSUM(mass(:,:)dtlscale(:,:))
[253]	1103	enddo
[728]	1104	dEtot=cppSUM(massarea(:,:)(dtlscale(:,:)))/totarea
	1105	if(isnan(dEtot)) then
	1106	print*,'Nan in largescale, abort'
	1107	STOP
	1108	endif
	1109	print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
	1110
	1111	! test water conservation
	1112	dWtot = SUM(massarea(:,:)dqvaplscale(:,:))ptimestep/totarea
	1113	dWtot = dWtot + SUM(massarea(:,:)dqcldlscale(:,:))ptimestep/totarea
	1114	print*,'In largescale atmospheric water change =',dWtot,' kg m-2'
[253]	1115	endif
	1116	!-------------------------
	1117
	1118	! compute cloud fraction
	1119	do l = 1, nlayer
[787]	1120	do ig=1,ngrid
[253]	1121	cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
	1122	enddo
	1123	enddo
	1124
[728]	1125	endif ! of if (watercondense)
[253]	1126
	1127
	1128	! --------------------------------
	1129	! Water ice / liquid precipitation
	1130	! --------------------------------
[728]	1131	if(waterrain)then
[253]	1132
[787]	1133	zdqrain(1:ngrid,1:nlayermx,1:nq) = 0.0
	1134	zdqsrain(1:ngrid) = 0.0
	1135	zdqssnow(1:ngrid) = 0.0
[253]	1136
[787]	1137	call rain(ngrid,nq,ptimestep,pplev,pplay,pt,pdt,pq,pdq, &
[253]	1138	zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
	1139
[787]	1140	pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) = pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) &
	1141	+zdqrain(1:ngrid,1:nlayermx,igcm_h2o_vap)
	1142	pdq(1:ngrid,1:nlayermx,igcm_h2o_ice) = pdq(1:ngrid,1:nlayermx,igcm_h2o_ice) &
	1143	+zdqrain(1:ngrid,1:nlayermx,igcm_h2o_ice)
	1144	pdt(1:ngrid,1:nlayermx) = pdt(1:ngrid,1:nlayermx)+zdtrain(1:ngrid,1:nlayermx)
	1145	dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid) ! a bug was here
	1146	dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
	1147	rainout(1:ngrid) = zdqsrain(1:ngrid)+zdqssnow(1:ngrid) ! diagnostic
[253]	1148
	1149
[651]	1150	!-------------------------
	1151	! test energy conservation
	1152	if(enertest)then
	1153	dEtot=cppSUM(massarea(:,:)zdtrain(:,:))/totarea
[253]	1154	print*,'In rain atmospheric T energy change =',dEtot,' W m-2'
[651]	1155	dItot = SUM(massarea(:,:)zdqrain(:,:,igcm_h2o_ice))/totareaRLVTT/cpp
	1156	dItot = dItot + SUM(area(:)zdqssnow(:))/totareaRLVTT/cpp
	1157	dVtot = SUM(massarea(:,:)zdqrain(:,:,igcm_h2o_vap))ptimestep/totarea
[728]	1158	dVtot = dVtot + SUM(area(:)zdqsrain(:))/totareaRLVTT/cpp
[651]	1159	dEtot = dItot + dVtot
	1160	print*,'In rain dItot =',dItot,' W m-2'
	1161	print*,'In rain dVtot =',dVtot,' W m-2'
[253]	1162	print*,'In rain atmospheric L energy change =',dEtot,' W m-2'
	1163
[651]	1164	! test water conservation
	1165	dWtot = SUM(massarea(:,:)zdqrain(:,:,igcm_h2o_vap))ptimestep/totarea
	1166	dWtot = dWtot + SUM(massarea(:,:)zdqrain(:,:,igcm_h2o_ice))ptimestep/totarea
	1167	dWtots = SUM((zdqsrain(:)+zdqssnow(:))area(:))ptimestep/totarea
[253]	1168	print*,'In rain atmospheric water change =',dWtot,' kg m-2'
	1169	print*,'In rain surface water change =',dWtots,' kg m-2'
	1170	print*,'In rain non-cons factor =',dWtot+dWtots,' kg m-2'
	1171	endif
	1172	!-------------------------
	1173
[728]	1174	end if ! of if (waterrain)
	1175	end if ! of if (water)
[253]	1176
	1177
	1178	! 7c. Aerosol particles
	1179	! -------------------
	1180	! -------------
	1181	! Sedimentation
	1182	! -------------
	1183	if (sedimentation) then
[787]	1184	zdqsed(1:ngrid,1:nlayermx,1:nq) = 0.0
	1185	zdqssed(1:ngrid,1:nq) = 0.0
[253]	1186
	1187
	1188	!-------------------------
	1189	! find qtot
	1190	if(watertest)then
[959]	1191	iq=igcm_h2o_ice
[651]	1192	dWtot = SUM(massarea(:,:)pq(:,:,iq))ptimestep/totarea
	1193	dWtots = SUM(massarea(:,:)pdq(:,:,iq))ptimestep/totarea
[253]	1194	print*,'Before sedim pq =',dWtot,' kg m-2'
	1195	print*,'Before sedim pdq =',dWtots,' kg m-2'
	1196	endif
	1197	!-------------------------
	1198
	1199	call callsedim(ngrid,nlayer,ptimestep, &
[858]	1200	pplev,zzlev,pt,pdt,pq,pdq,zdqsed,zdqssed,nq)
[253]	1201
	1202	!-------------------------
	1203	! find qtot
	1204	if(watertest)then
[959]	1205	iq=igcm_h2o_ice
[651]	1206	dWtot = SUM(massarea(:,:)pq(:,:,iq))ptimestep/totarea
	1207	dWtots = SUM(massarea(:,:)pdq(:,:,iq))ptimestep/totarea
[253]	1208	print*,'After sedim pq =',dWtot,' kg m-2'
	1209	print*,'After sedim pdq =',dWtots,' kg m-2'
	1210	endif
	1211	!-------------------------
	1212
[728]	1213	! for now, we only allow H2O ice to sediment
[253]	1214	! and as in rain.F90, whether it falls as rain or snow depends
	1215	! only on the surface temperature
[787]	1216	pdq(1:ngrid,1:nlayermx,1:nq) = pdq(1:ngrid,1:nlayermx,1:nq) + zdqsed(1:ngrid,1:nlayermx,1:nq)
	1217	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
[253]	1218
	1219	!-------------------------
	1220	! test water conservation
	1221	if(watertest)then
[651]	1222	dWtot = SUM(massarea(:,:)(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice)))ptimestep/totarea
	1223	dWtots = SUM((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))area(:))ptimestep/totarea
[253]	1224	print*,'In sedim atmospheric ice change =',dWtot,' kg m-2'
	1225	print*,'In sedim surface ice change =',dWtots,' kg m-2'
	1226	print*,'In sedim non-cons factor =',dWtot+dWtots,' kg m-2'
	1227	endif
	1228	!-------------------------
	1229
	1230	end if ! of if (sedimentation)
	1231
	1232
	1233	! 7d. Updates
	1234	! ---------
	1235
[728]	1236	! -----------------------------------
	1237	! Updating atm mass and tracer budget
	1238	! -----------------------------------
	1239
	1240	if(mass_redistrib) then
	1241
[787]	1242	zdmassmr(1:ngrid,1:nlayermx) = mass(1:ngrid,1:nlayermx) * &
	1243	( zdqevap(1:ngrid,1:nlayermx) &
	1244	+ zdqrain(1:ngrid,1:nlayermx,igcm_h2o_vap) &
	1245	+ dqmoist(1:ngrid,1:nlayermx,igcm_h2o_vap) &
	1246	+ dqvaplscale(1:ngrid,1:nlayermx) )
[863]	1247
	1248	do ig = 1, ngrid
	1249	zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayermx))
	1250	enddo
[728]	1251
[787]	1252	call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
[863]	1253	call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
[787]	1254	call writediagfi(ngrid,"mass","mass"," ",3,mass)
[728]	1255
	1256	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
	1257	rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
	1258	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
	1259	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
	1260
	1261
[787]	1262	pdq(1:ngrid,1:nlayermx,1:nq) = pdq(1:ngrid,1:nlayermx,1:nq) + zdqmr(1:ngrid,1:nlayermx,1:nq)
	1263	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
	1264	pdt(1:ngrid,1:nlayermx) = pdt(1:ngrid,1:nlayermx) + zdtmr(1:ngrid,1:nlayermx)
	1265	pdu(1:ngrid,1:nlayermx) = pdu(1:ngrid,1:nlayermx) + zdumr(1:ngrid,1:nlayermx)
	1266	pdv(1:ngrid,1:nlayermx) = pdv(1:ngrid,1:nlayermx) + zdvmr(1:ngrid,1:nlayermx)
	1267	pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
	1268	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
[728]	1269
	1270	! print,'after mass redistrib, q=',pq(211,1:nlayermx,igcm_h2o_vap)+ptimesteppdq(211,1:nlayermx,igcm_h2o_vap)
	1271	endif
	1272
	1273
	1274
[253]	1275	! ---------------------------------
	1276	! Updating tracer budget on surface
	1277	! ---------------------------------
	1278
	1279	if(hydrology)then
	1280
[787]	1281	call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, &
[253]	1282	capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice)
	1283	! note: for now, also changes albedo in the subroutine
	1284
[787]	1285	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurf_hyd(1:ngrid)
	1286	qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqs_hyd(1:ngrid,1:nq)
[253]	1287	! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
	1288
	1289	!-------------------------
	1290	! test energy conservation
	1291	if(enertest)then
[651]	1292	dEtots = SUM(area(:)capcal(:)zdtsurf_hyd(:))/totarea
	1293	print*,'In hydrol surface energy change =',dEtots,' W m-2'
[253]	1294	endif
	1295	!-------------------------
	1296
	1297	!-------------------------
	1298	! test water conservation
	1299	if(watertest)then
[651]	1300	dWtots = SUM(dqs_hyd(:,igcm_h2o_ice)area(:))ptimestep/totarea
[253]	1301	print*,'In hydrol surface ice change =',dWtots,' kg m-2'
[651]	1302	dWtots = SUM(dqs_hyd(:,igcm_h2o_vap)area(:))ptimestep/totarea
[253]	1303	print*,'In hydrol surface water change =',dWtots,' kg m-2'
	1304	print*,'---------------------------------------------------------------'
	1305	endif
	1306	!-------------------------
	1307
	1308	ELSE ! of if (hydrology)
	1309
[787]	1310	qsurf(1:ngrid,1:nq)=qsurf(1:ngrid,1:nq)+ptimestep*dqsurf(1:ngrid,1:nq)
[253]	1311
	1312	END IF ! of if (hydrology)
	1313
	1314	! Add qsurf to qsurf_hist, which is what we save in
	1315	! diagfi.nc etc. At the same time, we set the water
	1316	! content of ocean gridpoints back to zero, in order
	1317	! to avoid rounding errors in vdifc, rain
[622]	1318	qsurf_hist(:,:) = qsurf(:,:)
[253]	1319
	1320	if(ice_update)then
[787]	1321	ice_min(1:ngrid)=min(ice_min(1:ngrid),qsurf(1:ngrid,igcm_h2o_ice))
[253]	1322	endif
	1323
	1324	endif ! of if (tracer)
	1325
	1326	!-----------------------------------------------------------------------
	1327	! 9. Surface and sub-surface soil temperature
	1328	!-----------------------------------------------------------------------
	1329
	1330	! Increment surface temperature
[787]	1331	tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
[253]	1332
	1333	! Compute soil temperatures and subsurface heat flux
	1334	if (callsoil) then
[787]	1335	call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, &
[253]	1336	ptimestep,tsurf,tsoil,capcal,fluxgrd)
	1337	endif
	1338
	1339	!-------------------------
	1340	! test energy conservation
	1341	if(enertest)then
[728]	1342	dEtots = SUM(area(:)capcal(:)zdtsurf(:))/totarea
[597]	1343	print*,'Surface energy change =',dEtots,' W m-2'
[253]	1344	endif
	1345	!-------------------------
	1346
	1347	!-----------------------------------------------------------------------
	1348	! 10. Perform diagnostics and write output files
	1349	!-----------------------------------------------------------------------
	1350
	1351	! -------------------------------
	1352	! Dynamical fields incrementation
	1353	! -------------------------------
	1354	! For output only: the actual model integration is performed in the dynamics
	1355
	1356	! temperature, zonal and meridional wind
[787]	1357	zt(1:ngrid,1:nlayermx) = pt(1:ngrid,1:nlayermx) + pdt(1:ngrid,1:nlayermx)*ptimestep
	1358	zu(1:ngrid,1:nlayermx) = pu(1:ngrid,1:nlayermx) + pdu(1:ngrid,1:nlayermx)*ptimestep
	1359	zv(1:ngrid,1:nlayermx) = pv(1:ngrid,1:nlayermx) + pdv(1:ngrid,1:nlayermx)*ptimestep
[253]	1360
[728]	1361	! diagnostic
[787]	1362	zdtdyn(1:ngrid,1:nlayermx) = pt(1:ngrid,1:nlayermx)-ztprevious(1:ngrid,1:nlayermx)
	1363	ztprevious(1:ngrid,1:nlayermx) = zt(1:ngrid,1:nlayermx)
[253]	1364
	1365	if(firstcall)then
[787]	1366	zdtdyn(1:ngrid,1:nlayermx)=0.0
[253]	1367	endif
	1368
	1369	! dynamical heating diagnostic
	1370	do ig=1,ngrid
[728]	1371	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp/ptimestep
[253]	1372	enddo
	1373
	1374	! tracers
[787]	1375	zq(1:ngrid,1:nlayermx,1:nq) = pq(1:ngrid,1:nlayermx,1:nq) + pdq(1:ngrid,1:nlayermx,1:nq)*ptimestep
[253]	1376
	1377	! surface pressure
[787]	1378	ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
[253]	1379
	1380	! pressure
	1381	do l=1,nlayer
[787]	1382	zplev(1:ngrid,l) = pplev(1:ngrid,l)/pplev(1:ngrid,1)*ps(:)
	1383	zplay(1:ngrid,l) = pplay(1:ngrid,l)/pplev(1:ngrid,1)*ps(1:ngrid)
[253]	1384	enddo
	1385
	1386	! ---------------------------------------------------------
	1387	! Surface and soil temperature information
	1388	! ---------------------------------------------------------
	1389
[651]	1390	Ts1 = SUM(area(:)*tsurf(:))/totarea
	1391	Ts2 = MINVAL(tsurf(:))
	1392	Ts3 = MAXVAL(tsurf(:))
[253]	1393	if(callsoil)then
[651]	1394	TsS = SUM(area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
[880]	1395	if (cursor == 1) then
	1396	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
	1397	print*,Ts1,Ts2,Ts3,TsS
	1398	endif
[959]	1399	else
	1400	if (cursor == 1) then
	1401	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
	1402	print*,Ts1,Ts2,Ts3
	1403	endif
	1404	end if
[253]	1405
	1406	! ---------------------------------------------------------
	1407	! Check the energy balance of the simulation during the run
	1408	! ---------------------------------------------------------
	1409
	1410	if(corrk)then
	1411
[651]	1412	ISR = SUM(area(:)*fluxtop_dn(:))/totarea
	1413	ASR = SUM(area(:)*fluxabs_sw(:))/totarea
	1414	OLR = SUM(area(:)*fluxtop_lw(:))/totarea
	1415	GND = SUM(area(:)*fluxgrd(:))/totarea
	1416	DYN = SUM(area(:)*fluxdyn(:))/totarea
[787]	1417	do ig=1,ngrid
[253]	1418	if(fluxtop_dn(ig).lt.0.0)then
	1419	print*,'fluxtop_dn has gone crazy'
	1420	print*,'fluxtop_dn=',fluxtop_dn(ig)
	1421	print*,'tau_col=',tau_col(ig)
	1422	print*,'aerosol=',aerosol(ig,:,:)
	1423	print*,'temp= ',pt(ig,:)
	1424	print*,'pplay= ',pplay(ig,:)
	1425	call abort
	1426	endif
	1427	end do
	1428
[787]	1429	if(ngrid.eq.1)then
[253]	1430	DYN=0.0
	1431	endif
	1432
[959]	1433	print*,' ISR ASR OLR GND DYN [W m^-2]'
	1434	print*, ISR,ASR,OLR,GND,DYN
[961]	1435
[253]	1436	if(enertest)then
[651]	1437	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
	1438	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
	1439	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
[253]	1440	endif
	1441
	1442	if(meanOLR)then
[787]	1443	if((ngrid.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
[253]	1444	! to record global radiative balance
[588]	1445	open(92,file="rad_bal.out",form='formatted',position='append')
[651]	1446	write(92,*) zday,ISR,ASR,OLR
[253]	1447	close(92)
[588]	1448	open(93,file="tem_bal.out",form='formatted',position='append')
[253]	1449	write(93,*) zday,Ts1,Ts2,Ts3,TsS
	1450	close(93)
	1451	endif
	1452	endif
	1453
	1454	endif
	1455
[651]	1456
[253]	1457	! ------------------------------------------------------------------
	1458	! Diagnostic to test radiative-convective timescales in code
	1459	! ------------------------------------------------------------------
	1460	if(testradtimes)then
[588]	1461	open(38,file="tau_phys.out",form='formatted',position='append')
[253]	1462	ig=1
	1463	do l=1,nlayer
	1464	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
	1465	enddo
	1466	close(38)
[726]	1467	print*,'As testradtimes enabled,'
	1468	print*,'exiting physics on first call'
[253]	1469	call abort
	1470	endif
	1471
	1472	! ---------------------------------------------------------
	1473	! Compute column amounts (kg m-2) if tracers are enabled
	1474	! ---------------------------------------------------------
	1475	if(tracer)then
[787]	1476	qcol(1:ngrid,1:nq)=0.0
[253]	1477	do iq=1,nq
[787]	1478	do ig=1,ngrid
[728]	1479	qcol(ig,iq) = SUM( zq(ig,1:nlayermx,iq) * mass(ig,1:nlayermx))
	1480	enddo
[253]	1481	enddo
	1482
[726]	1483	! Generalised for arbitrary aerosols now. (LK)
[787]	1484	reffcol(1:ngrid,1:naerkind)=0.0
[728]	1485	if(co2cond.and.(iaero_co2.ne.0))then
[858]	1486	call co2_reffrad(ngrid,nq,zq,reffrad(1,1,iaero_co2))
[787]	1487	do ig=1,ngrid
[728]	1488	reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayermx,igcm_co2_ice)reffrad(ig,1:nlayermx,iaero_co2)mass(ig,1:nlayermx))
[253]	1489	enddo
[728]	1490	endif
	1491	if(water.and.(iaero_h2o.ne.0))then
[858]	1492	call h2o_reffrad(ngrid,zq(1,1,igcm_h2o_ice),zt, &
	1493	reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
[787]	1494	do ig=1,ngrid
[728]	1495	reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayermx,igcm_h2o_ice)reffrad(ig,1:nlayermx,iaero_h2o)mass(ig,1:nlayermx))
	1496	enddo
	1497	endif
[253]	1498
	1499	endif
	1500
	1501	! ---------------------------------------------------------
	1502	! Test for water conservation if water is enabled
	1503	! ---------------------------------------------------------
	1504
	1505	if(water)then
	1506
[651]	1507	icesrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_ice))/totarea
	1508	liqsrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_vap))/totarea
	1509	icecol = SUM(area(:)*qcol(:,igcm_h2o_ice))/totarea
	1510	vapcol = SUM(area(:)*qcol(:,igcm_h2o_vap))/totarea
[253]	1511
[651]	1512	h2otot = icesrf + liqsrf + icecol + vapcol
[253]	1513
[651]	1514	print*,' Total water amount [kg m^-2]: ',h2otot
[253]	1515	print*,' Surface ice Surface liq. Atmos. con. Atmos. vap. [kg m^-2] '
[651]	1516	print*, icesrf,liqsrf,icecol,vapcol
[253]	1517
	1518	if(meanOLR)then
[787]	1519	if((ngrid.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
[253]	1520	! to record global water balance
[588]	1521	open(98,file="h2o_bal.out",form='formatted',position='append')
[651]	1522	write(98,*) zday,icesrf,liqsrf,icecol,vapcol
[253]	1523	close(98)
	1524	endif
	1525	endif
	1526
	1527	endif
	1528
	1529	! ---------------------------------------------------------
	1530	! Calculate RH for diagnostic if water is enabled
	1531	! ---------------------------------------------------------
	1532
	1533	if(water)then
	1534	do l = 1, nlayer
[787]	1535	do ig=1,ngrid
[728]	1536	! call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
	1537	call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l))
	1538
[253]	1539	RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
	1540	enddo
	1541	enddo
	1542
	1543	! compute maximum possible H2O column amount (100% saturation)
	1544	do ig=1,ngrid
[728]	1545	H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
[253]	1546	enddo
	1547
	1548	endif
	1549
[880]	1550	if (cursor == 1) then
	1551	print,'--> Ls =',zls180./pi
	1552	endif
[253]	1553	! -------------------------------------------------------------------
	1554	! Writing NetCDF file "RESTARTFI" at the end of the run
	1555	! -------------------------------------------------------------------
	1556	! Note: 'restartfi' is stored just before dynamics are stored
	1557	! in 'restart'. Between now and the writting of 'restart',
	1558	! there will have been the itau=itau+1 instruction and
	1559	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	1560	! thus we store for time=time+dtvr
	1561
	1562	if(lastcall) then
	1563	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
	1564
	1565
	1566	! Update surface ice distribution to iterate to steady state if requested
	1567	if(ice_update)then
[305]	1568
[787]	1569	do ig=1,ngrid
[253]	1570
[305]	1571	delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
	1572
[365]	1573	! add multiple years of evolution
[728]	1574	qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep
[305]	1575
	1576	! if ice has gone -ve, set to zero
	1577	if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
	1578	qsurf_hist(ig,igcm_h2o_ice) = 0.0
[253]	1579	endif
[305]	1580
[365]	1581	! if ice is seasonal, set to zero (NEW)
	1582	if(ice_min(ig).lt.0.01)then
	1583	qsurf_hist(ig,igcm_h2o_ice) = 0.0
[253]	1584	endif
	1585
	1586	enddo
[305]	1587
	1588	! enforce ice conservation
[728]	1589	ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*area(:) )
	1590	qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
[305]	1591
[253]	1592	endif
	1593
	1594	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
[907]	1595	#ifdef CPP_PARA
	1596	! for now in parallel we use a different routine to write restartfi.nc
	1597	call phyredem(ngrid,"restartfi.nc", &
	1598	ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
	1599	cloudfrac,totcloudfrac,hice)
	1600	#else
[787]	1601	call physdem1(ngrid,"restartfi.nc",long,lati,nsoilmx,nq, &
[253]	1602	ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
	1603	area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
[787]	1604	cloudfrac,totcloudfrac,hice,noms)
[907]	1605	#endif
[253]	1606	endif
	1607
[861]	1608
[253]	1609	! -----------------------------------------------------------------
	1610	! Saving statistics :
	1611	! -----------------------------------------------------------------
	1612	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
	1613	! which can later be used to make the statistic files of the run:
	1614	! "stats") only possible in 3D runs !
	1615
	1616
	1617	if (callstats) then
	1618
	1619	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
	1620	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1621	call wstats(ngrid,"fluxsurf_lw", &
	1622	"Thermal IR radiative flux to surface","W.m-2",2, &
	1623	fluxsurf_lw)
	1624	! call wstats(ngrid,"fluxsurf_sw", &
	1625	! "Solar radiative flux to surface","W.m-2",2, &
	1626	! fluxsurf_sw_tot)
	1627	call wstats(ngrid,"fluxtop_lw", &
	1628	"Thermal IR radiative flux to space","W.m-2",2, &
	1629	fluxtop_lw)
	1630	! call wstats(ngrid,"fluxtop_sw", &
	1631	! "Solar radiative flux to space","W.m-2",2, &
	1632	! fluxtop_sw_tot)
[526]	1633
	1634	call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
	1635	call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
	1636	call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
	1637
[253]	1638	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
	1639	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
	1640	call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
	1641	call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
	1642	call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
	1643
	1644	if (tracer) then
[526]	1645	do iq=1,nq
	1646	call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
[787]	1647	call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
[526]	1648	'kg m^-2',2,qsurf(1,iq) )
	1649
[787]	1650	call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
[526]	1651	'kg m^-2',2,qcol(1,iq) )
[787]	1652	! call wstats(ngrid,trim(noms(iq))//'_reff', &
[726]	1653	! trim(noms(iq))//'_reff', &
	1654	! 'm',3,reffrad(1,1,iq))
[526]	1655	end do
[253]	1656	if (water) then
[787]	1657	vmr=zq(1:ngrid,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
[253]	1658	call wstats(ngrid,"vmr_h2ovapor", &
	1659	"H2O vapour volume mixing ratio","mol/mol", &
	1660	3,vmr)
	1661	endif ! of if (water)
	1662
	1663	endif !tracer
	1664
	1665	if(lastcall) then
	1666	write (,) "Writing stats..."
	1667	call mkstats(ierr)
	1668	endif
	1669	endif !if callstats
	1670
	1671
	1672	! ----------------------------------------------------------------------
	1673	! output in netcdf file "DIAGFI", containing any variable for diagnostic
	1674	! (output with period "ecritphy", set in "run.def")
	1675	! ----------------------------------------------------------------------
	1676	! writediagfi can also be called from any other subroutine for any variable.
	1677	! but its preferable to keep all the calls in one place...
	1678
	1679	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1680	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
	1681	call writediagfi(ngrid,"temp","temperature","K",3,zt)
[597]	1682	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
[253]	1683	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1684	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
	1685	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
[731]	1686	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
[253]	1687
[965]	1688	! Subsurface temperatures
[969]	1689	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1690	! call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
[965]	1691
[253]	1692	! Total energy balance diagnostics
	1693	if(callrad.and.(.not.newtonian))then
[731]	1694	call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo)
[253]	1695	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
	1696	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
	1697	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
	1698	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
	1699	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
	1700	endif
[594]	1701
	1702	if(enertest) then
[622]	1703	if (calldifv) then
[787]	1704	call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
	1705	! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
	1706	! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
	1707	! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
	1708	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
[622]	1709	endif
[596]	1710	if (corrk) then
[918]	1711	call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
	1712	call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
[596]	1713	endif
[594]	1714	if(watercond) then
	1715	call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
[622]	1716	call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
[787]	1717	call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
	1718	! call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
	1719	! call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
[594]	1720	endif
	1721	endif
[253]	1722
	1723	! Temporary inclusions for heating diagnostics
	1724	! call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
[863]	1725	call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
	1726	call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
	1727	call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
[253]	1728
	1729	! debugging
[368]	1730	!call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
[253]	1731	!call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
	1732
	1733	! Output aerosols
[728]	1734	if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
[787]	1735	call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
[728]	1736	if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
[787]	1737	call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
[728]	1738	if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
[787]	1739	call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
[728]	1740	if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
[787]	1741	call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
[253]	1742
	1743	! Output tracers
	1744	if (tracer) then
	1745	do iq=1,nq
[368]	1746	call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
[787]	1747	! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
[253]	1748	! 'kg m^-2',2,qsurf(1,iq) )
[787]	1749	call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
[253]	1750	'kg m^-2',2,qsurf_hist(1,iq) )
[787]	1751	call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
[253]	1752	'kg m^-2',2,qcol(1,iq) )
	1753
[594]	1754	if(watercond.or.CLFvarying)then
[918]	1755	call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
	1756	call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
	1757	call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
[253]	1758	call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
	1759	endif
	1760
	1761	if(waterrain)then
[787]	1762	call writediagfi(ngrid,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
	1763	call writediagfi(ngrid,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
[253]	1764	endif
	1765
	1766	if(hydrology)then
[787]	1767	call writediagfi(ngrid,"hice","oceanic ice height","m",2,hice)
[253]	1768	endif
	1769
	1770	if(ice_update)then
[787]	1771	call writediagfi(ngrid,"ice_min","min annual ice","m",2,ice_min)
	1772	call writediagfi(ngrid,"ice_ini","initial annual ice","m",2,ice_initial)
[253]	1773	endif
	1774
	1775	do ig=1,ngrid
	1776	if(tau_col(ig).gt.1.e3)then
[959]	1777	! print*,'WARNING: tau_col=',tau_col(ig)
	1778	! print*,'at ig=',ig,'in PHYSIQ'
[253]	1779	endif
	1780	end do
	1781
[787]	1782	call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
[253]	1783
	1784	enddo
	1785	endif
	1786
[526]	1787	! output spectrum
	1788	if(specOLR.and.corrk)then
[728]	1789	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
	1790	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
[526]	1791	endif
[253]	1792
	1793	icount=icount+1
	1794
[471]	1795	if (lastcall) then
[787]	1796
	1797	! deallocate gas variables
[716]	1798	IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )
	1799	IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac ) ! both allocated in su_gases.F90
[787]	1800
	1801	! deallocate saved arrays
	1802	! this is probably not that necessary
	1803	! ... but probably a good idea to clean the place before we leave
	1804	IF ( ALLOCATED(tsurf)) DEALLOCATE(tsurf)
	1805	IF ( ALLOCATED(tsoil)) DEALLOCATE(tsoil)
	1806	IF ( ALLOCATED(albedo)) DEALLOCATE(albedo)
	1807	IF ( ALLOCATED(albedo0)) DEALLOCATE(albedo0)
	1808	IF ( ALLOCATED(rnat)) DEALLOCATE(rnat)
	1809	IF ( ALLOCATED(emis)) DEALLOCATE(emis)
	1810	IF ( ALLOCATED(dtrad)) DEALLOCATE(dtrad)
	1811	IF ( ALLOCATED(fluxrad_sky)) DEALLOCATE(fluxrad_sky)
	1812	IF ( ALLOCATED(fluxrad)) DEALLOCATE(fluxrad)
	1813	IF ( ALLOCATED(capcal)) DEALLOCATE(capcal)
	1814	IF ( ALLOCATED(fluxgrd)) DEALLOCATE(fluxgrd)
	1815	IF ( ALLOCATED(qsurf)) DEALLOCATE(qsurf)
	1816	IF ( ALLOCATED(q2)) DEALLOCATE(q2)
	1817	IF ( ALLOCATED(ztprevious)) DEALLOCATE(ztprevious)
	1818	IF ( ALLOCATED(cloudfrac)) DEALLOCATE(cloudfrac)
	1819	IF ( ALLOCATED(totcloudfrac)) DEALLOCATE(totcloudfrac)
	1820	IF ( ALLOCATED(qsurf_hist)) DEALLOCATE(qsurf_hist)
	1821	IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
	1822	IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
	1823	IF ( ALLOCATED(ice_initial)) DEALLOCATE(ice_initial)
	1824	IF ( ALLOCATED(ice_min)) DEALLOCATE(ice_min)
	1825
	1826	IF ( ALLOCATED(fluxsurf_lw)) DEALLOCATE(fluxsurf_lw)
	1827	IF ( ALLOCATED(fluxsurf_sw)) DEALLOCATE(fluxsurf_sw)
	1828	IF ( ALLOCATED(fluxtop_lw)) DEALLOCATE(fluxtop_lw)
	1829	IF ( ALLOCATED(fluxabs_sw)) DEALLOCATE(fluxabs_sw)
	1830	IF ( ALLOCATED(fluxtop_dn)) DEALLOCATE(fluxtop_dn)
	1831	IF ( ALLOCATED(fluxdyn)) DEALLOCATE(fluxdyn)
	1832	IF ( ALLOCATED(OLR_nu)) DEALLOCATE(OLR_nu)
	1833	IF ( ALLOCATED(OSR_nu)) DEALLOCATE(OSR_nu)
	1834	IF ( ALLOCATED(sensibFlux)) DEALLOCATE(sensibFlux)
	1835	IF ( ALLOCATED(zdtlw)) DEALLOCATE(zdtlw)
	1836	IF ( ALLOCATED(zdtsw)) DEALLOCATE(zdtsw)
	1837	IF ( ALLOCATED(tau_col)) DEALLOCATE(tau_col)
	1838
	1839	!! this is defined in comsaison_h
	1840	IF ( ALLOCATED(mu0)) DEALLOCATE(mu0)
	1841	IF ( ALLOCATED(fract)) DEALLOCATE(fract)
	1842
	1843	!! this is defined in comsoil_h
	1844	IF ( ALLOCATED(layer)) DEALLOCATE(layer)
	1845	IF ( ALLOCATED(mlayer)) DEALLOCATE(mlayer)
	1846	IF ( ALLOCATED(inertiedat)) DEALLOCATE(inertiedat)
	1847
	1848	!! this is defined in surfdat_h
	1849	IF ( ALLOCATED(albedodat)) DEALLOCATE(albedodat)
	1850	IF ( ALLOCATED(phisfi)) DEALLOCATE(phisfi)
	1851	IF ( ALLOCATED(zmea)) DEALLOCATE(zmea)
	1852	IF ( ALLOCATED(zstd)) DEALLOCATE(zstd)
	1853	IF ( ALLOCATED(zsig)) DEALLOCATE(zsig)
	1854	IF ( ALLOCATED(zgam)) DEALLOCATE(zgam)
	1855	IF ( ALLOCATED(zthe)) DEALLOCATE(zthe)
	1856	IF ( ALLOCATED(dryness)) DEALLOCATE(dryness)
	1857	IF ( ALLOCATED(watercaptag)) DEALLOCATE(watercaptag)
	1858
	1859	!! this is defined in tracer_h
	1860	IF ( ALLOCATED(noms)) DEALLOCATE(noms)
	1861	IF ( ALLOCATED(mmol)) DEALLOCATE(mmol)
	1862	IF ( ALLOCATED(radius)) DEALLOCATE(radius)
	1863	IF ( ALLOCATED(rho_q)) DEALLOCATE(rho_q)
	1864	IF ( ALLOCATED(qext)) DEALLOCATE(qext)
	1865	IF ( ALLOCATED(alpha_lift)) DEALLOCATE(alpha_lift)
	1866	IF ( ALLOCATED(alpha_devil)) DEALLOCATE(alpha_devil)
	1867	IF ( ALLOCATED(qextrhor)) DEALLOCATE(qextrhor)
	1868	IF ( ALLOCATED(igcm_dustbin)) DEALLOCATE(igcm_dustbin)
	1869
	1870	!! this is defined in comgeomfi_h
	1871	IF ( ALLOCATED(lati)) DEALLOCATE(lati)
	1872	IF ( ALLOCATED(long)) DEALLOCATE(long)
	1873	IF ( ALLOCATED(area)) DEALLOCATE(area)
	1874	IF ( ALLOCATED(sinlat)) DEALLOCATE(sinlat)
	1875	IF ( ALLOCATED(coslat)) DEALLOCATE(coslat)
	1876	IF ( ALLOCATED(sinlon)) DEALLOCATE(sinlon)
	1877	IF ( ALLOCATED(coslon)) DEALLOCATE(coslon)
	1878
[471]	1879	endif
	1880
[253]	1881	return
[751]	1882	end subroutine physiq

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

Download in other formats:

Original Format