Context Navigation

physiq.F90 @ 1255

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

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

Download in other formats:

Original Format