Context Navigation

physiq.F90 @ 775

Last change on this file since 775 was 775, checked in by jleconte, 13 years ago

05/09/2012 == JL

Corrects typo in previous commit to account for possible eccentricity in

tlocked case

Property svn:executable set to *

File size: 70.0 KB

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

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

Download in other formats:

Original Format