Context Navigation

physiq.F90 @ 1243

Last change on this file since 1243 was 1216, checked in by emillour, 11 years ago

Generic model:
Major cleanup, in order to ease the use of LMDZ.GENERIC with (parallel) dynamics
in LMDZ.COMMON: (NB: this will break LMDZ.UNIVERSAL, which should be thrashed
in the near future)

Updated makegcm_* scripts (and makdim) and added the "-full" (to enforce full recomputation of the model) option
In dyn3d: converted control.h to module control_mod.F90 and converted iniadvtrac.F to module infotrac.F90
Added module mod_const_mpi.F90 in dyn3d (not used in serial mode)
Rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallelism)
added created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the max and min of a field over the whole planet. This should be further imroved with computation of means (possibly area weighed), etc.

Property svn:executable set to *

File size: 82.6 KB

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

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

Download in other formats:

Original Format