Context Navigation

physiq.F90 @ 1243

Last change on this file since 1243 was 1216, checked in by emillour, 11 years ago

Generic model:
Major cleanup, in order to ease the use of LMDZ.GENERIC with (parallel) dynamics
in LMDZ.COMMON: (NB: this will break LMDZ.UNIVERSAL, which should be thrashed
in the near future)

Updated makegcm_* scripts (and makdim) and added the "-full" (to enforce full recomputation of the model) option
In dyn3d: converted control.h to module control_mod.F90 and converted iniadvtrac.F to module infotrac.F90
Added module mod_const_mpi.F90 in dyn3d (not used in serial mode)
Rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallelism)
added created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the max and min of a field over the whole planet. This should be further imroved with computation of means (possibly area weighed), etc.

Property svn:executable set to *

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

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

Download in other formats:

Original Format