Context Navigation

physiq.F90 @ 799

Last change on this file since 799 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
RUN RCM1D : perfect match in normal mode.
RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

Additional arguments to some subroutines (ngrid and nq)
F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this

Property svn:executable set to *

File size: 76.0 KB

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

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

Download in other formats:

Original Format