Context Navigation

physiq.F90 @ 527

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

13/02/2012 == JL + AS

All outputs are now in netCDF format. Even in 1D (No more G1D)
Clean up of the call to callcorrk when CLFvarying=true
Corrects a bug in writediagspecIR/VI. Output are now in W/m2/cm-1 as a function of the wavenumber in cm-1
Enable writediagspecIR/V to work in the CLFvarying=true case (output now done in Physiq after writediagfi)
Add a simple treatment for the supersaturation of CO2 (see forget et al 2012)
corrects a small bug when no clouds are present in aeropacity

Property svn:executable set to *

File size: 76.7 KB

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

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

Download in other formats:

Original Format