Context Navigation

physiq.F90 @ 1371

Last change on this file since 1371 was 1346, checked in by emillour, 11 years ago

Mars and Generic GCM:

Bug correction on the computation of vertical velocity.

Property svn:executable set to *

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

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

Download in other formats:

Original Format