Context Navigation

physiq.F90 @ 1149

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

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

Download in other formats:

Original Format