Context Navigation

physiq.F90 @ 1202

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

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

Download in other formats:

Original Format