Context Navigation

physiq.F90 @ 1202

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

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

Download in other formats:

Original Format