Context Navigation

physiq.F90 @ 724

Last change on this file since 724 was 716, checked in by rwordsworth, 13 years ago

Mainly updates to radiative transfer and gas management scheme.
Most CIA data now read from standard HITRAN datafiles. For the H2O
continuum, two options have been added: the standard CKD continuum,
and the empirical formula in PPC (Pierrehumbert 2010). Use the toggle
'H2Ocont_simple' in callphys.def to choose.

Note to Martians: I've changed the default values of 'sedimentation' and
'co2cond' in inifis to false. Both these are defined in the standard deftank
callphys.def file, so there shouldn't be any compatibility problems.

Property svn:executable set to *

File size: 69.5 KB

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

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

Download in other formats:

Original Format