Context Navigation

physiq.F90 @ 879

Last change on this file since 879 was 867, checked in by aslmd, 13 years ago
LMDZ.GENERIC. Important note to developers: no more ngridmx otherwise parallel model is broken.
Property svn:executable set to ``*
File size: 76.0 KB

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

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

Download in other formats:

Original Format