Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/physiq.F90 @ 310

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

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

Download in other formats: