Context Navigation

physiq.F90 @ 1152

Last change on this file since 1152 was 1152, checked in by milmd, 11 years ago
LMDZ.GENERIC. added a condition to avoid a small bug in mucorr when running in 1D. but 'correction de rotondite' seems not generic at all in mucorr although the impact is small because hardcoded correction is pretty negligible (AS).
Property svn:executable set to ``*
File size: 78.5 KB

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

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

Download in other formats:

Original Format