Context Navigation

physiq.F90 @ 739

Last change on this file since 739 was 731, checked in by jleconte, 13 years ago
Corrected precipitation evaporation scheme + snow fall (JL+BC) Commented some unecessary writediagfi calls that where used for test.
Property svn:executable set to ``*
File size: 70.6 KB

Rev	Line
[597]	1	subroutine physiq(ngrid,nlayer,nq, &
[253]	2	firstcall,lastcall, &
	3	pday,ptime,ptimestep, &
	4	pplev,pplay,pphi, &
	5	pu,pv,pt,pq, &
	6	pw, &
	7	pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
	8
[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 :: reffrad1(ngridmx,nlayermx,naerkind) ! for clearsky call to callcorrk
	399	real,save :: nueffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius variance
	400	real :: reffh2oliq(ngridmx,nlayermx) ! liquid water particles effective radius (m)
	401	real :: reffh2oice(ngridmx,nlayermx) ! water ice particles effective radius (m)
[253]	402	real reffH2O(ngridmx,nlayermx)
[726]	403	real reffcol(ngridmx,naerkind)
[253]	404
	405	! included by RW for sourceevol
[728]	406	real, save :: ice_initial(ngridmx)!, isoil
[305]	407	real delta_ice,ice_tot
[728]	408	real, save :: ice_min(ngridmx)
	409
[253]	410	integer num_run
[728]	411	logical,save :: ice_update
	412
[253]	413	!=======================================================================
	414
	415
	416	! 1. Initialisation
	417	! -----------------
	418
	419	! 1.1 Initialisation only at first call
	420	! ---------------------------------------
	421	if (firstcall) then
	422
	423
	424	! variables set to 0
	425	! ~~~~~~~~~~~~~~~~~~
	426	dtrad(:,:) = 0.0
	427	fluxrad(:) = 0.0
	428	tau_col(:) = 0.0
	429	zdtsw(:,:) = 0.0
	430	zdtlw(:,:) = 0.0
[726]	431
	432	! initialize aerosol indexes
	433	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	434	call iniaerosol()
	435
[253]	436
	437	! initialize tracer names, indexes and properties
	438	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	439	tracerdyn=tracer
	440	if (tracer) then
	441	call initracer()
	442	endif ! end tracer
	443
[726]	444	!
	445
[253]	446	! read startfi (initial parameters)
	447	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	448	call phyetat0("startfi.nc",0,0,nsoilmx,nq, &
	449	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, &
	450	cloudfrac,totcloudfrac,hice)
	451
	452	if (pday.ne.day_ini) then
	453	write(,) "ERROR in physiq.F90:"
	454	write(,) "bad synchronization between physics and dynamics"
	455	write(,) "dynamics day: ",pday
	456	write(,) "physics day: ",day_ini
	457	stop
	458	endif
	459
	460	write (,) 'In physiq day_ini =', day_ini
	461
	462	! Initialize albedo and orbital calculation
	463	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	464	call surfini(ngrid,qsurf,albedo0)
	465	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
	466
[728]	467	albedo(:)=albedo0(:)
[253]	468
	469	if(tlocked)then
	470	print*,'Planet is tidally locked at resonance n=',nres
	471	print*,'Make sure you have the right rotation rate!!!'
	472	endif
	473
	474	! initialize soil
	475	! ~~~~~~~~~~~~~~~
	476	if (callsoil) then
	477	call soil(ngrid,nsoilmx,firstcall,inertiedat, &
	478	ptimestep,tsurf,tsoil,capcal,fluxgrd)
	479	else
	480	print*,'WARNING! Thermal conduction in the soil turned off'
[728]	481	capcal(:)=1.e6
	482	fluxgrd(:)=0.
	483	if(noradsurf)then
	484	fluxgrd(:)=10.0
	485	endif
[253]	486	print*,'Flux from ground = ',fluxgrd,' W m^-2'
	487	endif
	488	icount=1
	489
	490	! decide whether to update ice at end of run
	491	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	492	ice_update=.false.
	493	if(sourceevol)then
	494	open(128,file='num_run',form='formatted')
	495	read(128,*) num_run
	496	close(128)
	497
[365]	498	if(num_run.ne.0.and.mod(num_run,2).eq.0)then
	499	!if(num_run.ne.0.and.mod(num_run,3).eq.0)then
[253]	500	print*,'Updating ice at end of this year!'
	501	ice_update=.true.
	502	ice_min(:)=1.e4
	503	endif
	504	endif
	505
	506	! define surface as continent or ocean
	507	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[728]	508	rnat(:)=1
[253]	509	do ig=1,ngridmx
	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	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[728]	522	qsurf_hist(:,:)=qsurf(:,:)
[253]	523
	524	! initialise variable for dynamical heating diagnostic
	525	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	526	ztprevious(:,:)=pt(:,:)
	527
	528	! Set temperature just above condensation temperature (for Early Mars)
	529	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	530	if(nearco2cond) then
	531	write(,)' WARNING! Starting at Tcond+1K'
	532	do l=1, nlayer
	533	do ig=1,ngrid
	534	pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4 &
	535	-pt(ig,l)) / ptimestep
	536	enddo
	537	enddo
	538	endif
	539
	540	if(meanOLR)then
	541	! to record global radiative balance
	542	call system('rm -f rad_bal.out')
	543	! to record global mean/max/min temperatures
	544	call system('rm -f tem_bal.out')
	545	! to record global hydrological balance
	546	call system('rm -f h2o_bal.out')
	547	endif
	548
	549	watertest=.false.
	550	if(water)then
	551	! initialise variables for water cycle
	552
[365]	553	if(enertest)then
	554	watertest = .true.
	555	endif
	556
[728]	557	if(ice_update)then
	558	ice_initial(:)=qsurf(:,igcm_h2o_ice)
	559	endif
[253]	560
	561	endif
	562	call su_watercycle ! even if we don't have a water cycle, we might
	563	! need epsi for the wvp definitions in callcorrk.F
	564
	565	endif ! (end of "if firstcall")
	566
	567	! ---------------------------------------------------
	568	! 1.2 Initializations done at every physical timestep:
	569	! ---------------------------------------------------
	570
	571	if (ngrid.NE.ngridmx) then
	572	print*,'STOP in PHYSIQ'
	573	print*,'Probleme de dimensions :'
	574	print*,'ngrid = ',ngrid
	575	print*,'ngridmx = ',ngridmx
	576	stop
	577	endif
	578
	579	! Initialize various variables
	580	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	581
[728]	582	pdu(1:ngridmx,1:nlayermx) = 0.0
	583	pdv(1:ngridmx,1:nlayermx) = 0.0
[253]	584	if ( .not.nearco2cond ) then
[728]	585	pdt(1:ngridmx,1:nlayermx) = 0.0
	586	endif
	587	pdq(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
	588	pdpsrf(1:ngridmx) = 0.0
	589	zflubid(1:ngridmx) = 0.0
	590	zdtsurf(1:ngridmx) = 0.0
	591	dqsurf(1:ngridmx,1:nqmx)= 0.0
[253]	592
	593	zday=pday+ptime ! compute time, in sols (and fraction thereof)
	594
	595	! Compute Stellar Longitude (Ls)
	596	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	597	if (season) then
	598	call stellarlong(zday,zls)
	599	else
	600	call stellarlong(float(day_ini),zls)
	601	end if
	602
	603	! Compute geopotential between layers
	604	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	605
[728]	606	zzlay(1:ngridmx,1:nlayermx)=pphi(1:ngridmx,1:nlayermx)/g
	607
	608	zzlev(1:ngridmx,1)=0.
	609	zzlev(1:ngridmx,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
	610
[253]	611	do l=2,nlayer
	612	do ig=1,ngrid
	613	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
	614	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
	615	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	616	enddo
	617	enddo
	618	! Potential temperature calculation may not be the same in physiq and dynamic...
	619
	620	! Compute potential temperature
	621	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	622
[597]	623	do l=1,nlayer
[253]	624	do ig=1,ngrid
	625	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
[597]	626	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
[651]	627	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/g
	628	massarea(ig,l)=mass(ig,l)*area(ig)
[253]	629	enddo
	630	enddo
	631
	632	!-----------------------------------------------------------------------
	633	! 2. Compute radiative tendencies
	634	!-----------------------------------------------------------------------
	635
	636	if (callrad) then
[526]	637	if( mod(icount-1,iradia).eq.0.or.lastcall) then
[253]	638
	639	! Compute local stellar zenith angles
	640	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	641	call orbite(zls,dist_star,declin)
	642
	643	if (tlocked) then
	644	ztim1=SIN(declin)
	645	ztim2=COS(declin)COS(2.pi(zday/year_day) - zlsnres)
	646	ztim3=-COS(declin)SIN(2.pi(zday/year_day) - zlsnres)
	647
[728]	648	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
[253]	649	ztim1,ztim2,ztim3,mu0,fract)
	650
	651	elseif (diurnal) then
	652	ztim1=SIN(declin)
	653	ztim2=COS(declin)COS(2.pi*(zday-.5))
	654	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
	655
	656	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
	657	ztim1,ztim2,ztim3,mu0,fract)
	658
	659	else
	660
	661	call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
	662	! WARNING: this function appears not to work in 1D
	663
	664	endif
	665
	666	if (corrk) then
	667
	668	! a) Call correlated-k radiative transfer scheme
	669	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	670
	671	if(kastprof)then
[305]	672	print*,'kastprof should not = true here'
	673	call abort
[253]	674	endif
[538]	675	muvar(:,:)=0.0 ! only used for climate evolution studies (kcm1d) for now
	676
[526]	677	! standard callcorrk
	678	clearsky=.false.
[538]	679	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
[526]	680	albedo,emis,mu0,pplev,pplay,pt, &
[586]	681	tsurf,fract,dist_star,aerosol,muvar, &
[526]	682	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, &
[538]	683	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, &
[526]	684	reffrad,tau_col,cloudfrac,totcloudfrac, &
	685	clearsky,firstcall,lastcall)
[253]	686
[526]	687	! Option to call scheme once more for clear regions
[253]	688	if(CLFvarying)then
	689
[716]	690	! ---> PROBLEMS WITH ALLOCATED ARRAYS
	691	! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
[253]	692	clearsky=.true.
[538]	693	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
[253]	694	albedo,emis,mu0,pplev,pplay,pt, &
[586]	695	tsurf,fract,dist_star,aerosol,muvar, &
[253]	696	zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
[526]	697	fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1, &
[728]	698	reffrad1,tau_col1,cloudfrac,totcloudfrac, &
[538]	699	clearsky,firstcall,lastcall)
[716]	700	clearsky = .false. ! just in case.
[253]	701
	702	! Sum the fluxes and heating rates from cloudy/clear cases
	703	do ig=1,ngrid
	704	tf=totcloudfrac(ig)
	705	ntf=1.-tf
	706
[526]	707	fluxsurf_lw(ig) = ntffluxsurf_lw1(ig) + tffluxsurf_lw(ig)
	708	fluxsurf_sw(ig) = ntffluxsurf_sw1(ig) + tffluxsurf_sw(ig)
	709	fluxtop_lw(ig) = ntffluxtop_lw1(ig) + tffluxtop_lw(ig)
	710	fluxabs_sw(ig) = ntffluxabs_sw1(ig) + tffluxabs_sw(ig)
	711	tau_col(ig) = ntftau_col1(ig) + tftau_col(ig)
[253]	712
[728]	713	zdtlw(ig,1:nlayermx) = ntfzdtlw1(ig,1:nlayermx) + tfzdtlw(ig,1:nlayermx)
	714	zdtsw(ig,1:nlayermx) = ntfzdtsw1(ig,1:nlayermx) + tfzdtsw(ig,1:nlayermx)
[253]	715
[728]	716	OSR_nu(ig,1:L_NSPECTV) = ntfOSR_nu1(ig,1:L_NSPECTV) + tfOSR_nu(ig,1:L_NSPECTV)
	717	OLR_nu(ig,1:L_NSPECTV) = ntfOLR_nu1(ig,1:L_NSPECTV) + tfOLR_nu(ig,1:L_NSPECTV)
[253]	718
[526]	719	enddo
[253]	720
[526]	721	endif !CLFvarying
[253]	722
	723	! Radiative flux from the sky absorbed by the surface (W.m-2)
	724	GSR=0.0
[728]	725	fluxrad_sky(1:ngridmx)=emis(1:ngridmx)fluxsurf_lw(1:ngridmx)+fluxsurf_sw(1:ngridmx)(1.-albedo(1:ngridmx))
[253]	726
[728]	727	if(noradsurf)then ! no lower surface; SW flux just disappears
	728	GSR = SUM(fluxsurf_sw(1:ngridmx)*area(1:ngridmx))/totarea
	729	fluxrad_sky(1:ngridmx)=emis(1:ngridmx)*fluxsurf_lw(1:ngridmx)
	730	print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
[253]	731	endif
	732
	733	! Net atmospheric radiative heating rate (K.s-1)
[728]	734	dtrad(1:ngridmx,1:nlayermx)=zdtsw(1:ngridmx,1:nlayermx)+zdtlw(1:ngridmx,1:nlayermx)
[253]	735
	736	elseif(newtonian)then
	737
	738	! b) Call Newtonian cooling scheme
	739	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	740	call newtrelax(mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
	741
[728]	742	zdtsurf(1:ngridmx) = +(pt(1:ngridmx,1)-tsurf(1:ngridmx))/ptimestep
[253]	743	! e.g. surface becomes proxy for 1st atmospheric layer ?
	744
	745	else
	746
	747	! c) Atmosphere has no radiative effect
	748	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[728]	749	fluxtop_dn(1:ngridmx) = fract(1:ngridmx)mu0(1:ngridmx)Fat1AU/dist_star**2
	750	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
	751	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
	752	endif
	753	fluxsurf_sw(1:ngridmx) = fluxtop_dn(1:ngridmx)
	754	fluxrad_sky(1:ngridmx) = fluxtop_dn(1:ngridmx)*(1.-albedo(1:ngridmx))
	755	fluxtop_lw(1:ngridmx) = emis(1:ngridmx)sigmatsurf(1:ngridmx)**4
	756	! radiation skips the atmosphere entirely
[253]	757
	758
[728]	759	dtrad(1:ngridmx,1:nlayermx)=0.0
	760	! hence no atmospheric radiative heating
	761
[253]	762	endif ! if corrk
	763
	764	endif ! of if(mod(icount-1,iradia).eq.0)
	765
	766
	767	! Transformation of the radiative tendencies
	768	! ------------------------------------------
	769
[728]	770	zplanck(1:ngridmx)=tsurf(1:ngridmx)*tsurf(1:ngridmx)
	771	zplanck(1:ngridmx)=emis(1:ngridmx)sigmazplanck(1:ngridmx)*zplanck(1:ngridmx)
	772	fluxrad(1:ngridmx)=fluxrad_sky(1:ngridmx)-zplanck(1:ngridmx)
	773	pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+dtrad(1:ngridmx,1:nlayermx)
[253]	774
	775	!-------------------------
	776	! test energy conservation
	777	if(enertest)then
[651]	778	dEtotSW = cppSUM(massarea(:,:)zdtsw(:,:))/totarea
	779	dEtotLW = cppSUM(massarea(:,:)zdtlw(:,:))/totarea
	780	dEtotsSW = SUM(fluxsurf_sw(:)(1.-albedo(:))area(:))/totarea
[726]	781	dEtotsLW = SUM((fluxsurf_lw(:)emis(:)-zplanck(:))area(:))/totarea
[651]	782	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
	783	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
[253]	784	print*,'---------------------------------------------------------------'
[594]	785	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
	786	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
	787	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
	788	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
	789	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
	790	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
[253]	791	endif
	792	!-------------------------
	793
	794	endif ! of if (callrad)
	795
	796	!-----------------------------------------------------------------------
	797	! 4. Vertical diffusion (turbulent mixing):
	798	! -----------------------------------------
	799
	800	if (calldifv) then
[526]	801
[728]	802	zflubid(1:ngridmx)=fluxrad(1:ngridmx)+fluxgrd(1:ngridmx)
[253]	803
[728]	804	zdum1(1:ngridmx,1:nlayermx)=0.0
	805	zdum2(1:ngridmx,1:nlayermx)=0.0
[253]	806
[594]	807
	808	!JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff
	809	if (UseTurbDiff) then
	810
	811	call turbdiff(ngrid,nlayer,nq,rnat, &
[253]	812	ptimestep,capcal,lwrite, &
	813	pplay,pplev,zzlay,zzlev,z0, &
[594]	814	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
	815	zdum1,zdum2,pdt,pdq,zflubid, &
	816	zdudif,zdvdif,zdtdif,zdtsdif, &
[728]	817	sensibFlux,q2,zdqdif,zdqevap,zdqsdif,lastcall)
[594]	818
	819	else
	820
[728]	821	zdh(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)/zpopsk(1:ngridmx,1:nlayermx)
[594]	822
	823	call vdifc(ngrid,nlayer,nq,rnat,zpopsk, &
	824	ptimestep,capcal,lwrite, &
	825	pplay,pplev,zzlay,zzlev,z0, &
[253]	826	pu,pv,zh,pq,tsurf,emis,qsurf, &
	827	zdum1,zdum2,zdh,pdq,zflubid, &
[594]	828	zdudif,zdvdif,zdhdif,zdtsdif, &
	829	sensibFlux,q2,zdqdif,zdqsdif,lastcall)
[253]	830
[728]	831	zdtdif(1:ngridmx,1:nlayermx)=zdhdif(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx) ! for diagnostic only
	832	zdqevap(1:ngridmx,1:nlayermx)=0.
[594]	833
	834	end if !turbdiff
	835
[728]	836	pdv(1:ngridmx,1:nlayermx)=pdv(1:ngridmx,1:nlayermx)+zdvdif(1:ngridmx,1:nlayermx)
	837	pdu(1:ngridmx,1:nlayermx)=pdu(1:ngridmx,1:nlayermx)+zdudif(1:ngridmx,1:nlayermx)
	838	pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+zdtdif(1:ngridmx,1:nlayermx)
	839	zdtsurf(1:ngridmx)=zdtsurf(1:ngridmx)+zdtsdif(1:ngridmx)
[253]	840	if (tracer) then
[728]	841	pdq(1:ngridmx,1:nlayermx,1:nqmx)=pdq(1:ngridmx,1:nlayermx,1:nqmx)+ zdqdif(1:ngridmx,1:nlayermx,1:nqmx)
	842	dqsurf(1:ngridmx,1:nqmx)=dqsurf(1:ngridmx,1:nqmx) + zdqsdif(1:ngridmx,1:nqmx)
[253]	843	end if ! of if (tracer)
	844
	845	!-------------------------
	846	! test energy conservation
	847	if(enertest)then
[651]	848	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
[253]	849	do ig = 1, ngrid
[651]	850	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
[622]	851	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
[253]	852	enddo
[651]	853	dEtot = SUM(dEdiff(:)*area(:))/totarea
	854	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
	855	dEtots = SUM(dEdiffs(:)*area(:))/totarea
	856	AtmToSurf_TurbFlux=SUM(sensibFlux(:)*area(:))/totarea
[597]	857	if (UseTurbDiff) then
	858	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
	859	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
	860	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
	861	else
	862	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
	863	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
	864	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
	865	end if
[594]	866	! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
	867	! but not given back elsewhere
[253]	868	endif
	869	!-------------------------
	870
	871	!-------------------------
	872	! test water conservation
	873	if(watertest.and.water)then
[651]	874	dWtot = SUM(massarea(:,:)zdqdif(:,:,igcm_h2o_vap))ptimestep/totarea
	875	dWtots = SUM(zdqsdif(:,igcm_h2o_vap)area(:))ptimestep/totarea
[253]	876	do ig = 1, ngrid
[651]	877	vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
	878	Enddo
	879	dWtot = dWtot + SUM(massarea(:,:)zdqdif(:,:,igcm_h2o_ice))ptimestep/totarea
	880	dWtots = dWtots + SUM(zdqsdif(:,igcm_h2o_ice)area(:))ptimestep/totarea
	881	do ig = 1, ngrid
	882	vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
	883	Enddo
	884	nconsMAX=MAXVAL(vdifcncons(:))
[253]	885
	886	print*,'---------------------------------------------------------------'
	887	print*,'In difv atmospheric water change =',dWtot,' kg m-2'
	888	print*,'In difv surface water change =',dWtots,' kg m-2'
	889	print*,'In difv non-cons factor =',dWtot+dWtots,' kg m-2'
	890	print*,'In difv MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
	891
	892	endif
	893	!-------------------------
	894
	895	else
	896
	897	if(.not.newtonian)then
	898
[728]	899	zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + (fluxrad(1:ngridmx) + fluxgrd(1:ngridmx))/capcal(1:ngridmx)
[253]	900
	901	endif
	902
	903	endif ! of if (calldifv)
	904
	905
	906	!-----------------------------------------------------------------------
	907	! 5. Dry convective adjustment:
	908	! -----------------------------
	909
	910	if(calladj) then
	911
[728]	912	zdh(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)/zpopsk(1:ngridmx,1:nlayermx)
	913	zduadj(1:ngridmx,1:nlayermx)=0.0
	914	zdvadj(1:ngridmx,1:nlayermx)=0.0
	915	zdhadj(1:ngridmx,1:nlayermx)=0.0
[253]	916
	917
[586]	918	call convadj(ngrid,nlayer,nq,ptimestep, &
	919	pplay,pplev,zpopsk, &
	920	pu,pv,zh,pq, &
	921	pdu,pdv,zdh,pdq, &
	922	zduadj,zdvadj,zdhadj, &
	923	zdqadj)
[253]	924
[728]	925	pdu(1:ngridmx,1:nlayermx) = pdu(1:ngridmx,1:nlayermx) + zduadj(1:ngridmx,1:nlayermx)
	926	pdv(1:ngridmx,1:nlayermx) = pdv(1:ngridmx,1:nlayermx) + zdvadj(1:ngridmx,1:nlayermx)
	927	pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx) + zdhadj(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx)
	928	zdtadj(1:ngridmx,1:nlayermx) = zdhadj(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx) ! for diagnostic only
	929
[253]	930	if(tracer) then
[728]	931	pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqadj(1:ngridmx,1:nlayermx,1:nqmx)
[253]	932	end if
	933
	934	!-------------------------
	935	! test energy conservation
	936	if(enertest)then
[651]	937	dEtot=cppSUM(massarea(:,:)zdtadj(:,:))/totarea
[594]	938	print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
[253]	939	endif
	940	!-------------------------
	941
	942	!-------------------------
	943	! test water conservation
	944	if(watertest)then
[651]	945	dWtot = SUM(massarea(:,:)zdqadj(:,:,igcm_h2o_vap))ptimestep/totarea
[253]	946	do ig = 1, ngrid
[651]	947	cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
	948	Enddo
	949	dWtot = dWtot + SUM(massarea(:,:)zdqadj(:,:,igcm_h2o_ice))ptimestep/totarea
	950	do ig = 1, ngrid
	951	cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
	952	Enddo
	953	nconsMAX=MAXVAL(cadjncons(:))
[253]	954
	955	print*,'In convadj atmospheric water change =',dWtot,' kg m-2'
[651]	956	print*,'In convadj MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
[253]	957	endif
	958	!-------------------------
	959
	960	endif ! of if(calladj)
	961
	962	!-----------------------------------------------------------------------
	963	! 6. Carbon dioxide condensation-sublimation:
	964	! -------------------------------------------
	965
	966	if (co2cond) then
	967	if (.not.tracer) then
	968	print*,'We need a CO2 ice tracer to condense CO2'
	969	call abort
	970	endif
[726]	971	print*, 'we are in co2cond!!!!!'
[305]	972	call condense_cloud(ngrid,nlayer,nq,ptimestep, &
[253]	973	capcal,pplay,pplev,tsurf,pt, &
	974	pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, &
	975	qsurf(1,igcm_co2_ice),albedo,emis, &
	976	zdtc,zdtsurfc,pdpsrf,zduc,zdvc, &
[586]	977	zdqc,reffrad)
[253]	978
	979
[728]	980	pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+zdtc(1:ngridmx,1:nlayermx)
	981	pdv(1:ngridmx,1:nlayermx)=pdv(1:ngridmx,1:nlayermx)+zdvc(1:ngridmx,1:nlayermx)
	982	pdu(1:ngridmx,1:nlayermx)=pdu(1:ngridmx,1:nlayermx)+zduc(1:ngridmx,1:nlayermx)
	983	zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurfc(1:ngridmx)
	984
	985	pdq(1:ngridmx,1:nlayermx,1:nqmx)=pdq(1:ngridmx,1:nlayermx,1:nqmx)+ zdqc(1:ngridmx,1:nlayermx,1:nqmx)
[253]	986	! Note: we do not add surface co2ice tendency
	987	! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
	988
	989	!-------------------------
	990	! test energy conservation
	991	if(enertest)then
[651]	992	dEtot = cppSUM(massarea(:,:)zdtc(:,:))/totarea
	993	dEtots = SUM(capcal(:)zdtsurfc(:)area(:))/totarea
[253]	994	print*,'In co2cloud atmospheric energy change =',dEtot,' W m-2'
	995	print*,'In co2cloud surface energy change =',dEtots,' W m-2'
	996	endif
	997	!-------------------------
	998
	999	endif ! of if (co2cond)
	1000
	1001
	1002	!-----------------------------------------------------------------------
	1003	! 7. Specific parameterizations for tracers
	1004	! -----------------------------------------
	1005
	1006	if (tracer) then
	1007
	1008	! 7a. Water and ice
	1009	! ---------------
	1010	if (water) then
	1011
	1012	! ----------------------------------------
	1013	! Water ice condensation in the atmosphere
	1014	! ----------------------------------------
[728]	1015	if(watercond.and.(RLVTT.gt.1.e-8))then
[253]	1016
[728]	1017	! ----------------
	1018	! Moist convection
	1019	! ----------------
[253]	1020	! Re-evaporate cloud water/ice
[728]	1021	! dqevap1(1:ngridmx,1:nlayermx)=0.
	1022	! dtevap1(1:ngridmx,1:nlayermx)=0.
	1023	! call evap(ptimestep,pt,pq,pdq,pdt,dqevap1,dtevap1,qevap,tevap)
	1024
	1025	dqmoist(1:ngridmx,1:nlayermx,1:nqmx)=0.
	1026	dtmoist(1:ngridmx,1:nlayermx)=0.
[253]	1027
[728]	1028	call moistadj(pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
[253]	1029
[728]	1030	pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) &
	1031	+dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_vap)
	1032	pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) =pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) &
	1033	+dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_ice)
	1034	pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtmoist(1:ngridmx,1:nlayermx)
	1035
	1036	! do l=1,nlayermx
	1037	! do ig=1,ngridmx
	1038	! if(isnan(dqmoist(ig,l,igcm_h2o_ice))) print*,'Nan in dqmoist ice, abort,ig,l',ig,l
	1039	! if(isnan(dqmoist(ig,l,igcm_h2o_vap))) print*,'Nan in dqmoist vap, abort,ig,l',ig,l
	1040	! if(isnan(dtmoist(ig,l))) print*,'Nan in dtmoist, abort,ig,l',ig,l
	1041	! if(isnan(dqmoist(ig,l,igcm_h2o_ice)).or.isnan(dqmoist(ig,l,igcm_h2o_vap)).or.isnan(dtmoist(ig,l))) STOP
	1042	! enddo
	1043	! enddo
	1044	! print,'after moist',dqmoist(185,1,igcm_h2o_ice),dqmoist(185,1,igcm_h2o_vap),dtmoist(185,1)86400.
	1045
[253]	1046	!-------------------------
	1047	! test energy conservation
	1048	if(enertest)then
[728]	1049	dEtot=cppSUM(massarea(:,:)dtmoist(:,:))/totarea
[253]	1050	do ig = 1, ngrid
[728]	1051	madjdE(ig) = cppSUM(mass(:,:)dtmoist(:,:))
[253]	1052	enddo
	1053	print*,'In moistadj atmospheric energy change =',dEtot,' W m-2'
[728]	1054	print,'In moistadj MAX atmospheric energy change =',MAXVAL(pdt(:,:))86400.,'K/day'
[651]	1055
	1056	! test energy conservation
	1057	dWtot = SUM(massarea(:,:)dqmoist(:,:,igcm_h2o_vap))ptimestep/totarea
	1058	dWtot = dWtot + SUM(massarea(:,:)dqmoist(:,:,igcm_h2o_ice))ptimestep/totarea
[253]	1059	print*,'In moistadj atmospheric water change =',dWtot,' kg m-2'
	1060	endif
	1061	!-------------------------
	1062
	1063
[728]	1064	! --------------------------------
	1065	! Large scale condensation/evaporation
	1066	! --------------------------------
[253]	1067
[728]	1068	call largescale(ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
[253]	1069
[728]	1070	pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtlscale(1:ngridmx,1:nlayermx)
	1071	pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap)+dqvaplscale(1:ngridmx,1:nlayermx)
	1072	pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice)+dqcldlscale(1:ngridmx,1:nlayermx)
[253]	1073
	1074	!-------------------------
	1075	! test energy conservation
	1076	if(enertest)then
	1077	do ig = 1, ngrid
[728]	1078	lscaledE(ig) = cppSUM(mass(:,:)dtlscale(:,:))
[253]	1079	enddo
[728]	1080	dEtot=cppSUM(massarea(:,:)(dtlscale(:,:)))/totarea
	1081	if(isnan(dEtot)) then
	1082	print*,'Nan in largescale, abort'
	1083	STOP
	1084	endif
	1085	print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
	1086
	1087	! test water conservation
	1088	dWtot = SUM(massarea(:,:)dqvaplscale(:,:))ptimestep/totarea
	1089	dWtot = dWtot + SUM(massarea(:,:)dqcldlscale(:,:))ptimestep/totarea
	1090	print*,'In largescale atmospheric water change =',dWtot,' kg m-2'
[253]	1091	endif
	1092	!-------------------------
	1093
	1094	! compute cloud fraction
	1095	do l = 1, nlayer
	1096	do ig = 1,ngrid
	1097	cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
	1098	enddo
	1099	enddo
	1100
[728]	1101	endif ! of if (watercondense)
[253]	1102
	1103
	1104	! --------------------------------
	1105	! Water ice / liquid precipitation
	1106	! --------------------------------
[728]	1107	if(waterrain)then
[253]	1108
[728]	1109	zdqrain(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
	1110	zdqsrain(1:ngridmx) = 0.0
	1111	zdqssnow(1:ngridmx) = 0.0
[253]	1112
[728]	1113	call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq, &
[253]	1114	zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
	1115
[728]	1116	pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) &
	1117	+zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_vap)
	1118	pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) &
	1119	+zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_ice)
	1120	pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+zdtrain(1:ngridmx,1:nlayermx)
	1121	dqsurf(1:ngridmx,igcm_h2o_vap) = dqsurf(1:ngridmx,igcm_h2o_vap)+zdqsrain(1:ngridmx) ! a bug was here
	1122	dqsurf(1:ngridmx,igcm_h2o_ice) = dqsurf(1:ngridmx,igcm_h2o_ice)+zdqssnow(1:ngridmx)
	1123	rainout(1:ngridmx) = zdqsrain(1:ngridmx)+zdqssnow(1:ngridmx) ! diagnostic
[253]	1124
	1125
[651]	1126	!-------------------------
	1127	! test energy conservation
	1128	if(enertest)then
	1129	dEtot=cppSUM(massarea(:,:)zdtrain(:,:))/totarea
[253]	1130	print*,'In rain atmospheric T energy change =',dEtot,' W m-2'
[651]	1131	dItot = SUM(massarea(:,:)zdqrain(:,:,igcm_h2o_ice))/totareaRLVTT/cpp
	1132	dItot = dItot + SUM(area(:)zdqssnow(:))/totareaRLVTT/cpp
	1133	dVtot = SUM(massarea(:,:)zdqrain(:,:,igcm_h2o_vap))ptimestep/totarea
[728]	1134	dVtot = dVtot + SUM(area(:)zdqsrain(:))/totareaRLVTT/cpp
[651]	1135	dEtot = dItot + dVtot
	1136	print*,'In rain dItot =',dItot,' W m-2'
	1137	print*,'In rain dVtot =',dVtot,' W m-2'
[253]	1138	print*,'In rain atmospheric L energy change =',dEtot,' W m-2'
	1139
[651]	1140	! test water conservation
	1141	dWtot = SUM(massarea(:,:)zdqrain(:,:,igcm_h2o_vap))ptimestep/totarea
	1142	dWtot = dWtot + SUM(massarea(:,:)zdqrain(:,:,igcm_h2o_ice))ptimestep/totarea
	1143	dWtots = SUM((zdqsrain(:)+zdqssnow(:))area(:))ptimestep/totarea
[253]	1144	print*,'In rain atmospheric water change =',dWtot,' kg m-2'
	1145	print*,'In rain surface water change =',dWtots,' kg m-2'
	1146	print*,'In rain non-cons factor =',dWtot+dWtots,' kg m-2'
	1147	endif
	1148	!-------------------------
	1149
[728]	1150	end if ! of if (waterrain)
	1151	end if ! of if (water)
[253]	1152
	1153
	1154	! 7c. Aerosol particles
	1155	! -------------------
	1156	! -------------
	1157	! Sedimentation
	1158	! -------------
	1159	if (sedimentation) then
[728]	1160	zdqsed(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
	1161	zdqssed(1:ngridmx,1:nqmx) = 0.0
[253]	1162
	1163
	1164	!-------------------------
	1165	! find qtot
	1166	if(watertest)then
	1167	iq=3
[651]	1168	dWtot = SUM(massarea(:,:)pq(:,:,iq))ptimestep/totarea
	1169	dWtots = SUM(massarea(:,:)pdq(:,:,iq))ptimestep/totarea
[253]	1170	print*,'Before sedim pq =',dWtot,' kg m-2'
	1171	print*,'Before sedim pdq =',dWtots,' kg m-2'
	1172	endif
	1173	!-------------------------
	1174
	1175	call callsedim(ngrid,nlayer,ptimestep, &
	1176	pplev,zzlev,pt,pq,pdq,zdqsed,zdqssed,nq,reffH2O)
	1177
	1178	!-------------------------
	1179	! find qtot
	1180	if(watertest)then
	1181	iq=3
[651]	1182	dWtot = SUM(massarea(:,:)pq(:,:,iq))ptimestep/totarea
	1183	dWtots = SUM(massarea(:,:)pdq(:,:,iq))ptimestep/totarea
[253]	1184	print*,'After sedim pq =',dWtot,' kg m-2'
	1185	print*,'After sedim pdq =',dWtots,' kg m-2'
	1186	endif
	1187	!-------------------------
	1188
[728]	1189	! for now, we only allow H2O ice to sediment
[253]	1190	! and as in rain.F90, whether it falls as rain or snow depends
	1191	! only on the surface temperature
[728]	1192	pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqsed(1:ngridmx,1:nlayermx,1:nqmx)
	1193	dqsurf(1:ngridmx,1:nqmx) = dqsurf(1:ngridmx,1:nqmx) + zdqssed(1:ngridmx,1:nqmx)
[253]	1194
	1195	!-------------------------
	1196	! test water conservation
	1197	if(watertest)then
[651]	1198	dWtot = SUM(massarea(:,:)(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice)))ptimestep/totarea
	1199	dWtots = SUM((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))area(:))ptimestep/totarea
[253]	1200	print*,'In sedim atmospheric ice change =',dWtot,' kg m-2'
	1201	print*,'In sedim surface ice change =',dWtots,' kg m-2'
	1202	print*,'In sedim non-cons factor =',dWtot+dWtots,' kg m-2'
	1203	endif
	1204	!-------------------------
	1205
	1206	end if ! of if (sedimentation)
	1207
	1208
	1209	! 7d. Updates
	1210	! ---------
	1211
[728]	1212	! -----------------------------------
	1213	! Updating atm mass and tracer budget
	1214	! -----------------------------------
	1215
	1216	if(mass_redistrib) then
	1217
	1218	zdmassmr(1:ngridmx,1:nlayermx) = mass(1:ngridmx,1:nlayermx) * &
	1219	( zdqevap(1:ngridmx,1:nlayermx) &
	1220	+ zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_vap) &
	1221	+ dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_vap) &
	1222	+ dqvaplscale(1:ngridmx,1:nlayermx) )
	1223
	1224
	1225	call writediagfi(ngridmx,"mass_evap","mass gain"," ",3,zdmassmr)
	1226	call writediagfi(ngridmx,"mass","mass"," ",3,mass)
	1227
	1228	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
	1229	rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
	1230	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
	1231	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
	1232
	1233
	1234	pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqmr(1:ngridmx,1:nlayermx,1:nqmx)
	1235	dqsurf(1:ngridmx,1:nqmx) = dqsurf(1:ngridmx,1:nqmx) + zdqsurfmr(1:ngridmx,1:nqmx)
	1236	pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx) + zdtmr(1:ngridmx,1:nlayermx)
	1237	pdu(1:ngridmx,1:nlayermx) = pdu(1:ngridmx,1:nlayermx) + zdumr(1:ngridmx,1:nlayermx)
	1238	pdv(1:ngridmx,1:nlayermx) = pdv(1:ngridmx,1:nlayermx) + zdvmr(1:ngridmx,1:nlayermx)
	1239	pdpsrf(1:ngridmx) = pdpsrf(1:ngridmx) + zdpsrfmr(1:ngridmx)
	1240	zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurfmr(1:ngridmx)
	1241
	1242	! print,'after mass redistrib, q=',pq(211,1:nlayermx,igcm_h2o_vap)+ptimesteppdq(211,1:nlayermx,igcm_h2o_vap)
	1243	endif
	1244
	1245
	1246
[253]	1247	! ---------------------------------
	1248	! Updating tracer budget on surface
	1249	! ---------------------------------
	1250
	1251	if(hydrology)then
	1252
	1253	call hydrol(ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, &
	1254	capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice)
	1255	! note: for now, also changes albedo in the subroutine
	1256
[728]	1257	zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurf_hyd(1:ngridmx)
	1258	qsurf(1:ngridmx,1:nqmx) = qsurf(1:ngridmx,1:nqmx) + ptimestep*dqs_hyd(1:ngridmx,1:nqmx)
[253]	1259	! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
	1260
	1261	!-------------------------
	1262	! test energy conservation
	1263	if(enertest)then
[651]	1264	dEtots = SUM(area(:)capcal(:)zdtsurf_hyd(:))/totarea
	1265	print*,'In hydrol surface energy change =',dEtots,' W m-2'
[253]	1266	endif
	1267	!-------------------------
	1268
	1269	!-------------------------
	1270	! test water conservation
	1271	if(watertest)then
[651]	1272	dWtots = SUM(dqs_hyd(:,igcm_h2o_ice)area(:))ptimestep/totarea
[253]	1273	print*,'In hydrol surface ice change =',dWtots,' kg m-2'
[651]	1274	dWtots = SUM(dqs_hyd(:,igcm_h2o_vap)area(:))ptimestep/totarea
[253]	1275	print*,'In hydrol surface water change =',dWtots,' kg m-2'
	1276	print*,'---------------------------------------------------------------'
	1277	endif
	1278	!-------------------------
	1279
	1280	ELSE ! of if (hydrology)
	1281
[728]	1282	qsurf(1:ngridmx,1:nqmx)=qsurf(1:ngridmx,1:nqmx)+ptimestep*dqsurf(1:ngridmx,1:nqmx)
[253]	1283
	1284	END IF ! of if (hydrology)
	1285
	1286	! Add qsurf to qsurf_hist, which is what we save in
	1287	! diagfi.nc etc. At the same time, we set the water
	1288	! content of ocean gridpoints back to zero, in order
	1289	! to avoid rounding errors in vdifc, rain
[622]	1290	qsurf_hist(:,:) = qsurf(:,:)
[253]	1291
	1292	if(ice_update)then
[728]	1293	ice_min(1:ngridmx)=min(ice_min(1:ngridmx),qsurf(1:ngridmx,igcm_h2o_ice))
[253]	1294	endif
	1295
	1296	endif ! of if (tracer)
	1297
	1298	!-----------------------------------------------------------------------
	1299	! 9. Surface and sub-surface soil temperature
	1300	!-----------------------------------------------------------------------
	1301
	1302
	1303	! Increment surface temperature
[728]	1304	tsurf(1:ngridmx)=tsurf(1:ngridmx)+ptimestep*zdtsurf(1:ngridmx)
[253]	1305
	1306	! Compute soil temperatures and subsurface heat flux
	1307	if (callsoil) then
	1308	call soil(ngrid,nsoilmx,.false.,inertiedat, &
	1309	ptimestep,tsurf,tsoil,capcal,fluxgrd)
	1310	endif
	1311
	1312	!-------------------------
	1313	! test energy conservation
	1314	if(enertest)then
[728]	1315	dEtots = SUM(area(:)capcal(:)zdtsurf(:))/totarea
[597]	1316	print*,'Surface energy change =',dEtots,' W m-2'
[253]	1317	endif
	1318	!-------------------------
	1319
	1320	!-----------------------------------------------------------------------
	1321	! 10. Perform diagnostics and write output files
	1322	!-----------------------------------------------------------------------
	1323
	1324	! -------------------------------
	1325	! Dynamical fields incrementation
	1326	! -------------------------------
	1327	! For output only: the actual model integration is performed in the dynamics
	1328
	1329	! temperature, zonal and meridional wind
[728]	1330	zt(1:ngridmx,1:nlayermx) = pt(1:ngridmx,1:nlayermx) + pdt(1:ngridmx,1:nlayermx)*ptimestep
	1331	zu(1:ngridmx,1:nlayermx) = pu(1:ngridmx,1:nlayermx) + pdu(1:ngridmx,1:nlayermx)*ptimestep
	1332	zv(1:ngridmx,1:nlayermx) = pv(1:ngridmx,1:nlayermx) + pdv(1:ngridmx,1:nlayermx)*ptimestep
[253]	1333
[728]	1334	! diagnostic
	1335	zdtdyn(1:ngridmx,1:nlayermx) = pt(1:ngridmx,1:nlayermx)-ztprevious(1:ngridmx,1:nlayermx)
	1336	ztprevious(1:ngridmx,1:nlayermx) = zt(1:ngridmx,1:nlayermx)
[253]	1337
	1338	if(firstcall)then
[728]	1339	zdtdyn(1:ngridmx,1:nlayermx)=0.0
[253]	1340	endif
	1341
	1342	! dynamical heating diagnostic
	1343	do ig=1,ngrid
[728]	1344	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp/ptimestep
[253]	1345	enddo
	1346
	1347	! tracers
[728]	1348	zq(1:ngridmx,1:nlayermx,1:nqmx) = pq(1:ngridmx,1:nlayermx,1:nqmx) + pdq(1:ngridmx,1:nlayermx,1:nqmx)*ptimestep
[253]	1349
	1350	! surface pressure
[728]	1351	ps(1:ngridmx) = pplev(1:ngridmx,1) + pdpsrf(1:ngridmx)*ptimestep
[253]	1352
	1353	! pressure
	1354	do l=1,nlayer
[728]	1355	zplev(1:ngridmx,l) = pplev(1:ngridmx,l)/pplev(1:ngridmx,1)*ps(:)
	1356	zplay(1:ngridmx,l) = pplay(1:ngridmx,l)/pplev(1:ngridmx,1)*ps(1:ngridmx)
[253]	1357	enddo
	1358
	1359	! ---------------------------------------------------------
	1360	! Surface and soil temperature information
	1361	! ---------------------------------------------------------
	1362
[651]	1363	Ts1 = SUM(area(:)*tsurf(:))/totarea
	1364	Ts2 = MINVAL(tsurf(:))
	1365	Ts3 = MAXVAL(tsurf(:))
[253]	1366	if(callsoil)then
[651]	1367	TsS = SUM(area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
[253]	1368	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
	1369	print*,Ts1,Ts2,Ts3,TsS
	1370	else
	1371	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
	1372	print*,Ts1,Ts2,Ts3
	1373	endif
	1374
	1375	! ---------------------------------------------------------
	1376	! Check the energy balance of the simulation during the run
	1377	! ---------------------------------------------------------
	1378
	1379	if(corrk)then
	1380
[651]	1381	ISR = SUM(area(:)*fluxtop_dn(:))/totarea
	1382	ASR = SUM(area(:)*fluxabs_sw(:))/totarea
	1383	OLR = SUM(area(:)*fluxtop_lw(:))/totarea
	1384	GND = SUM(area(:)*fluxgrd(:))/totarea
	1385	DYN = SUM(area(:)*fluxdyn(:))/totarea
	1386	do ig=1,ngrid
[253]	1387	if(fluxtop_dn(ig).lt.0.0)then
	1388	print*,'fluxtop_dn has gone crazy'
	1389	print*,'fluxtop_dn=',fluxtop_dn(ig)
	1390	print*,'tau_col=',tau_col(ig)
	1391	print*,'aerosol=',aerosol(ig,:,:)
	1392	print*,'temp= ',pt(ig,:)
	1393	print*,'pplay= ',pplay(ig,:)
	1394	call abort
	1395	endif
	1396	end do
	1397
	1398	if(ngridmx.eq.1)then
	1399	DYN=0.0
	1400	endif
	1401
	1402	print*,' ISR ASR OLR GND DYN [W m^-2]'
[651]	1403	print*, ISR,ASR,OLR,GND,DYN
[253]	1404
	1405	if(enertest)then
[651]	1406	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
	1407	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
	1408	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
[253]	1409	endif
	1410
	1411	if(meanOLR)then
[526]	1412	if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
[253]	1413	! to record global radiative balance
[588]	1414	open(92,file="rad_bal.out",form='formatted',position='append')
[651]	1415	write(92,*) zday,ISR,ASR,OLR
[253]	1416	close(92)
[588]	1417	open(93,file="tem_bal.out",form='formatted',position='append')
[253]	1418	write(93,*) zday,Ts1,Ts2,Ts3,TsS
	1419	close(93)
	1420	endif
	1421	endif
	1422
	1423	endif
	1424
[651]	1425
[253]	1426	! ------------------------------------------------------------------
	1427	! Diagnostic to test radiative-convective timescales in code
	1428	! ------------------------------------------------------------------
	1429	if(testradtimes)then
[588]	1430	open(38,file="tau_phys.out",form='formatted',position='append')
[253]	1431	ig=1
	1432	do l=1,nlayer
	1433	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
	1434	enddo
	1435	close(38)
[726]	1436	print*,'As testradtimes enabled,'
	1437	print*,'exiting physics on first call'
[253]	1438	call abort
	1439	endif
	1440
	1441	! ---------------------------------------------------------
	1442	! Compute column amounts (kg m-2) if tracers are enabled
	1443	! ---------------------------------------------------------
	1444	if(tracer)then
[728]	1445	qcol(1:ngridmx,1:nqmx)=0.0
[253]	1446	do iq=1,nq
[728]	1447	do ig=1,ngridmx
	1448	qcol(ig,iq) = SUM( zq(ig,1:nlayermx,iq) * mass(ig,1:nlayermx))
	1449	enddo
[253]	1450	enddo
	1451
[726]	1452	! Generalised for arbitrary aerosols now. (LK)
[728]	1453	reffcol(1:ngridmx,1:naerkind)=0.0
	1454	if(co2cond.and.(iaero_co2.ne.0))then
	1455	call co2_reffrad(zq,reffrad)
	1456	do ig=1,ngridmx
	1457	reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayermx,igcm_co2_ice)reffrad(ig,1:nlayermx,iaero_co2)mass(ig,1:nlayermx))
[253]	1458	enddo
[728]	1459	endif
	1460	if(water.and.(iaero_h2o.ne.0))then
	1461	call h2o_reffrad(zq,zt,reffrad,nueffrad)
	1462	do ig=1,ngridmx
	1463	reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayermx,igcm_h2o_ice)reffrad(ig,1:nlayermx,iaero_h2o)mass(ig,1:nlayermx))
	1464	enddo
	1465	endif
[253]	1466
	1467	endif
	1468
	1469	! ---------------------------------------------------------
	1470	! Test for water conservation if water is enabled
	1471	! ---------------------------------------------------------
	1472
	1473	if(water)then
	1474
[651]	1475	icesrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_ice))/totarea
	1476	liqsrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_vap))/totarea
	1477	icecol = SUM(area(:)*qcol(:,igcm_h2o_ice))/totarea
	1478	vapcol = SUM(area(:)*qcol(:,igcm_h2o_vap))/totarea
[253]	1479
[651]	1480	h2otot = icesrf + liqsrf + icecol + vapcol
[253]	1481
[651]	1482	print*,' Total water amount [kg m^-2]: ',h2otot
[253]	1483	print*,' Surface ice Surface liq. Atmos. con. Atmos. vap. [kg m^-2] '
[651]	1484	print*, icesrf,liqsrf,icecol,vapcol
[253]	1485
	1486	if(meanOLR)then
[526]	1487	if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
[253]	1488	! to record global water balance
[588]	1489	open(98,file="h2o_bal.out",form='formatted',position='append')
[651]	1490	write(98,*) zday,icesrf,liqsrf,icecol,vapcol
[253]	1491	close(98)
	1492	endif
	1493	endif
	1494
	1495	endif
	1496
	1497	! ---------------------------------------------------------
	1498	! Calculate RH for diagnostic if water is enabled
	1499	! ---------------------------------------------------------
	1500
	1501	if(water)then
	1502	do l = 1, nlayer
	1503	do ig = 1, ngrid
[728]	1504	! call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
	1505	call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l))
	1506
[253]	1507	RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
	1508	enddo
	1509	enddo
	1510
	1511	! compute maximum possible H2O column amount (100% saturation)
	1512	do ig=1,ngrid
[728]	1513	H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
[253]	1514	enddo
	1515
	1516	endif
	1517
	1518
	1519	print*,''
	1520	print,'--> Ls =',zls180./pi
	1521	! -------------------------------------------------------------------
	1522	! Writing NetCDF file "RESTARTFI" at the end of the run
	1523	! -------------------------------------------------------------------
	1524	! Note: 'restartfi' is stored just before dynamics are stored
	1525	! in 'restart'. Between now and the writting of 'restart',
	1526	! there will have been the itau=itau+1 instruction and
	1527	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	1528	! thus we store for time=time+dtvr
	1529
	1530	if(lastcall) then
	1531	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
	1532
	1533
	1534	! Update surface ice distribution to iterate to steady state if requested
	1535	if(ice_update)then
[305]	1536
[253]	1537	do ig = 1, ngrid
	1538
[305]	1539	delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
	1540
[365]	1541	! add multiple years of evolution
[728]	1542	qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep
[305]	1543
	1544	! if ice has gone -ve, set to zero
	1545	if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
	1546	qsurf_hist(ig,igcm_h2o_ice) = 0.0
[253]	1547	endif
[305]	1548
[365]	1549	! if ice is seasonal, set to zero (NEW)
	1550	if(ice_min(ig).lt.0.01)then
	1551	qsurf_hist(ig,igcm_h2o_ice) = 0.0
[253]	1552	endif
	1553
	1554	enddo
[305]	1555
	1556	! enforce ice conservation
[728]	1557	ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*area(:) )
	1558	qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
[305]	1559
[253]	1560	endif
	1561
	1562	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
	1563	call physdem1("restartfi.nc",long,lati,nsoilmx,nq, &
	1564	ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
	1565	area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
	1566	cloudfrac,totcloudfrac,hice)
	1567	endif
	1568
	1569	! -----------------------------------------------------------------
	1570	! Saving statistics :
	1571	! -----------------------------------------------------------------
	1572	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
	1573	! which can later be used to make the statistic files of the run:
	1574	! "stats") only possible in 3D runs !
	1575
	1576
	1577	if (callstats) then
	1578
	1579	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
	1580	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1581	call wstats(ngrid,"fluxsurf_lw", &
	1582	"Thermal IR radiative flux to surface","W.m-2",2, &
	1583	fluxsurf_lw)
	1584	! call wstats(ngrid,"fluxsurf_sw", &
	1585	! "Solar radiative flux to surface","W.m-2",2, &
	1586	! fluxsurf_sw_tot)
	1587	call wstats(ngrid,"fluxtop_lw", &
	1588	"Thermal IR radiative flux to space","W.m-2",2, &
	1589	fluxtop_lw)
	1590	! call wstats(ngrid,"fluxtop_sw", &
	1591	! "Solar radiative flux to space","W.m-2",2, &
	1592	! fluxtop_sw_tot)
[526]	1593
	1594	call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
	1595	call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
	1596	call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
	1597
[253]	1598	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
	1599	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
	1600	call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
	1601	call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
	1602	call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
	1603
	1604	if (tracer) then
[526]	1605	do iq=1,nq
	1606	call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
	1607	call wstats(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
	1608	'kg m^-2',2,qsurf(1,iq) )
	1609
	1610	call wstats(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
	1611	'kg m^-2',2,qcol(1,iq) )
[726]	1612	! call wstats(ngridmx,trim(noms(iq))//'_reff', &
	1613	! trim(noms(iq))//'_reff', &
	1614	! 'm',3,reffrad(1,1,iq))
[526]	1615	end do
[253]	1616	if (water) then
	1617	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
	1618	call wstats(ngrid,"vmr_h2ovapor", &
	1619	"H2O vapour volume mixing ratio","mol/mol", &
	1620	3,vmr)
	1621	endif ! of if (water)
	1622
	1623	endif !tracer
	1624
	1625	if(lastcall) then
	1626	write (,) "Writing stats..."
	1627	call mkstats(ierr)
	1628	endif
	1629	endif !if callstats
	1630
	1631
	1632	! ----------------------------------------------------------------------
	1633	! output in netcdf file "DIAGFI", containing any variable for diagnostic
	1634	! (output with period "ecritphy", set in "run.def")
	1635	! ----------------------------------------------------------------------
	1636	! writediagfi can also be called from any other subroutine for any variable.
	1637	! but its preferable to keep all the calls in one place...
	1638
	1639	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1640	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
	1641	call writediagfi(ngrid,"temp","temperature","K",3,zt)
[597]	1642	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
[253]	1643	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1644	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
	1645	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
[731]	1646	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
[253]	1647
	1648	! Total energy balance diagnostics
	1649	if(callrad.and.(.not.newtonian))then
[731]	1650	call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo)
[253]	1651	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
	1652	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
	1653	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
	1654	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
	1655	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
	1656	endif
[594]	1657
	1658	if(enertest) then
[622]	1659	if (calldifv) then
	1660	call writediagfi(ngridmx,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
[731]	1661	! call writediagfi(ngridmx,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
	1662	! call writediagfi(ngridmx,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
	1663	! call writediagfi(ngridmx,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
[622]	1664	call writediagfi(ngridmx,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
	1665	endif
[596]	1666	if (corrk) then
[731]	1667	! call writediagfi(ngridmx,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
	1668	! call writediagfi(ngridmx,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
[596]	1669	endif
[594]	1670	if(watercond) then
	1671	call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
[622]	1672	call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
	1673	call writediagfi(ngridmx,"RH","relative humidity"," ",3,RH)
[731]	1674	! call writediagfi(ngridmx,"qsatatm","atm qsat"," ",3,qsat)
	1675	! call writediagfi(ngridmx,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
[594]	1676	endif
	1677	endif
[253]	1678
	1679	! Temporary inclusions for heating diagnostics
	1680	! call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
	1681	! call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
	1682	! call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
	1683	! call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
	1684
	1685	! debugging
[368]	1686	!call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
[253]	1687	!call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
	1688
	1689	! Output aerosols
[728]	1690	if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
	1691	call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
	1692	if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
	1693	call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
	1694	if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
	1695	call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
	1696	if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
	1697	call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
[253]	1698
	1699	! Output tracers
	1700	if (tracer) then
	1701	do iq=1,nq
[368]	1702	call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
[253]	1703	! call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
	1704	! 'kg m^-2',2,qsurf(1,iq) )
	1705	call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
	1706	'kg m^-2',2,qsurf_hist(1,iq) )
	1707	call writediagfi(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
	1708	'kg m^-2',2,qcol(1,iq) )
	1709
[594]	1710	if(watercond.or.CLFvarying)then
[731]	1711	! call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
	1712	! call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
	1713	! call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
[253]	1714	call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
	1715	endif
	1716
	1717	if(waterrain)then
	1718	call writediagfi(ngridmx,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
	1719	call writediagfi(ngridmx,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
	1720	endif
	1721
	1722	if(hydrology)then
	1723	call writediagfi(ngridmx,"hice","oceanic ice height","m",2,hice)
	1724	endif
	1725
	1726	if(ice_update)then
	1727	call writediagfi(ngridmx,"ice_min","min annual ice","m",2,ice_min)
	1728	call writediagfi(ngridmx,"ice_ini","initial annual ice","m",2,ice_initial)
	1729	endif
	1730
	1731	do ig=1,ngrid
	1732	if(tau_col(ig).gt.1.e3)then
	1733	print*,'WARNING: tau_col=',tau_col(ig)
	1734	print*,'at ig=',ig,'in PHYSIQ'
	1735	endif
	1736	end do
	1737
	1738	call writediagfi(ngridmx,"tau_col","Total aerosol optical depth","[]",2,tau_col)
	1739
	1740	enddo
	1741	endif
	1742
[526]	1743	! output spectrum
	1744	if(specOLR.and.corrk)then
[728]	1745	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
	1746	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
[526]	1747	endif
[253]	1748
	1749
	1750	icount=icount+1
	1751
[716]	1752	! deallocate gas variables
[471]	1753	if (lastcall) then
[716]	1754	IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )
	1755	IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac ) ! both allocated in su_gases.F90
[471]	1756	endif
	1757
	1758
[253]	1759	return
	1760	end subroutine physiq

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

Download in other formats:

Original Format