Context Navigation

physiq.F90 @ 1115

Last change on this file since 1115 was 1016, checked in by jleconte, 12 years ago

07/08/2013 == JL

Water cycle in double precision (largescale+moistadj)
Improved wate rayleigh.
First step for rayleigh with variable species. Now, just need to change optcv.
changed some interpolation indices in callcorrk to limit dependency of OLR on the number of layers

Property svn:executable set to *

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

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

Download in other formats:

Original Format