Context Navigation

physiq.F90 @ 797

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

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
RUN RCM1D : perfect match in normal mode.
RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

Additional arguments to some subroutines (ngrid and nq)
F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this

Property svn:executable set to *

File size: 76.0 KB

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

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

Download in other formats:

Original Format