Context Navigation

physiq.F90 @ 1057

Last change on this file since 1057 was 317, checked in by millour, 9 years ago
Add possibility to output Ls (as a surface field). EM.
Property svn:executable set to ``*
File size: 96.4 KB

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

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

Context Navigation

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

Download in other formats: