Context Navigation

physiq.F90 @ 862

Last change on this file since 862 was 862, checked in by aslmd, 13 years ago

LMDZ.UNIVERSAL
LMDZ.GENERIC

Added calls to Frederic Hourdin's subroutines to output physical fields in parallel. This is under precompiling flags CPP_PARA.
This allows for LMDZ.UNIVERSAL users to use writediagfi with parallel computations.
These lines are not compiled by casual users of LMDZ.GENERIC (or users of LMDZ.UNIVERSAL in sequential mode).

The precompiling flag NOWRITEDIAGFI is obsolete and has been deleted.

NOTES:

A better cleaner method might be proposed later
Added subroutines

A 0 LMDZ.UNIVERSAL/libf/phygeneric/iotd_ecrit.F90
A 0 LMDZ.UNIVERSAL/libf/phygeneric/iophys.F90
A 0 LMDZ.UNIVERSAL/libf/phygeneric/iotd.h
A 0 LMDZ.UNIVERSAL/libf/phygeneric/iotd_ini.F90
A 0 LMDZ.UNIVERSAL/libf/phygeneric/iotd_fin.F90
The iotd_* subroutines are actually supposed to be put in bibio in LMDZ.COMMON
This will be done later once agreed in the team

Property svn:executable set to *

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

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

Download in other formats:

Original Format