Context Navigation

physiq_mod.F90 @ 3882

Last change on this file since 3882 was 3882, checked in by debatzbr, 5 months ago
Pluto PCM: Clean outputs for photolysis and mufi BBT
Property svn:executable set to ``*
File size: 101.2 KB

Line
1	module physiq_mod
2
3	implicit none
4
5	contains
6
7	subroutine physiq(ngrid,nlayer,nq, &
8	firstcall,lastcall, &
9	pday,ptime,ptimestep, &
10	pplev,pplay,pphi, &
11	pu,pv,pt,pq, &
12	flxw, &
13	pdu,pdv,pdt,pdq,pdpsrf)
14
15	!!
16	use write_field_phy, only: Writefield_phy
17	!!
18	use ioipsl_getin_p_mod, only: getin_p
19	use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind, corrkdir, banddir
20	use gases_h, only: gnom, gfrac
21	use radcommon_h, only: sigma, glat, grav, BWNV, WNOI, DWNI, DWNV, WNOV
22	use suaer_corrk_mod, only: suaer_corrk
23	use aerosol_mod, only: i_haze, haze_prof
24	use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, &
25	dryness
26	use comdiurn_h, only: coslat, sinlat, coslon, sinlon
27	use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
28	use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
29	use geometry_mod, only: latitude, longitude, cell_area, &
30	cell_area_for_lonlat_outputs
31	USE comgeomfi_h, only: totarea, totarea_planet
32	USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
33	igcm_n2,igcm_ch4_gas,igcm_ch4_ice,igcm_haze,&
34	igcm_co_gas,igcm_co_ice,igcm_prec_haze,lw_n2,lw_ch4,lw_co,&
35	alpha_lift, alpha_devil, qextrhor, &
36	nesp, is_chim
37	use time_phylmdz_mod, only: diagfi_output_rate, startfi_output_rate, nday
38	use phyetat0_mod, only: phyetat0
39	use tabfi_mod, only: tab_cntrl_mod
40	use wstats_mod, only: callstats, wstats, mkstats
41	use phyredem, only: physdem0, physdem1
42	use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
43	use mod_phys_lmdz_para, only : is_master
44	use planete_mod, only: apoastr, periastr, year_day, peri_day, &
45	obliquit, z0, adjust, tpal
46	use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
47	use calldrag_noro_mod, only: calldrag_noro
48	use time_phylmdz_mod, only: daysec
49	use callkeys_mod, only: albedo_spectral_mode, calladj, calldifv, &
50	calllott, callrad, callsoil, nosurf, &
51	callconduct,callmolvis,callmoldiff, &
52	corrk, &
53	diurnal, enertest, fat1au, &
54	icetstep, intheat, iradia, kastprof, &
55	lwrite, mass_redistrib, meanOLR, &
56	fast,fasthaze,haze,metcloud,monoxcloud, &
57	n2cond,noseason_day,conservn2,conservch4, &
58	convergeps,kbo,triton,paleo,paleoyears,glaflow, &
59	carbox, methane,condmetsurf,condcosurf, &
60	oldplutovdifc,oldplutocorrk,oldplutosedim, &
61	optichaze,haze_proffix,haze_radproffix, &
62	source_haze, tsurfmax, albmin_ch4, &
63	season, sedimentation, &
64	specOLR, &
65	startphy_file, testradtimes, &
66	tracer, UseTurbDiff, &
67	global1d, szangle, &
68	callmufi
69	use check_fields_mod, only: check_physics_fields
70	use conc_mod, only: rnew, cpnew, ini_conc_mod
71	use phys_state_var_mod
72	use callcorrk_mod, only: callcorrk
73	use callcorrk_pluto_mod, only: callcorrk_pluto
74	use surfprop_mod, only: surfprop
75	use vdifc_mod, only: vdifc
76	use vdifc_pluto_mod, only: vdifc_pluto
77	use turbdiff_mod, only: turbdiff
78	use turb_mod, only : q2,sensibFlux,turb_resolved
79	use mass_redistribution_mod, only: mass_redistribution
80	use datafile_mod, only: datadir
81	USE vertical_layers_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt
82	use mod_phys_lmdz_omp_data, ONLY: is_omp_master
83	USE mod_grid_phy_lmdz, ONLY: regular_lonlat, grid_type, unstructured
84	! Microphysical model (mp2m)
85	use mp2m_calmufi
86	use mp2m_diagnostics
87
88	#ifdef CPP_XIOS
89	use xios_output_mod, only: initialize_xios_output, &
90	update_xios_timestep, &
91	send_xios_field
92	use wxios, only: wxios_context_init, xios_context_finalize
93	#endif
94	USE mod_grid_phy_lmdz, ONLY: grid_type,unstructured,regular_lonlat
95	use write_output_mod, only: write_output
96
97	implicit none
98
99
100	!==================================================================
101	!
102	! Purpose
103	! -------
104	! Central subroutine for all the physics parameterisations in the
105	! Pluto model. Originally adapted from the Mars LMDZ model.
106	!
107	! The model can be run with 1 (N2) or more tracer transport
108	! depending on the value of "tracer" in file "callphys.def".
109	!
110	! It includes:
111	!
112	! I. Initialization :
113	! I.1 Firstcall initializations.
114	! I.2 Initialization for every call to physiq.
115	!
116	! II.1 Thermosphere
117	! II.2 Compute radiative transfer tendencies (longwave and shortwave) :
118	! II.2.a Option 1 : Call correlated-k radiative transfer scheme.
119	! II.2.b Option 2 : Atmosphere has no radiative effect.
120	!
121	! II.3 Gravity wave and subgrid scale topography drag :
122	!
123	! III. Vertical diffusion (turbulent mixing)
124	!
125	! IV. Convection :
126	! IV.a Dry convective adjusment
127	!
128	! V. Condensation and sublimation of gases.
129	!
130	! VI. Tracers
131	! VI.1. Microphysics / Aerosols and particles.
132	! VI.2. Updates (pressure variations, surface budget).
133	! VI.3. Surface Tracer Update.
134	!
135	! VII. Surface and sub-surface soil temperature.
136	!
137	! VIII. Perform diagnostics and write output files.
138	!
139	!
140	! arguments
141	! ---------
142	!
143	! INPUT
144	! -----
145	!
146	! ngrid Size of the horizontal grid.
147	! nlayer Number of vertical layers.
148	! nq Number of advected fields.
149	!
150	! firstcall True at the first call.
151	! lastcall True at the last call.
152	!
153	! pday Number of days counted from the North. Spring equinoxe.
154	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
155	! ptimestep timestep (s).
156	!
157	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa).
158	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
159	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2.s-2).
160	!
161	! pu(ngrid,nlayer) u, zonal component of the wind (ms-1).
162	! pv(ngrid,nlayer) v, meridional component of the wind (ms-1).
163	!
164	! pt(ngrid,nlayer) Temperature (K).
165	!
166	! pq(ngrid,nlayer,nq) Advected fields.
167	!
168	! pudyn(ngrid,nlayer) \
169	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
170	! ptdyn(ngrid,nlayer) / corresponding variables.
171	! pqdyn(ngrid,nlayer,nq) /
172	! flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
173	!
174	! OUTPUT
175	! ------
176	!
177	! pdu(ngrid,nlayer) \
178	! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding
179	! pdt(ngrid,nlayer) / variables due to physical processes.
180	! pdq(ngrid,nlayer) /
181	! pdpsrf(ngrid) /
182	!
183	!
184	! Authors
185	! -------
186	! Frederic Hourdin 15/10/93
187	! Francois Forget 1994
188	! Christophe Hourdin 02/1997
189	! Subroutine completely rewritten by F. Forget (01/2000)
190	! Water ice clouds: Franck Montmessin (update 06/2003)
191	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
192	! New correlated-k radiative scheme: R. Wordsworth (2009)
193	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
194	! Improved water cycle: R. Wordsworth / B. Charnay (2010)
195	! To F90: R. Wordsworth (2010)
196	! New turbulent diffusion scheme: J. Leconte (2012)
197	! Loops converted to F90 matrix format: J. Leconte (2012)
198	! No more ngrid/nq, F90 commons and adaptation to parallel: A. Spiga (2012)
199	! Purge of the code : M. Turbet (2015)
200	! Photochemical core developped by F. Lefevre: B. Charnay (2017)
201	! Purge for Pluto model : A. Falco (2024)
202	! Adapting to Pluto : A. Falco, T. Bertrand (2024)
203	! Microphysical moment model: B. de Batz de Trenquelléon (2024)
204	!==================================================================
205
206
207	! 0. Declarations :
208	! ------------------
209
210	include "netcdf.inc"
211
212	! Arguments :
213	! -----------
214
215	! INPUTS:
216	! -------
217
218	integer,intent(in) :: ngrid ! Number of atmospheric columns.
219	integer,intent(in) :: nlayer ! Number of atmospheric layers.
220	integer,intent(in) :: nq ! Number of tracers.
221
222	logical,intent(in) :: firstcall ! Signals first call to physics.
223	logical,intent(in) :: lastcall ! Signals last call to physics.
224
225	real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0.
226	real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
227	real,intent(in) :: ptimestep ! Physics timestep (s).
228	real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
229	real,intent(in) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
230	real,intent(in) :: pphi(ngrid,nlayer) ! Geopotential at mid-layer (m2s-2).
231	real,intent(in) :: pu(ngrid,nlayer) ! Zonal wind component (m/s).
232	real,intent(in) :: pv(ngrid,nlayer) ! Meridional wind component (m/s).
233	real,intent(in) :: pt(ngrid,nlayer) ! Temperature (K).
234	real,intent(in) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
235	real,intent(in) :: flxw(ngrid,nlayer) ! Vertical mass flux (ks/s) at lower boundary of layer
236
237	! OUTPUTS:
238	! --------
239
240	! Physical tendencies :
241
242	real,intent(out) :: pdu(ngrid,nlayer) ! Zonal wind tendencies (m/s/s).
243	real,intent(out) :: pdv(ngrid,nlayer) ! Meridional wind tendencies (m/s/s).
244	real,intent(out) :: pdt(ngrid,nlayer) ! Temperature tendencies (K/s).
245	real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
246	real,intent(out) :: pdpsrf(ngrid) ! Surface pressure tendency (Pa/s).
247
248	! Local saved variables:
249	! ----------------------
250	integer,save :: day_ini ! Initial date of the run (sol since Ls=0).
251	integer,save :: icount ! Counter of calls to physiq during the run.
252	!$OMP THREADPRIVATE(day_ini,icount)
253
254	!Pluto specific
255	REAL,save :: acond,bcond
256	REAL,save :: tcond1p4Pa
257	DATA tcond1p4Pa/38/
258
259	! Local variables :
260	! -----------------
261	! Tendencies for the paleoclimate mode
262	REAL qsurfyear(ngrid,nq) ! kg.m-2 averaged mass of ice lost/gained in the last Pluto year of the run
263	REAL phisfinew(ngrid) ! geopotential of the bedrock (= phisfi-qsurf/1000*g)
264	REAL qsurfpal(ngrid,nq) ! qsurf after a paleoclimate step : for physdem1 and restartfi
265	REAL phisfipal(ngrid) ! geopotential after a paleoclimate step : for physdem1 and restartfi
266	REAL oblipal ! change of obliquity
267	REAL peri_daypal ! new periday
268	REAL eccpal ! change of eccentricity
269	REAL tpalnew ! change of time
270	REAL adjustnew ! change in N2 ice albedo
271	REAL pdaypal ! new pday = day_ini + step
272	REAL zdt_tot ! time range corresponding to the flux of qsurfyear
273	REAL massacc(nq) ! accumulated mass
274	REAL masslost(nq) ! accumulated mass
275
276	REAL globave ! globalaverage 2D ps
277	REAL globaveice(nq) ! globalaverage 2D ice
278	REAL globavenewice(nq) ! globalaverage 2D ice
279	INTEGER lecttsoil ! lecture of tsoil from proftsoil
280	REAL qsurf1(ngrid,nq) ! saving qsurf to calculate flux over long timescales kg.m-2
281	REAL flusurf(ngrid,nq) ! flux cond/sub kg.m-2.s-1
282	REAL flusurfold(ngrid,nq) ! old flux cond/sub kg.m-2.s-1
283	REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
284
285	REAL,SAVE :: ptime0 ! store the first time
286	REAL dstep
287	REAL,SAVE :: glastep=20 ! step in pluto day to spread glacier
288
289	! Aerosol (dust or ice) extinction optical depth at reference wavelength
290	! for the "naerkind" optically active aerosols:
291
292	real,save,allocatable :: dtau_aer(:,:,:) ! Aerosols
293	!$OMP THREADPRIVATE(dtau_aer)
294	real zh(ngrid,nlayer) ! Potential temperature (K).
295	real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
296	real omega(ngrid,nlayer) ! omega velocity (Pa/s, >0 when downward)
297
298	integer i,l,ig,ierr,iq,nw,isoil,iesp
299
300	real zls ! Solar longitude (radians).
301	real zlss ! Sub solar point longitude (radians).
302	real zday ! Date (time since Ls=0, calculated in sols).
303	REAL,save :: saveday ! saved date
304	REAL,save :: savedeclin ! saved declin
305	real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers.
306	real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries.
307
308	! TENDENCIES due to various processes :
309
310	! For Surface Temperature : (K/s)
311	real zdtsurf(ngrid) ! Cumulated tendencies.
312	real zdtsurfmr(ngrid) ! Mass_redistribution routine.
313	real zdtsurfc(ngrid) ! Condense_n2 routine.
314	real zdtsdif(ngrid) ! Turbdiff/vdifc routines.
315
316	! For Atmospheric Temperatures : (K/s)
317	real dtlscale(ngrid,nlayer) ! Largescale routine.
318	real zdtc(ngrid,nlayer) ! Condense_n2 routine.
319	real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
320	real zdtmr(ngrid,nlayer) ! Mass_redistribution routine.
321	real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
322	real zdtchim(ngrid,nlayer) ! Calchim routine.
323
324	! For Surface Tracers : (kg/m2/s)
325	real dqsurf(ngrid,nq) ! Cumulated tendencies.
326	!real zdqsurfc(ngrid) ! Condense_n2 routine.
327	REAL zdqsc(ngrid,nq) ! Condense_n2 routine.
328	real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines.
329	real zdqssed(ngrid,nq) ! Callsedim routine.
330	real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine.
331
332	! For Tracers : (kg/kg_of_air/s)
333	real zdqc(ngrid,nlayer,nq) ! Condense_n2 routine.
334	real zdqadj(ngrid,nlayer,nq) ! Convadj routine.
335	real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines.
336	real zdqevap(ngrid,nlayer) ! Turbdiff routine.
337	real zdqsed(ngrid,nlayer,nq) ! Callsedim routine.
338	real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine.
339	REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine
340	REAL,allocatable,save :: zdqschim(:,:) ! Calchim_asis routine
341	!$OMP THREADPRIVATE(zdqchim,zdqschim)
342
343	!! PLUTO variables
344	REAL zdqch4cloud(ngrid,nlayer,nq)
345	REAL zdqsch4cloud(ngrid,nq)
346	REAL zdtch4cloud(ngrid,nlayer)
347	REAL zdqcocloud(ngrid,nlayer,nq)
348	REAL zdqscocloud(ngrid,nq)
349	REAL zdtcocloud(ngrid,nlayer)
350	REAL rice_ch4(ngrid,nlayer) ! Methane ice geometric mean radius (m)
351	REAL rice_co(ngrid,nlayer) ! CO ice geometric mean radius (m)
352
353	REAL zdqsch4fast(ngrid) ! used only for fast model nogcm
354	REAL zdqch4fast(ngrid) ! used only for fast model nogcm
355	REAL zdqscofast(ngrid) ! used only for fast model nogcm
356	REAL zdqcofast(ngrid) ! used only for fast model nogcm
357	REAL zdqflow(ngrid,nq)
358
359	REAL zdtconduc(ngrid,nlayer) ! (K/s)
360	REAL zdumolvis(ngrid,nlayer)
361	REAL zdvmolvis(ngrid,nlayer)
362	real zdqmoldiff(ngrid,nlayer,nq)
363
364	! Haze relatated tendancies
365	REAL zdqhaze(ngrid,nlayer,nq)
366	REAL zdqprodhaze(ngrid,nq)
367	REAL zdqsprodhaze(ngrid)
368	REAL zdqhaze_col(ngrid)
369	REAL zdqphot_prec(ngrid,nlayer)
370	REAL zdqphot_ch4(ngrid,nlayer)
371	REAL zdqconv_prec(ngrid,nlayer)
372	REAL zdq_source(ngrid,nlayer,nq)
373	! Fast Haze relatated tendancies
374	REAL fluxbot(ngrid)
375	REAL gradflux(ngrid)
376	REAL fluxlym_sol_bot(ngrid) ! Solar flux Lyman alpha ph.m-2.s-1 reaching the surface
377	REAL fluxlym_ipm_bot(ngrid) ! IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 reaching the surface
378	REAL flym_sol(ngrid) ! Incident Solar flux Lyman alpha ph.m-2.s-1
379	REAL flym_ipm(ngrid) ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1
380	REAL zfluxuv ! Lyman alpha flux at 1AU
381
382	REAL array_zero1(ngrid)
383	REAL array_zero2(ngrid,nlayer)
384
385	! For Winds : (m/s/s)
386	real zdvadj(ngrid,nlayer), zduadj(ngrid,nlayer) ! Convadj routine.
387	real zdumr(ngrid,nlayer), zdvmr(ngrid,nlayer) ! Mass_redistribution routine.
388	real zdvdif(ngrid,nlayer), zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
389	real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
390	real zdhadj(ngrid,nlayer) ! Convadj routine.
391	REAL zdtgw(ngrid,nlayer) ! Gravity waves (K/s)
392	REAL zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer) ! Gravity waves (m.s-2)
393	REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer) ! condense_n2 routine.
394
395	! For Pressure and Mass :
396	real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
397	real zdmassmr_col(ngrid) ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
398	real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s).
399
400	! Local variables for MICROPHYSICS:
401	! ---------------------------------
402	real gzlat(ngrid,nlayer) ! Altitude-Latitude-dependent gravity (this should be stored elsewhere...).
403	real pdqmufi(ngrid,nlayer,nq) ! Microphysical tendency (X/kg_of_air/s).
404	real pdqmufi_prod(ngrid,nlayer,nq) ! Aerosols production tendency (kg/kg_of_air/s).
405	real int2ext(ngrid,nlayer) ! Intensive to extensive factor (kg_air/m3: X/kg_air --> X/m3).
406
407	! Local variables for LOCAL CALCULATIONS:
408	! ---------------------------------------
409	real zflubid(ngrid)
410	real zplanck(ngrid),zpopsk(ngrid,nlayer)
411	REAL zdum1(ngrid,nlayer)
412	REAL zdum2(ngrid,nlayer)
413	real ztim1,ztim2,ztim3, z1,z2
414	real ztime_restart
415	real zdh(ngrid,nlayer)
416	real gmplanet
417	real taux(ngrid),tauy(ngrid)
418
419	! local variables for DIAGNOSTICS : (diagfi & stat)
420	! -------------------------------------------------
421	real ps(ngrid) ! Surface Pressure.
422	real zt(ngrid,nlayer) ! Atmospheric Temperature.
423	real zu(ngrid,nlayer),zv(ngrid,nlayer) ! Zonal and Meridional Winds.
424	real zq(ngrid,nlayer,nq) ! Atmospheric Tracers.
425	real zdtadj(ngrid,nlayer) ! Convadj Diagnostic.
426	real zdtdyn(ngrid,nlayer) ! Dynamical Heating (K/s).
427	real zdudyn(ngrid,nlayer) ! Dynamical Zonal Wind tendency (m.s-2).
428
429	real reff(ngrid,nlayer) ! Effective dust radius (used if doubleq=T).
430	real vmr(ngrid,nlayer) ! volume mixing ratio
431	real time_phys
432
433	real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
434
435	real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
436
437	! Pluto adding variables
438	real vmr_ch4(ngrid) ! vmr ch4
439	real vmr_co(ngrid) ! vmr co
440	real rho(ngrid,nlayer) ! density
441	real zrho_ch4(ngrid,nlayer) ! density methane kg.m-3
442	real zrho_co(ngrid,nlayer) ! density CO kg.m-3
443	real zrho_haze(ngrid,nlayer) ! density haze kg.m-3
444	real zdqrho_photprec(ngrid,nlayer) !photolysis rate kg.m-3.s-1
445	real zq1temp_ch4(ngrid) !
446	real qsat_ch4(ngrid) !
447	real qsat_ch4_l1(ngrid) !
448	! CHARACTER(LEN=20) :: txt ! to temporarly store text for eddy tracers
449	real profmmr(ngrid,nlayer) ! fixed profile of haze if haze_proffix
450	real sensiblehf1(ngrid) ! sensible heat flux
451	real sensiblehf2(ngrid) ! sensible heat flux
452
453	! included by RW for H2O Manabe scheme
454	real rneb_man(ngrid,nlayer) ! H2O cloud fraction (moistadj).
455	real rneb_lsc(ngrid,nlayer) ! H2O cloud fraction (large scale).
456
457	! to test energy conservation (RW+JL)
458	real mass(ngrid,nlayer),massarea(ngrid,nlayer)
459	real dEtot, dEtots, AtmToSurf_TurbFlux
460	real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
461	!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
462
463	!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
464	real dtmoist_max,dtmoist_min
465	real dItot, dItot_tmp, dVtot, dVtot_tmp
466	real dWtot, dWtot_tmp, dWtots, dWtots_tmp
467
468	real nconsMAX, vdifcncons(ngrid), cadjncons(ngrid) ! Vdfic water conservation test. By RW
469
470	real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW
471
472	! Non-oro GW tendencies
473	REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
474	REAL d_t_hin(ngrid,nlayer)
475	! Diagnostics 2D of gw_nonoro
476	REAL zustrhi(ngrid), zvstrhi(ngrid)
477
478	real :: tsurf2(ngrid)
479	!! real :: flux_o(ngrid),flux_g(ngrid)
480	real :: flux_g(ngrid)
481	real :: flux_sens_lat(ngrid)
482	real :: qsurfint(ngrid,nq)
483
484	! local variables for skin depth check
485	real :: therm_inertia(ngrid,nsoilmx)
486	real :: inertia_min,inertia_max
487	real :: diurnal_skin ! diurnal skin depth (m)
488	real :: annual_skin ! anuual skin depth (m)
489
490	! when no startfi file is asked for init
491	real alpha,lay1 ! coefficients for building layers
492	integer iloop
493
494	! flags to trigger extra sanity checks
495	logical, save :: check_physics_inputs=.false.
496	logical, save :: check_physics_outputs=.false.
497	!$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)
498	LOGICAL write_restartfi
499
500	! Misc
501	character*2 :: str2
502	character(len=10) :: tmp1
503	character(len=10) :: tmp2
504	!==================================================================================================
505
506	! -----------------
507	! I. INITIALISATION
508	! -----------------
509
510	! --------------------------------
511	! I.1 First Call Initialisation.
512	! --------------------------------
513	if (firstcall) then
514	call getin_p("check_physics_inputs", check_physics_inputs)
515	call getin_p("check_physics_outputs", check_physics_outputs)
516
517	! Allocate saved arrays (except for 1D model, where this has already
518	! been done)
519	if (ngrid>1) call phys_state_var_init(nq)
520
521	! Variables set to 0
522	! ~~~~~~~~~~~~~~~~~~
523	dtrad(:,:) = 0.0
524	fluxrad(:) = 0.0
525	tau_col(:) = 0.0
526	zdtsw(:,:) = 0.0
527	zdtlw(:,:) = 0.0
528	zdqc(:,:,:)=0.
529	zdqsc(:,:)=0.
530
531	! Initialize tracer names, indexes and properties.
532	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
533	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
534	if (tracer) then
535	call initracer(ngrid,nq)
536	! if(photochem) then !AF24: removed
537	endif
538	! Initialize aerosol indexes.
539	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
540	! call iniaerosol
541	! allocate related local arrays
542	! (need be allocated instead of automatic because of "naerkind")
543	allocate(dtau_aer(ngrid,nlayer,naerkind))
544
545	#ifdef CPP_XIOS
546	! Initialize XIOS context
547	write(,) "physiq: call wxios_context_init"
548	CALL wxios_context_init
549	#endif
550
551	! Read 'startfi.nc' file.
552	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
553	call phyetat0(startphy_file, &
554	ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
555	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,inertiedat)
556
557	if (.not.startphy_file) then
558	! starting without startfi.nc and with callsoil
559	! is not yet possible as soildepth default is not defined
560	if (callsoil) then
561	! default mlayer distribution, following a power law:
562	! mlayer(k)=lay1alpha*(k-1/2)
563	lay1=2.e-4
564	alpha=2
565	do iloop=0,nsoilmx-1
566	mlayer(iloop)=lay1(alpha*(iloop-0.5))
567	enddo
568	lay1=sqrt(mlayer(0)*mlayer(1))
569	alpha=mlayer(1)/mlayer(0)
570	do iloop=1,nsoilmx
571	layer(iloop)=lay1(alpha*(iloop-1))
572	enddo
573	endif
574	! additionnal "academic" initialization of physics
575	if (is_master) write(,) "Physiq: initializing tsurf(:) to pt(:,1) !!"
576	tsurf(:)=pt(:,1)
577	if (is_master) write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
578	do isoil=1,nsoilmx
579	tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
580	enddo
581	if (is_master) write(,) "Physiq: initializing day_ini to pday !"
582	day_ini=pday
583	endif
584
585	if (pday.ne.day_ini) then
586	write(,) "ERROR in physiq.F90:"
587	write(,) "bad synchronization between physics and dynamics"
588	write(,) "dynamics day: ",pday
589	write(,) "physics day: ",day_ini
590	stop
591	endif
592
593	write (,) 'In physiq day_ini =', day_ini
594
595	! Initialize albedo calculation.
596	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
597	albedo(:,:)=0.0
598	albedo_bareground(:)=0.0
599	albedo_snow_SPECTV(:)=0.0
600	albedo_n2_ice_SPECTV(:)=0.0
601
602	ptime0=ptime
603	write (,) 'In physiq ptime0 =', ptime
604
605	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV)
606
607	! Initialize orbital calculation.
608	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
609	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
610
611	savedeclin=0.
612	saveday=pday
613	adjust=0. ! albedo adjustment for convergeps
614
615	! Initialize soil.
616	! ~~~~~~~~~~~~~~~~
617	if (callsoil) then
618	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
619	ptimestep,tsurf,tsoil,capcal,fluxgrd)
620	else ! else of 'callsoil'.
621	print*,'WARNING! Thermal conduction in the soil turned off'
622	capcal(:)=1.e6
623	fluxgrd(:)=intheat
624	print*,'Flux from ground = ',intheat,' W m^-2'
625	endif ! end of 'callsoil'.
626
627	icount=1
628
629	! Initialize variable for dynamical heating and zonal wind tendency diagnostic
630	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
631	ztprevious(:,:)=pt(:,:)
632	zuprevious(:,:)=pu(:,:)
633
634	if(meanOLR)then
635	call system('rm -f rad_bal.out') ! to record global radiative balance.
636	call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
637	call system('rm -f h2o_bal.out') ! to record global hydrological balance.
638	endif
639
640	if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
641	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
642	ptimestep,pday+nday,time_phys,cell_area, &
643	albedo_bareground,zmea,zstd,zsig,zgam,zthe)
644	endif
645
646	! Initialize correlated-k.
647	! ~~~~~~~~~~~~~~~~~~~~~~~~
648	if (corrk) then
649	! We initialise the spectral grid here instead of
650	! at firstcall of callcorrk so we can output XspecIR, XspecVI
651	! when using Dynamico
652	if (is_master) print*, "physiq_mod: Correlated-k data base folder:",trim(datadir)
653	call getin_p("corrkdir",corrkdir)
654	if (is_master) print*,"corrkdir = ", corrkdir
655	write (tmp1, '(i4)') L_NSPECTI
656	write (tmp2, '(i4)') L_NSPECTV
657	banddir=trim(trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2)))
658	banddir=trim(trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir)))
659	call setspi !Basic infrared properties.
660	call setspv ! Basic visible properties.
661	call sugas_corrk ! Set up gaseous absorption properties.
662	if (optichaze) then
663	call suaer_corrk ! Set up aerosol optical properties.
664	endif
665	endif
666
667	! Initialize microphysics.
668	! ~~~~~~~~~~~~~~~~~~~~~~~~
669	IF (callmufi) THEN
670	! Initialize microphysics arrays.
671	call inimufi(ptimestep)
672	ENDIF ! end callmufi
673
674	!! call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
675	! XIOS outputs
676	#ifdef CPP_XIOS
677
678	if (is_master) write(,) "physiq: call initialize_xios_output"
679	call initialize_xios_output(pday,ptime,ptimestep,daysec,year_day, &
680	presnivs,pseudoalt,mlayer,WNOI,WNOV)
681	#endif
682
683	!! call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
684
685	write(,) "physiq: end of firstcall"
686	endif ! end of 'firstcall'
687
688	!! call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
689	!! call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
690
691	if (check_physics_inputs) then
692	!check the validity of input fields coming from the dynamics
693	call check_physics_fields("begin physiq:", pt, pu, pv, pplev, pq)
694	endif
695
696	! call writediagfi(ngrid,"pre_physical_rnat"," "," ",2,rnat)
697	! call writediagfi(ngrid,"pre_physical_capcal"," "," ",2,capcal)
698
699	! ------------------------------------------------------
700	! I.2 Initializations done at every physical timestep:
701	! ------------------------------------------------------
702
703	#ifdef CPP_XIOS
704	! update XIOS time/calendar
705	call update_xios_timestep
706	#endif
707
708	! Initialize various variables
709	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
710
711	pdt(1:ngrid,1:nlayer) = 0.0
712	zdtsurf(1:ngrid) = 0.0
713	pdq(1:ngrid,1:nlayer,1:nq) = 0.0
714	dqsurf(1:ngrid,1:nq)= 0.0
715	pdu(1:ngrid,1:nlayer) = 0.0
716	pdv(1:ngrid,1:nlayer) = 0.0
717	pdpsrf(1:ngrid) = 0.0
718	zflubid(1:ngrid) = 0.0
719	flux_sens_lat(1:ngrid) = 0.0
720	taux(1:ngrid) = 0.0
721	tauy(1:ngrid) = 0.0
722
723	if (conservn2) then
724	write(,) 'conservn2 iniloop'
725	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
726	endif
727
728	zday=pday+ptime ! Compute time, in sols (and fraction thereof).
729
730	! Compute Stellar Longitude (Ls), and orbital parameters.
731	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
732	if (season) then
733	call stellarlong(zday,zls)
734	else
735	call stellarlong(noseason_day,zls)
736	end if
737
738	! Get Lyman alpha flux at specific Ls
739	if (callmufi.or.haze) then
740	call lymalpha(zls,zfluxuv)
741	if (is_master) then
742	print*, 'Haze lyman-alpha zls,zfluxuv=',zls,zfluxuv
743	endif
744	end if
745
746	IF (triton) then
747	CALL orbitetriton(zls,zday,dist_star,declin)
748	ELSE
749	call orbite(zls,dist_star,declin,right_ascen)
750	ENDIF
751
752	if (diurnal) then
753	zlss=-2.pi(zday-.5)
754	else if(diurnal .eqv. .false.) then
755	zlss=9999.
756	endif
757
758	glat(:) = g !AF24: removed oblateness
759
760	! Compute geopotential between layers.
761	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
762	zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
763	do l=1,nlayer
764	zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
765	enddo
766
767	zzlev(1:ngrid,1)=0.
768
769	do l=2,nlayer
770	do ig=1,ngrid
771	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
772	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
773	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
774	enddo
775	enddo
776
777	!Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
778
779	zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
780
781	! Compute potential temperature
782	! Note : Potential temperature calculation may not be the same in physiq and dynamic...
783	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
784	do l=1,nlayer
785	do ig=1,ngrid
786	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
787	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
788	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
789	massarea(ig,l)=mass(ig,l)*cell_area(ig)
790	enddo
791	enddo
792
793	! Compute vertical velocity (m/s) from vertical mass flux
794	! w = F / (rhoarea) and rho = P/(rT)
795	! But first linearly interpolate mass flux to mid-layers
796	if (.not.fast) then
797	do l=1,nlayer-1
798	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
799	enddo
800	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
801	do l=1,nlayer
802	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) / &
803	(pplay(1:ngrid,l)*cell_area(1:ngrid))
804	enddo
805	! omega in Pa/s
806	do l=1,nlayer-1
807	omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
808	enddo
809	omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
810	do l=1,nlayer
811	omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid)
812	enddo
813	endif
814
815	if (conservn2) then
816	write(,) 'conservn2 thermo'
817	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
818	endif
819
820	! Compute variations of g with latitude (to do).
821	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
822	gzlat(:,:) = g
823
824	! Initialize microphysical diagnostics.
825	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
826	IF (callmufi) THEN
827	! Initialize intensive to extensive factor (kg_air/m3: X/kg_air --> X/m3).
828	int2ext(:,:) = (pplev(:,1:nlayer)-pplev(:,2:nlayer+1)) / gzlat(:,1:nlayer) / (zzlev(:,2:nlayer+1)-zzlev(:,1:nlayer))
829
830	! Initialize microphysics diagnostics arrays.
831	call inimufi_diag(ngrid,nlayer,nq,pq,int2ext)
832	ENDIF ! end callmufi
833
834	! --------------------------------------------------------
835	! II.1 Thermosphere
836	! --------------------------------------------------------
837
838	! ajout de la conduction depuis la thermosphere
839	IF (callconduct) THEN
840
841	call conduction (ngrid,nlayer,ptimestep, &
842	pplay,pt,zzlay,zzlev,zdtconduc,tsurf)
843	DO l=1,nlayer
844	DO ig=1,ngrid
845	pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
846	ENDDO
847	ENDDO
848
849	ENDIF
850
851	! ajout de la viscosite moleculaire
852	IF (callmolvis) THEN
853	call molvis(ngrid,nlayer,ptimestep, &
854	pplay,pt,zzlay,zzlev, &
855	zdtconduc,pu,tsurf,zdumolvis)
856	call molvis(ngrid,nlayer,ptimestep, &
857	pplay,pt,zzlay,zzlev, &
858	zdtconduc,pv,tsurf,zdvmolvis)
859
860	DO l=1,nlayer
861	DO ig=1,ngrid
862	! pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
863	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
864	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
865	ENDDO
866	ENDDO
867	ENDIF
868
869	IF (callmoldiff) THEN
870	call moldiff_red(ngrid,nlayer,nq, &
871	pplay,pplev,pt,pdt,pq,pdq,ptimestep, &
872	zzlay,zdtconduc,zdqmoldiff)
873
874	DO l=1,nlayer
875	DO ig=1,ngrid
876	DO iq=1, nq
877	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
878	ENDDO
879	ENDDO
880	ENDDO
881	ENDIF
882
883	if (conservn2) then
884	write(,) 'conservn2 thermosphere'
885	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
886	endif
887
888
889	!---------------------------------
890	! II.2 Compute radiative tendencies
891	!---------------------------------
892	! Saving qsurf to compute paleo flux condensation/sublimation
893	DO iq=1, nq
894	DO ig=1,ngrid
895	IF (qsurf(ig,iq).lt.0.) then
896	qsurf(ig,iq)=0.
897	ENDIF
898	qsurf1(ig,iq)=qsurf(ig,iq)
899	ENDDO
900	ENDDO
901
902
903	! Compute local stellar zenith angles
904	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
905	fract = 0
906	if (diurnal) then
907	ztim1=SIN(declin)
908	ztim2=COS(declin)COS(2.pi*(zday-.5))
909	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
910
911	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
912	ztim1,ztim2,ztim3,mu0,fract)
913	else if(diurnal .eqv. .false.) then
914
915	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad)
916	! WARNING: this function appears not to work in 1D
917
918	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
919	mu0 = cos(pi*szangle/180.0)
920	fract= 1/(4*mu0) ! AF24: from pluto.old
921	endif
922
923	endif
924
925
926	! Pluto albedo /IT changes depending on surface ices (only in 3D)
927	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
928	if (ngrid.ne.1) then
929
930	!! Specific to change albedo of N2 so that Psurf
931	!! converges toward 1.4 Pa in "1989" seasons for Triton
932	!! converges toward 1.1 Pa in "2015" seasons for Pluto
933	if (convergeps) then
934	if (triton) then
935	! 1989 declination
936	if (declin180./pi.gt.-46..and.declin180./pi.lt.-45. &
937	.and.zday.gt.saveday+1000. &
938	.and.declin.lt.savedeclin) then
939	call planetwide_sumval(pplev(:,1)*cell_area(:)/totarea_planet,globave)
940	if (globave.gt.1.5) then
941	adjust=adjust+0.005
942	else if (globave.lt.1.3) then
943	adjust=adjust-0.005
944	endif
945	saveday=zday
946	endif
947	else
948	! Pluto : 2015 declination current epoch
949	if (declin180./pi.gt.50.and.declin180./pi.lt.51. &
950	.and.zday.gt.saveday+10000. &
951	.and.declin.gt.savedeclin) then
952	call planetwide_sumval(pplev(:,1)*cell_area(:)/totarea_planet,globave)
953	if (globave.gt.1.2) then
954	adjust=adjust+0.005
955	else if (globave.lt.1.) then
956	adjust=adjust-0.005
957	endif
958	saveday=zday
959	endif
960	endif
961	end if
962	end if ! if ngrid ne 1
963
964	call surfprop(ngrid,nq,fract,qsurf,tsurf, &
965	capcal,adjust,dist_star,flusurfold,ptimestep,zls,&
966	albedo,emis,therm_inertia)
967	! do i=2,ngrid
968	! albedo(i,:) = albedo(1,:)
969	! enddo
970	! AF24: TODO check albedo has been initialized here
971
972	if (firstcall.and.callsoil) then
973	! AF24 Originally in soil.F, but therm_inertia is modified by surfprop
974	! Additional checks: is the vertical discretization sufficient
975	! to resolve diurnal and annual waves?
976	do ig=1,ngrid
977	! extreme inertia for this column
978	inertia_min=minval(therm_inertia(ig,:))
979	inertia_max=maxval(therm_inertia(ig,:))
980	! diurnal and annual skin depth
981	diurnal_skin=(inertia_min/volcapa)*sqrt(daysec/pi)
982	annual_skin=(inertia_max/volcapa)sqrt(year_daydaysec/pi)
983	if (0.5*diurnal_skin<layer(1)) then
984	! one should have the fist layer be at least half of diurnal skin depth
985	write(,) "soil Error: grid point ig=",ig
986	write(,) " longitude=",longitude(ig)*(180./pi)
987	write(,) " latitude=",latitude(ig)*(180./pi)
988	write(,) " first soil layer depth ",layer(1)
989	write(,) " not small enough for a diurnal skin depth of ", &
990	diurnal_skin
991	write(,) " change soil layer distribution (comsoil_h.F90)"
992	call abort_physic("soil","change soil layer distribution (comsoil_h.F90)",1)
993	endif
994	if (2.*annual_skin>layer(nsoilmx)) then
995	! one should have the full soil be at least twice the diurnal skin depth
996	write(,) "soil Error: grid point ig=",ig
997	write(,) " longitude=",longitude(ig)*(180./pi)
998	write(,) " latitude=",latitude(ig)*(180./pi)
999	write(,) " total soil layer depth ",layer(nsoilmx)
1000	write(,) " not large enough for an annual skin depth of ", &
1001	annual_skin
1002	write(,) " change soil layer distribution (comsoil_h.F90)"
1003	call abort_physic("soil","change soil layer distribution (comsoil_h.F90)",1)
1004	endif
1005	enddo ! of do ig=1,ngrid
1006
1007	end if ! callsoil
1008
1009	if (callrad) then
1010	if( mod(icount-1,iradia).eq.0.or.lastcall) then
1011
1012	! Eclipse incoming sunlight !AF24: removed
1013
1014	!! call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1015	!! call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1016
1017
1018	if (corrk) then
1019
1020	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1021	! II.a Call correlated-k radiative transfer scheme
1022	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1023	if(kastprof)then
1024	print*,'kastprof should not = true here'
1025	call abort
1026	endif
1027
1028	! standard callcorrk
1029	if (oldplutocorrk) then
1030	call callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf, &
1031	albedo(:,1),emis,mu0,pplev,pplay,pt, &
1032	zzlay,zzlev,tsurf,fract,dist_star,dtau_aer, &
1033	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, &
1034	fluxtop_sw,fluxtop_dn,reffrad, &
1035	int_dtaui,int_dtauv,tau_col, &
1036	ptime,pday,firstcall,lastcall)
1037	albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
1038	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+ &
1039	fluxsurf_sw(1:ngrid)*(1.-albedo(1:ngrid,1))
1040	fluxabs_sw(1:ngrid)=fluxtop_dn(1:ngrid)-fluxtop_sw(1:ngrid)
1041	else
1042	muvar(1:ngrid,1:nlayer+1)=mugaz
1043	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
1044	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
1045	zzlay,zzlev,tsurf,fract,dist_star,dtau_aer,muvar, &
1046	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, &
1047	fluxsurfabs_sw,fluxtop_lw, &
1048	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,GSR_nu, &
1049	int_dtaui,int_dtauv, &
1050	tau_col,firstcall,lastcall)
1051	! Radiative flux from the sky absorbed by the surface (W.m-2).
1052	GSR=0.0
1053	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+ &
1054	fluxsurfabs_sw(1:ngrid)
1055	endif ! oldplutocorrk
1056	!GG (feb2021): Option to "artificially" decrease the raditive time scale in
1057	!the deep atmosphere press > 0.1 bar. Suggested by MT.
1058	!! COEFF_RAD to be "tuned" to facilitate convergence of tendency
1059
1060	!coeff_rad=0. ! 0 values, it doesn't accelerate the convergence
1061	!coeff_rad=0.5
1062	!coeff_rad=1.
1063	!do l=1, nlayer
1064	! do ig=1,ngrid
1065	! if(pplay(ig,l).ge.1.d4) then
1066	! zdtsw(ig,l)=zdtsw(ig,l)(pplay(ig,l)/1.d4)*coeff_rad
1067	! zdtlw(ig,l)=zdtlw(ig,l)(pplay(ig,l)/1.d4)*coeff_rad
1068	! endif
1069	! enddo
1070	!enddo
1071
1072	! AF24: removed CLFvarying Option
1073
1074
1075	!if(noradsurf)then ! no lower surface; SW flux just disappears
1076	! GSR = SUM(fluxsurf_sw(1:ngrid)*cell_area(1:ngrid))/totarea
1077	! fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
1078	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
1079	!endif
1080
1081	! Net atmospheric radiative heating rate (K.s-1)
1082	dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
1083
1084	! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that !
1085	if (firstcall .and. albedo_spectral_mode) then
1086	call spectral_albedo_calc(albedo_snow_SPECTV)
1087	endif
1088
1089	else
1090	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1091	! II.b Atmosphere has no radiative effect
1092	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1093	fluxtop_dn(1:ngrid) = fract(1:ngrid)mu0(1:ngrid)Fat1AU/dist_star**2
1094	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
1095	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
1096	endif
1097	fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
1098	print*,'------------WARNING---WARNING------------' ! by MT2015.
1099	print*,'You are in corrk=false mode, '
1100	print*,'and the surface albedo is taken equal to the first visible spectral value'
1101
1102	albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
1103	fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
1104	fluxabs_sw(1:ngrid)=fluxsurfabs_sw(1:ngrid)
1105	fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid)
1106	fluxtop_lw(1:ngrid) = emis(1:ngrid)sigmatsurf(1:ngrid)**4
1107
1108	dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating
1109
1110	endif ! end of corrk
1111
1112	endif ! of if(mod(icount-1,iradia).eq.0)
1113
1114
1115	! Transformation of the radiative tendencies
1116	! ------------------------------------------
1117	zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
1118	zplanck(1:ngrid)=emis(1:ngrid)sigmazplanck(1:ngrid)*zplanck(1:ngrid)
1119	fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
1120	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
1121
1122	! Test of energy conservation
1123	!----------------------------
1124	if(enertest)then
1125	call planetwide_sumval(cppmassarea(:,:)zdtsw(:,:)/totarea_planet,dEtotSW)
1126	call planetwide_sumval(cppmassarea(:,:)zdtlw(:,:)/totarea_planet,dEtotLW)
1127	!call planetwide_sumval(fluxsurf_sw(:)(1.-albedo_equivalent(:))cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
1128	call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
1129	call planetwide_sumval((fluxsurf_lw(:)emis(:)-zplanck(:))cell_area(:)/totarea_planet,dEtotsLW)
1130	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
1131	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
1132	if (is_master) then
1133	print*,'---------------------------------------------------------------'
1134	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
1135	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
1136	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
1137	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
1138	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
1139	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
1140	endif
1141	endif ! end of 'enertest'
1142
1143	endif ! of if (callrad)
1144
1145	!! call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1146	!! call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1147
1148	if (conservn2) then
1149	write(,) 'conservn2 radiat'
1150	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
1151	endif
1152
1153
1154	!-----------------------------------------------------------------------
1155	! II.3 Gravity wave and subgrid scale topography drag :
1156	! -------------------------------------------------
1157	IF(calllott)THEN
1158	CALL calldrag_noro(ngrid,nlayer,ptimestep, &
1159	pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
1160
1161	DO l=1,nlayer
1162	DO ig=1,ngrid
1163	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
1164	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
1165	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
1166	ENDDO
1167	ENDDO
1168	ENDIF
1169
1170
1171	! --------------------------------------------
1172	! III. Vertical diffusion (turbulent mixing) :
1173	! --------------------------------------------
1174
1175	if (calldifv) then
1176
1177	zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
1178
1179	if (oldplutovdifc) then
1180	zdum1(:,:) = 0.
1181	zdum2(:,:) = 0.
1182	zdh(:,:)=pdt(:,:)/zpopsk(:,:)
1183
1184	! Calling vdif (Martian version WITH N2 condensation)
1185	CALL vdifc_pluto(ngrid,nlayer,nq,zpopsk, &
1186	ptimestep,capcal,lwrite, &
1187	pplay,pplev,zzlay,zzlev,z0, &
1188	pu,pv,zh,pq,pt,tsurf,emis,qsurf, &
1189	zdum1,zdum2,zdh,pdq,pdt,zflubid, &
1190	zdudif,zdvdif,zdhdif,zdtsdif,q2, &
1191	zdqdif,zdqsdif,qsat_ch4,qsat_ch4_l1) !,zq1temp_ch4,qsat_ch4)
1192
1193	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
1194
1195	bcond=1./tcond1p4Pa
1196	acond=r/lw_n2
1197
1198	! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
1199	else if (UseTurbDiff) then
1200
1201	call turbdiff(ngrid,nlayer,nq, &
1202	ptimestep,capcal, &
1203	pplay,pplev,zzlay,zzlev,z0, &
1204	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
1205	pdt,pdq,zflubid, &
1206	zdudif,zdvdif,zdtdif,zdtsdif, &
1207	sensibFlux,q2,zdqdif,zdqevap,zdqsdif, &
1208	taux,tauy)
1209
1210	else ! if .not. (oldplutovdifc) .and. (UseTurbDiff)
1211
1212	zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1213
1214	call vdifc(ngrid,nlayer,nq,zpopsk, &
1215	ptimestep,capcal,lwrite, &
1216	pplay,pplev,zzlay,zzlev,z0, &
1217	pu,pv,zh,pq,tsurf,emis,qsurf, &
1218	zdh,pdq,zflubid, &
1219	zdudif,zdvdif,zdhdif,zdtsdif, &
1220	sensibFlux,q2,zdqdif,zdqsdif)
1221
1222	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1223	zdqevap(1:ngrid,1:nlayer)=0.
1224
1225	end if !end of 'UseTurbDiff'
1226
1227	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
1228
1229	!!! this is always done, except for turbulence-resolving simulations
1230	if (.not. turb_resolved) then
1231	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
1232	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
1233	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
1234	endif
1235
1236	if (tracer) then
1237	pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
1238	dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
1239	end if ! of if (tracer)
1240
1241	! test energy conservation
1242	!-------------------------
1243	if(enertest)then
1244
1245	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
1246	do ig = 1, ngrid
1247	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
1248	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
1249	enddo
1250
1251	call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
1252	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
1253	call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
1254	call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
1255
1256	if (is_master) then
1257
1258	if (UseTurbDiff) then
1259	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1260	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
1261	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1262	else
1263	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1264	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
1265	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1266	end if
1267	endif ! end of 'is_master'
1268
1269	! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
1270	endif ! end of 'enertest'
1271
1272	else ! calldifv
1273
1274	ztim1=4.sigmaptimestep
1275	DO ig=1,ngrid
1276	ztim2=ztim1emis(ig)tsurf(ig)**3
1277	z1=capcal(ig)*tsurf(ig)+ &
1278	ztim2tsurf(ig)+ (fluxrad(ig)+fluxgrd(ig))ptimestep
1279	z2= capcal(ig)+ztim2
1280	zdtsurf(ig)=(z1/z2 - tsurf(ig))/ptimestep
1281
1282	! for output:
1283	!dplanck(ig)=4.stephanptimestepemis(ig)tsurf(ig)**3
1284	ENDDO
1285
1286	! if(.not.newtonian)then
1287	!zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
1288
1289	! ------------------------------------------------------------------
1290	! Methane surface sublimation and condensation in fast model (nogcm)
1291	! ------------------------------------------------------------------
1292	if ((methane).and.(fast).and.condmetsurf) THEN
1293
1294	call ch4surf(ngrid,nlayer,nq,ptimestep,capcal, &
1295	tsurf,zdtsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf, &
1296	zdqch4fast,zdqsch4fast)
1297
1298	dqsurf(1:ngrid,igcm_ch4_ice)= dqsurf(1:ngrid,igcm_ch4_ice) + &
1299	zdqsch4fast(1:ngrid)
1300	pdq(1:ngrid,1,igcm_ch4_gas)= pdq(1:ngrid,1,igcm_ch4_gas) + &
1301	zdqch4fast(1:ngrid)
1302	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+lw_ch4*zdqsch4fast(1:ngrid)/capcal(1:ngrid)
1303	end if
1304	! ------------------------------------------------------------------
1305	! CO surface sublimation and condensation in fast model (nogcm)
1306	! ------------------------------------------------------------------
1307	if ((carbox).and.(fast).and.condcosurf) THEN
1308
1309	call cosurf(ngrid,nlayer,nq,ptimestep, &
1310	tsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf, &
1311	zdqcofast,zdqscofast)
1312
1313	dqsurf(1:ngrid,igcm_co_ice)= dqsurf(1:ngrid,igcm_co_ice) + &
1314	zdqscofast(1:ngrid)
1315	pdq(1:ngrid,1,igcm_co_gas)= pdq(1:ngrid,1,igcm_co_gas) + &
1316	zdqcofast(1:ngrid)
1317	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+lw_co*zdqscofast(1:ngrid)/capcal(1:ngrid)
1318	end if
1319
1320
1321	endif ! end of 'calldifv'
1322
1323	if (conservn2) then
1324	write(,) 'conservn2 calldifv'
1325	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)+ &
1326	dqsurf(:,1)*ptimestep)
1327	endif
1328	if (methane.and.conservch4) then
1329	write(,) 'conservch4 calldifv'
1330	if (fast) then
1331	call testconservfast(ngrid,nlayer,nq,pq(:,1,igcm_ch4_gas),pdq(:,1,igcm_ch4_gas), &
1332	qsurf(:,igcm_ch4_ice),dqsurf(:,igcm_ch4_ice), &
1333	ptimestep,pplev,zdqch4fast,zdqsch4fast,'CH4',' vdifc ')
1334	else
1335	call testconserv(ngrid,nlayer,nq,pq,pdq,qsurf,dqsurf, &
1336	igcm_ch4_gas,igcm_ch4_ice, &
1337	ptimestep,pplev,zdqdif,zdqsdif,'CH4',' vdifc ')
1338	endif
1339	endif
1340
1341	!-------------------
1342	! IV. Convection :
1343	!-------------------
1344
1345	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1346	! IV.a Dry convective adjustment :
1347	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1348
1349	if(calladj) then
1350
1351	zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1352	zduadj(1:ngrid,1:nlayer)=0.0
1353	zdvadj(1:ngrid,1:nlayer)=0.0
1354	zdhadj(1:ngrid,1:nlayer)=0.0
1355
1356
1357	call convadj(ngrid,nlayer,nq,ptimestep, &
1358	pplay,pplev,zpopsk, &
1359	pu,pv,zh,pq, &
1360	pdu,pdv,zdh,pdq, &
1361	zduadj,zdvadj,zdhadj, &
1362	zdqadj)
1363
1364	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
1365	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer)
1366	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
1367	zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1368
1369	if(tracer) then
1370	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
1371	end if
1372
1373	! Test energy conservation
1374	if(enertest)then
1375	call planetwide_sumval(cppmassarea(:,:)zdtadj(:,:)/totarea_planet,dEtot)
1376	if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
1377	endif
1378
1379	! ! Test water conservation !AF24: removed
1380
1381	endif ! end of 'calladj'
1382
1383	!-----------------------------------------------
1384	! V. Nitrogen condensation-sublimation :
1385	!-----------------------------------------------
1386
1387	if (n2cond) then
1388
1389	if (.not.tracer) then
1390	print*,'We need a N2 ice tracer to condense N2'
1391	call abort
1392	endif
1393
1394	call condense_n2(ngrid,nlayer,nq,ptimestep, &
1395	capcal,pplay,pplev,tsurf,pt, &
1396	pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, &
1397	qsurf(1,igcm_n2),albedo,emis, &
1398	zdtc,zdtsurfc,pdpsrf,zduc,zdvc, &
1399	zdqc,zdqsc(1,igcm_n2))
1400
1401	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
1402	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer)+zdvc(1:ngrid,1:nlayer)
1403	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer)+zduc(1:ngrid,1:nlayer)
1404	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
1405
1406	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq)
1407	dqsurf(1:ngrid,igcm_n2) = dqsurf(1:ngrid,igcm_n2) + zdqsc(1:ngrid,igcm_n2)
1408
1409	!! call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1410	!! call writediagfi(ngrid,"condense_n2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1411
1412	! test energy conservation
1413	if(enertest)then
1414	call planetwide_sumval(cppmassarea(:,:)zdtc(:,:)/totarea_planet,dEtot)
1415	call planetwide_sumval(capcal(:)zdtsurfc(:)cell_area(:)/totarea_planet,dEtots)
1416	if (is_master) then
1417	print*,'In n2cloud atmospheric energy change =',dEtot,' W m-2'
1418	print*,'In n2cloud surface energy change =',dEtots,' W m-2'
1419	endif
1420	endif
1421
1422	endif ! end of 'n2cond'
1423
1424	if (conservn2) then
1425	write(,) 'conservn2 n2cond'
1426	call testconservmass(ngrid,nlayer,pplev(:,1)+ &
1427	pdpsrf(:)ptimestep,qsurf(:,1)+dqsurf(:,1)ptimestep)
1428	endif
1429	if (methane.and.conservch4) then
1430	write(,) 'conservch4 n2cond'
1431	if (fast) then
1432	call testconservfast(ngrid,nlayer,nq,pq(:,1,igcm_ch4_gas),pdq(:,1,igcm_ch4_gas), &
1433	qsurf(:,igcm_ch4_ice),dqsurf(:,igcm_ch4_ice), &
1434	ptimestep,pplev,zdqch4fast,zdqsch4fast,'CH4',' n2cond')
1435	else
1436	call testconserv(ngrid,nlayer,nq,pq,pdq,qsurf,dqsurf, &
1437	igcm_ch4_gas,igcm_ch4_ice, &
1438	ptimestep,pplev,zdqc,zdqsc,'CH4',' n2cond')
1439	endif
1440	endif
1441
1442	!---------------------------------------------
1443	! VI. Specific parameterizations for tracers
1444	!---------------------------------------------
1445
1446	if (tracer) then
1447
1448	! ---------------------------------------
1449	! Methane ice condensation in the atmosphere
1450	! ----------------------------------------
1451	rice_ch4(:,:)=0 ! initialization needed for callsedim
1452	zdqch4cloud(:,:,:)=0.
1453	if ((methane).and.(metcloud).and.(.not.fast)) THEN
1454	call ch4cloud(ngrid,nlayer,ptimestep, &
1455	pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, &
1456	pq,pdq,zdqch4cloud,zdqsch4cloud,zdtch4cloud, &
1457	nq,rice_ch4)
1458
1459	DO l=1,nlayer
1460	DO ig=1,ngrid
1461	pdq(ig,l,igcm_ch4_gas)=pdq(ig,l,igcm_ch4_gas)+ &
1462	zdqch4cloud(ig,l,igcm_ch4_gas)
1463	pdq(ig,l,igcm_ch4_ice)=pdq(ig,l,igcm_ch4_ice)+ &
1464	zdqch4cloud(ig,l,igcm_ch4_ice)
1465	ENDDO
1466	ENDDO
1467
1468	! Increment methane ice surface tracer tendency
1469	DO ig=1,ngrid
1470	dqsurf(ig,igcm_ch4_ice)=dqsurf(ig,igcm_ch4_ice)+ &
1471	zdqsch4cloud(ig,igcm_ch4_ice)
1472	ENDDO
1473
1474	! update temperature tendancy
1475	DO ig=1,ngrid
1476	DO l=1,nlayer
1477	pdt(ig,l)=pdt(ig,l)+zdtch4cloud(ig,l)
1478	ENDDO
1479	ENDDO
1480	end if
1481
1482	! ---------------------------------------
1483	! CO ice condensation in the atmosphere
1484	! ----------------------------------------
1485	zdqcocloud(:,:,:)=0.
1486	IF ((carbox).and.(monoxcloud).and.(.not.fast)) THEN
1487	call cocloud(ngrid,nlayer,ptimestep, &
1488	pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, &
1489	pq,pdq,zdqcocloud,zdqscocloud,zdtcocloud, &
1490	nq,rice_co,qsurf(1,igcm_n2),dqsurf(1,igcm_n2))
1491
1492	DO l=1,nlayer
1493	DO ig=1,ngrid
1494	pdq(ig,l,igcm_co_gas)=pdq(ig,l,igcm_co_gas)+ &
1495	zdqcocloud(ig,l,igcm_co_gas)
1496	pdq(ig,l,igcm_co_ice)=pdq(ig,l,igcm_co_ice)+ &
1497	zdqcocloud(ig,l,igcm_co_ice)
1498	ENDDO
1499	ENDDO
1500
1501	! Increment CO ice surface tracer tendency
1502	DO ig=1,ngrid
1503	dqsurf(ig,igcm_co_ice)=dqsurf(ig,igcm_co_ice)+ &
1504	zdqscocloud(ig,igcm_co_ice)
1505	ENDDO
1506
1507	! update temperature tendancy
1508	DO ig=1,ngrid
1509	DO l=1,nlayer
1510	pdt(ig,l)=pdt(ig,l)+zdtcocloud(ig,l)
1511	ENDDO
1512	ENDDO
1513	ELSE
1514	rice_co(:,:)=0 ! initialization needed for callsedim
1515	END IF ! of IF (carbox)
1516
1517	! ----------------------------------------
1518	! VI.1. Microphysics / Aerosol particles
1519	! ----------------------------------------
1520	! Call of microphysics
1521	IF (callmufi) THEN
1522
1523	! Production for microphysics
1524	IF (call_haze_prod_pCH4) THEN
1525	zdqphot_prec(:,:) = 0.
1526	zdqphot_ch4(:,:) = 0.
1527	pdqmufi_prod(:,:,:) = 0.
1528	call hazecloud(ngrid,nlayer,nq,ptimestep,zday, &
1529	pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,pdqmufi_prod, &
1530	zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
1531	pdq(:,:,:) = pdq(:,:,:) + pdqmufi_prod(:,:,:) ! Should be updated
1532	ENDIF ! end call_haze_prod_pCH4
1533
1534	pdqmufi(:,:,:) = 0.
1535
1536	call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,pdqmufi_prod,pdqmufi)
1537
1538	pdq(:,:,:) = pdq(:,:,:) + pdqmufi(:,:,:)
1539
1540	ELSE
1541	IF (haze) THEN
1542	zdqphot_prec(:,:) = 0.
1543	zdqphot_ch4(:,:) = 0.
1544	zdqhaze(:,:,:) = 0.
1545	! Forcing to a fixed haze profile if haze_proffix
1546	if (haze_proffix.and.i_haze.gt.0.) then
1547	call haze_prof(ngrid,nlayer,zzlay,pplay,pt, &
1548	reffrad,profmmr)
1549	zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze))/ptimestep
1550	else
1551	call hazecloud(ngrid,nlayer,nq,ptimestep,zday, &
1552	pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,zdqhaze, &
1553	zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
1554	endif
1555	pdq(:,:,:) = pdq(:,:,:) + zdqhaze(:,:,:) ! Should be updated
1556	ENDIF ! end haze
1557
1558	IF (fast.and.fasthaze) THEN
1559	call prodhaze(ngrid,nlayer,nq,ptimestep,pplev,pq,pdq,dist_star, &
1560	mu0,declin,zdqprodhaze,zdqsprodhaze,gradflux,fluxbot, &
1561	fluxlym_sol_bot,fluxlym_ipm_bot,flym_sol,flym_ipm)
1562	DO ig=1,ngrid
1563	pdq(ig,1,igcm_ch4_gas)=pdq(ig,1,igcm_ch4_gas)+ &
1564	zdqprodhaze(ig,igcm_ch4_gas)
1565	pdq(ig,1,igcm_prec_haze)=pdq(ig,1,igcm_prec_haze)+ &
1566	zdqprodhaze(ig,igcm_prec_haze)
1567	pdq(ig,1,igcm_haze)=abs(pdq(ig,1,igcm_haze)+ &
1568	zdqprodhaze(ig,igcm_haze))
1569	qsurf(ig,igcm_haze)= qsurf(ig,igcm_haze)+ &
1570	zdqsprodhaze(ig)*ptimestep
1571	ENDDO
1572	ENDIF ! end fast.and.fasthaze
1573
1574	! Sedimentation.
1575	if (sedimentation) then
1576	zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
1577	zdqssed(1:ngrid,1:nq) = 0.0
1578	if (oldplutosedim)then
1579	call callsedim_pluto(ngrid,nlayer,ptimestep, &
1580	pplev,zzlev,pt,pdt,rice_ch4,rice_co, &
1581	pq,pdq,zdqsed,zdqssed,nq,pphi)
1582	else
1583	call callsedim(ngrid,nlayer,ptimestep, &
1584	pplev,zzlev,pt,pdt,pq,pdq, &
1585	zdqsed,zdqssed,nq,pphi)
1586	endif
1587	! Whether it falls as rain or snow depends only on the surface temperature
1588	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
1589	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
1590	end if ! end of 'sedimentation'
1591
1592	ENDIF ! end callmufi
1593
1594	! ---------------
1595	! VI.2. Updates
1596	! ---------------
1597
1598	! Updating Atmospheric Mass and Tracers budgets.
1599	if(mass_redistrib) then
1600
1601	zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * 0
1602	! ( zdqevap(1:ngrid,1:nlayer) &
1603	! ! + zdqrain(1:ngrid,1:nlayer,igcm_h2o_gas) &
1604	! ! + dqmoist(1:ngrid,1:nlayer,igcm_h2o_gas) &
1605	! + dqvaplscale(1:ngrid,1:nlayer) )
1606
1607	do ig = 1, ngrid
1608	zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
1609	enddo
1610
1611	! call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
1612	! call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
1613	call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
1614
1615	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
1616	capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
1617	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
1618	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
1619
1620	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
1621	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
1622	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
1623	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
1624	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
1625	pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
1626	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
1627
1628	endif
1629
1630	! call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1631
1632	! -----------------------------
1633	! VI.3. Surface Tracer Update
1634	! -----------------------------
1635
1636	qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
1637
1638	endif! end of if 'tracer'
1639
1640	if (conservn2) then
1641	write(,) 'conservn2 tracer'
1642	call testconservmass(ngrid,nlayer,pplev(:,1)+ &
1643	pdpsrf(:)*ptimestep,qsurf(:,1))
1644	endif
1645
1646	DO ig=1,ngrid
1647	flusurf(ig,igcm_n2)=(qsurf(ig,igcm_n2)- &
1648	qsurf1(ig,igcm_n2))/ptimestep
1649	flusurfold(ig,igcm_n2)=flusurf(ig,igcm_n2)
1650	if (methane) then
1651	flusurf(ig,igcm_ch4_ice)=(qsurf(ig,igcm_ch4_ice)- &
1652	qsurf1(ig,igcm_ch4_ice))/ptimestep
1653	flusurfold(ig,igcm_ch4_ice)=flusurf(ig,igcm_ch4_ice)
1654	endif
1655	if (carbox) then
1656	flusurf(ig,igcm_co_ice)=(qsurf(ig,igcm_co_ice)- &
1657	qsurf1(ig,igcm_co_ice))/ptimestep
1658	!flusurfold(ig,igcm_co_ice)=flusurf(ig,igcm_co_ice)
1659	endif
1660	ENDDO
1661
1662	!! Special source of haze particle !
1663	! todo: should be placed in haze and use tendency of n2 instead of flusurf
1664	IF (source_haze) THEN
1665	write(,) "Source haze not supported yet."
1666	stop
1667	! call hazesource(ngrid,nlayer,nq,ptimestep, &
1668	! pplev,flusurf,mu0,zdq_source)
1669
1670	DO iq=1, nq
1671	DO l=1,nlayer
1672	DO ig=1,ngrid
1673	pdq(ig,l,iq)=pdq(ig,l,iq)+zdq_source(ig,l,iq)
1674	ENDDO
1675	ENDDO
1676	ENDDO
1677	ENDIF
1678
1679	!------------------------------------------------
1680	! VII. Surface and sub-surface soil temperature
1681	!------------------------------------------------
1682
1683	! For diagnostic
1684	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1685	if (.not.fast) then
1686	DO ig=1,ngrid
1687	rho(ig,1) = pplay(ig,1)/(r*pt(ig,1))
1688	sensiblehf1(ig)=rho(ig,1)cpp(0.4/log(zzlay(ig,1)/z0))*2 &
1689	(pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1))*0.5 &
1690	(tsurf(ig)-pt(ig,1))
1691	if (calldifv) then
1692	sensiblehf2(ig)=zflubid(ig)-capcal(ig)*zdtsdif(ig)
1693	end if
1694	ENDDO
1695	endif
1696
1697
1698	! VII.1 Increment surface temperature
1699	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1700	tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
1701
1702	! Prevent surface (.e.g. non volatile ch4) to exceed max temperature
1703	! Lellouch et al., 2000,2011
1704	IF (tsurfmax) THEN
1705	DO ig=1,ngrid
1706	if (albedo_equivalent(ig).gt.albmin_ch4.and. &
1707	qsurf(ig,igcm_n2).eq.0.) then
1708	tsurf(ig)=min(tsurf(ig),54.)
1709	endif
1710	ENDDO
1711	ENDIF
1712
1713	! VII.2 Compute soil temperatures and subsurface heat flux.
1714	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1715	if (callsoil) then
1716	call soil(ngrid,nsoilmx,.false.,lastcall,therm_inertia, &
1717	ptimestep,tsurf,tsoil,capcal,fluxgrd)
1718	endif
1719
1720	! ! For output :
1721	! tidat_out(:,:)=0.
1722	! DO l=1,nsoilmx
1723	! tidat_out(:,l)=therm_inertia(:,l)
1724	! ENDDO
1725
1726	! Test energy conservation
1727	if(enertest)then
1728	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf(:)/totarea_planet,dEtots)
1729	if (is_master) print*,'Surface energy change =',dEtots,' W m-2'
1730	endif
1731
1732
1733
1734	! VII.3 multiply tendencies of cond/subli for paleo loop only in the
1735	! last Pluto year of the simulation
1736	! Year day must be adapted in the startfi for each object
1737	! Paleo uses year_day to calculate the annual mean tendancies
1738	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1739	IF (paleo) then
1740	if (zday.gt.day_ini+ptime0+nday-year_day) then
1741	DO iq=1,nq
1742	DO ig=1,ngrid
1743	qsurfyear(ig,iq)=qsurfyear(ig,iq)+ &
1744	(qsurf(ig,iq)-qsurf1(ig,iq)) !kg m-2 !ptimestep
1745	ENDDO
1746	ENDDO
1747	endif
1748	endif
1749
1750	! VII.4 Glacial flow at each timestep glastep or at lastcall
1751	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1752	IF (fast.and.glaflow) THEN
1753	if ((mod(zday-day_ini-ptime0,glastep)).lt.1. &
1754	.or.lastcall) then
1755	IF (lastcall) then
1756	dstep=mod(zday-day_ini-ptime0,glastep)*daysec
1757	else
1758	dstep=glastep*daysec
1759	endif
1760	zdqflow(:,:)=qsurf(:,:)
1761	IF (paleo) then
1762	call spreadglacier_paleo(ngrid,nq,qsurf, &
1763	phisfinew,dstep,tsurf)
1764	else
1765	call spreadglacier_simple(ngrid,nq,qsurf,dstep)
1766	endif
1767	zdqflow(:,:)=(zdqflow(:,:)-qsurf(:,:))/dstep
1768
1769	if (conservn2) then
1770	write(,) 'conservn2 glaflow'
1771	call testconservmass(ngrid,nlayer,pplev(:,1)+ &
1772	pdpsrf(:)*ptimestep,qsurf(:,1))
1773	endif
1774
1775	endif
1776	ENDIF
1777
1778	!---------------------------------------------------
1779	! VIII. Perform diagnostics and write output files
1780	!---------------------------------------------------
1781
1782	! Note : For output only: the actual model integration is performed in the dynamics.
1783
1784
1785	! Temperature, zonal and meridional winds.
1786	zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
1787	zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
1788	zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
1789
1790	!! Recast thermal plume vertical velocity array for outputs
1791	!! AF24: removed
1792
1793	! Diagnostic.
1794	zdtdyn(1:ngrid,1:nlayer) = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep
1795	ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
1796
1797	zdudyn(1:ngrid,1:nlayer) = (pu(1:ngrid,1:nlayer)-zuprevious(1:ngrid,1:nlayer)) / ptimestep
1798	zuprevious(1:ngrid,1:nlayer) = zu(1:ngrid,1:nlayer)
1799
1800	if(firstcall)then
1801	zdtdyn(1:ngrid,1:nlayer)=0.0
1802	zdudyn(1:ngrid,1:nlayer)=0.0
1803	endif
1804
1805	! Dynamical heating diagnostic
1806	fluxdyn(:)=0.0
1807	if (.not.fast) then
1808	do ig=1,ngrid
1809	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp
1810	enddo
1811	endif
1812
1813	! Tracers.
1814	zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
1815
1816	! Surface pressure.
1817	ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
1818	call planetwide_sumval(ps(:)*cell_area(:)/totarea_planet,globave)
1819
1820	! pressure density !pluto specific
1821	IF (.not.fast) then !
1822	do ig=1,ngrid
1823	do l=1,nlayer
1824	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1825	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1826	rho(ig,l) = zplay(ig,l)/(r*zt(ig,l))
1827	enddo
1828	zplev(ig,nlayer+1)=pplev(ig,nlayer+1)/pplev(ig,1)*ps(ig)
1829	enddo
1830	ENDIF
1831
1832
1833	! Surface and soil temperature information
1834	call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
1835	call planetwide_minval(tsurf(:),Ts2)
1836	call planetwide_maxval(tsurf(:),Ts3)
1837	if(callsoil)then
1838	TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
1839	if (is_master) then
1840	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
1841	print*,Ts1,Ts2,Ts3,TsS
1842	end if
1843	else
1844	if (is_master) then
1845	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
1846	print*,Ts1,Ts2,Ts3
1847	endif
1848	end if
1849
1850
1851	! Check the energy balance of the simulation during the run
1852	if(corrk)then
1853
1854	call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
1855	call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
1856	call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
1857	call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
1858	call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
1859	do ig=1,ngrid
1860	if(fluxtop_dn(ig).lt.0.0)then
1861	print*,'fluxtop_dn has gone crazy'
1862	print*,'fluxtop_dn=',fluxtop_dn(ig)
1863	print*,'tau_col=',tau_col(ig)
1864	print*,'dtau_aer=',dtau_aer(ig,:,:)
1865	print*,'temp= ',pt(ig,:)
1866	print*,'pplay= ',pplay(ig,:)
1867	call abort
1868	endif
1869	end do
1870
1871	if(ngrid.eq.1)then
1872	DYN=0.0
1873	endif
1874
1875	if (is_master) then
1876	print*,' ISR ASR OLR GND DYN [W m^-2]'
1877	print*, ISR,ASR,OLR,GND,DYN
1878	endif
1879
1880	if(enertest .and. is_master)then
1881	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
1882	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
1883	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
1884	endif
1885
1886	if(meanOLR .and. is_master)then
1887	if((ngrid.gt.1) .or. (mod(icount-1,diagfi_output_rate).eq.0))then
1888	! to record global radiative balance
1889	open(92,file="rad_bal.out",form='formatted',position='append')
1890	write(92,*) zday,ISR,ASR,OLR
1891	close(92)
1892	open(93,file="tem_bal.out",form='formatted',position='append')
1893	if(callsoil)then
1894	write(93,*) zday,Ts1,Ts2,Ts3,TsS
1895	else
1896	write(93,*) zday,Ts1,Ts2,Ts3
1897	endif
1898	close(93)
1899	endif
1900	endif
1901
1902	endif ! end of 'corrk'
1903
1904
1905	! Diagnostic to test radiative-convective timescales in code.
1906	if(testradtimes)then
1907	open(38,file="tau_phys.out",form='formatted',position='append')
1908	ig=1
1909	do l=1,nlayer
1910	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1911	enddo
1912	close(38)
1913	print*,'As testradtimes enabled,'
1914	print*,'exiting physics on first call'
1915	call abort
1916	endif
1917
1918
1919	! Compute column amounts (kg m-2) if tracers are enabled.
1920	if(tracer)then
1921	qcol(1:ngrid,1:nq)=0.0
1922	do iq=1,nq
1923	do ig=1,ngrid
1924	qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
1925	enddo
1926	enddo
1927
1928	endif ! end of 'tracer'
1929
1930	if (methane) then
1931	IF (fast) then ! zq is the mixing ratio supposingly mixed in all atmosphere
1932	DO ig=1,ngrid
1933	vmr_ch4(ig)=zq(ig,1,igcm_ch4_gas)* &
1934	mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.
1935	ENDDO
1936	ELSE
1937	DO ig=1,ngrid
1938	! compute vmr methane
1939	vmr_ch4(ig)=qcol(ig,igcm_ch4_gas)* &
1940	g/ps(ig)mmol(igcm_n2)/mmol(igcm_ch4_gas)100.
1941	! compute density methane
1942	DO l=1,nlayer
1943	zrho_ch4(ig,l)=zq(ig,l,igcm_ch4_gas)*rho(ig,l)
1944	ENDDO
1945	ENDDO
1946	ENDIF
1947	endif
1948
1949	if (carbox) then
1950	IF (fast) then
1951	DO ig=1,ngrid
1952	vmr_co(ig)=zq(ig,1,igcm_co_gas)* &
1953	mmol(igcm_n2)/mmol(igcm_co_gas)*100.
1954	ENDDO
1955	ELSE
1956	DO ig=1,ngrid
1957	! compute vmr CO
1958	vmr_co(ig)=qcol(ig,igcm_co_gas)* &
1959	g/ps(ig)mmol(igcm_n2)/mmol(igcm_co_gas)100.
1960	! compute density CO
1961	DO l=1,nlayer
1962	zrho_co(ig,l)=zq(ig,l,igcm_co_gas)*rho(ig,l)
1963	ENDDO
1964	ENDDO
1965	ENDIF
1966	endif
1967
1968	zrho_haze(:,:)=0.
1969	zdqrho_photprec(:,:)=0.
1970	IF (haze.or.callmufi) then
1971	DO ig=1,ngrid
1972	DO l=1,nlayer
1973	zrho_haze(ig,l)=zq(ig,l,igcm_haze)*rho(ig,l)
1974	zdqrho_photprec(ig,l)=zdqphot_prec(ig,l)*rho(ig,l)
1975	ENDDO
1976	ENDDO
1977	ENDIF
1978
1979	IF (fasthaze) then
1980	DO ig=1,ngrid
1981	qcol(ig,igcm_haze)=zq(ig,1,igcm_haze)*pplev(ig,1)/g
1982	qcol(ig,igcm_prec_haze)=zq(ig,1,igcm_prec_haze)*pplev(ig,1)/g
1983	ENDDO
1984	ENDIF
1985
1986	! Info about Ls, declin...
1987	IF (fast) THEN
1988	if (is_master) write(,)'Ls=',zls180./pi,' dec=',declin180./pi
1989	if (is_master) write(,)'zday=',zday,' ps=',globave
1990	IF (lastcall) then
1991	if (is_master) write(,)'lastcall'
1992	ENDIF
1993	ELSE
1994	if (is_master) write(,)'Ls=',zls180./pi,'decli=',declin180./pi,'zday=',zday
1995	ENDIF
1996
1997	lecttsoil=0 ! default value for lecttsoil
1998	call getin_p("lecttsoil",lecttsoil)
1999	IF (lastcall.and.(ngrid.EQ.1).and.(lecttsoil.eq.1)) THEN
2000	! save tsoil temperature profile for 1D profile
2001	OPEN(13,file='proftsoil.out',form='formatted')
2002	DO i=1,nsoilmx
2003	write(13,*) tsoil(1,i)
2004	ENDDO
2005	CLOSE(13)
2006	ENDIF
2007
2008	if (is_master) print,'--> Ls =',zls180./pi
2009
2010	write_restartfi = .false.
2011	if(lastcall) then
2012	write_restartfi = .true.
2013
2014	if (paleo) then
2015	! time range for tendencies of ice flux qsurfyear
2016	zdt_tot=year_day ! Last year of simulation
2017
2018	masslost(:)=0.
2019	massacc(:)=0.
2020
2021	DO ig=1,ngrid
2022	! update new reservoir of ice on the surface
2023	DO iq=1,nq
2024	! kg/m2 to be sublimed or condensed during paleoyears
2025	qsurfyear(ig,iq)=qsurfyear(ig,iq)* &
2026	paleoyears365.25/(zdt_totdaysec/86400.)
2027
2028	! special case if we sublime the entire reservoir
2029	!! AF: TODO : fix following lines (real_area), using line below:
2030	! call planetwide_sumval((-qsurfyear(:,iq)-qsurf(:,iq))*cell_area(:),masslost)
2031
2032	! IF (-qsurfyear(ig,iq).gt.qsurf(ig,iq)) THEN
2033	! masslost(iq)=masslost(iq)+real_area(ig)* &
2034	! (-qsurfyear(ig,iq)-qsurf(ig,iq))
2035	! qsurfyear(ig,iq)=-qsurf(ig,iq)
2036	! ENDIF
2037
2038	! IF (qsurfyear(ig,iq).gt.0.) THEN
2039	! massacc(iq)=massacc(iq)+real_area(ig)*qsurfyear(ig,iq)
2040	! ENDIF
2041
2042
2043	ENDDO
2044	ENDDO
2045
2046	DO ig=1,ngrid
2047	DO iq=1,nq
2048	qsurfpal(ig,iq)=qsurf(ig,iq)+qsurfyear(ig,iq)
2049	IF (qsurfyear(ig,iq).gt.0.) THEN
2050	qsurfpal(ig,iq)=qsurfpal(ig,iq)- &
2051	qsurfyear(ig,iq)*masslost(iq)/massacc(iq)
2052	ENDIF
2053	ENDDO
2054	ENDDO
2055	! Finally ensure conservation of qsurf
2056	DO iq=1,nq
2057	call planetwide_sumval(qsurf(:,iq)*cell_area(:)/totarea_planet,globaveice(iq))
2058	call planetwide_sumval(qsurfpal(:,iq)*cell_area(:)/totarea_planet,globavenewice(iq))
2059	IF (globavenewice(iq).gt.0.) THEN
2060	qsurfpal(:,iq)=qsurfpal(:,iq)* &
2061	globaveice(iq)/globavenewice(iq)
2062	ENDIF
2063	ENDDO
2064
2065	! update new geopotential depending on the ice reservoir
2066	phisfipal(:)=phisfinew(:)+qsurfpal(:,igcm_n2)*g/1000.
2067	!phisfipal(ig)=phisfi(ig)
2068
2069	if (kbo.or.triton) then ! case of Triton : we do not change the orbital parameters
2070	pdaypal=pday ! no increment of pdaypal to keep same evolution of the subsolar point
2071	eccpal=1.-periastr/((periastr+apoastr)/2.) !no change of ecc
2072	peri_daypal=peri_day ! no change
2073	oblipal=obliquit ! no change
2074	tpalnew=tpal
2075	adjustnew=adjust
2076
2077	else ! Pluto
2078	! update new pday and tpal (Myr) to be set in startfi controle
2079	pdaypal=int(day_ini+paleoyears*365.25/6.3872)
2080	tpalnew=tpal+paleoyears*1.e-6 ! Myrs
2081
2082	! update new N2 ice adjustment (not tested yet on Pluto)
2083	adjustnew=adjust
2084
2085	! update milankovitch parameters : obliquity,Lsp,ecc
2086	call calcmilank(tpalnew,oblipal,peri_daypal,eccpal)
2087	!peri_daypal=peri_day
2088	!eccpal=0.009
2089	endif
2090
2091	if (is_master) write(,) "Paleo peri=",peri_daypal," tpal=",tpalnew
2092	if (is_master) write(,) "Paleo eccpal=",eccpal," tpal=",tpalnew
2093
2094	! create restartfi
2095	if (ngrid.ne.1) then
2096	print*, "physdem1pal not yet implemented"
2097	stop
2098	!TODO: import this routine from pluto.old
2099	! call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq, &
2100	! ptimestep,pdaypal, &
2101	! ztime_restart,tsurf,tsoil,emis,q2,qsurfpal, &
2102	! cell_area,albedodat,therm_inertia,zmea,zstd,zsig, &
2103	! zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal, &
2104	! peri_daypal)
2105	endif
2106	else ! 'paleo'
2107
2108
2109	endif ! end of 'paleo'
2110	endif ! end of 'lastcall'
2111
2112	if(mod(icount,diagfi_output_rate).eq.0) then
2113	write_restartfi = .true.
2114	endif
2115	if(startfi_output_rate.ne.0)then
2116	if (mod(icount,startfi_output_rate).eq.0) then
2117	write_restartfi = .true.
2118	endif
2119	endif
2120	!----------------------------------------------------------------------
2121	! Writing NetCDF file "RESTARTFI" at the end of the run
2122	!----------------------------------------------------------------------
2123	! Note: 'restartfi' is stored just before dynamics are stored
2124	! in 'restart'. Between now and the writing of 'restart',
2125	! there will have been the itau=itau+1 instruction and
2126	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
2127	! thus we store for time=time+dtvr
2128
2129	if (is_omp_master.and.write_restartfi) then
2130	if (ngrid.ne.1) then
2131	IF (grid_type==unstructured) THEN !IF DYNAMICO
2132	! DYNAMICO: no need to add a dynamics time step to ztime_restart
2133	ztime_restart = ptime
2134	ELSE ! LMDZ
2135	ztime_restart = ptime + (ptimestep/daysec)
2136	ENDIF ! of IF (grid_type==unstructured)
2137
2138	if (is_master) write(,)'PHYSIQ: writing restartfi at time =',ztime_restart
2139
2140	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
2141	ptimestep,ztime_restart,tsurf, &
2142	tsoil,therm_inertia,emis,albedo,q2,qsurf)
2143	endif ! ngrid
2144	endif ! is_omp_master
2145
2146
2147
2148	!------------------------------------------------------------------------------
2149	! OUTPUT in netcdf file "DIAGFI.NC",
2150	! containing any variable for diagnostic
2151	!
2152	! Note 1 : output with period "diagfi_output_rate", set in "run.def"
2153	! Note 2 : writediagfi can also be called from any other subroutine
2154	! for any variable, but its preferable to keep all the
2155	! calls in one place ...
2156	!------------------------------------------------------------------------------
2157
2158	!-------- General 1D variables
2159
2160	call write_output("Ls","solar longitude","deg",zls*180./pi)
2161	! call write_output("Lss","sub solar longitude","deg",zlss*180./pi)
2162	call write_output("RA","right ascension","deg",right_ascen*180./pi)
2163	call write_output("Declin","solar declination","deg",declin*180./pi)
2164	call write_output("dist_star","dist_star","AU",dist_star)
2165	call write_output("globave","surf press","Pa",globave)
2166
2167	!-------- General 2D variables
2168
2169	call write_output("tsurf","Surface temperature","K",tsurf)
2170	call write_output("ps","Surface pressure","Pa",ps)
2171	call write_output("emis","Emissivity","",emis)
2172	call write_output('phisfi','Surface Geopotential','m2s-2',phisfi)
2173	if (ngrid.ne.1) then
2174	if (grid_type == regular_lonlat) then
2175	call write_output("area","Mesh area","m2", &
2176	cell_area_for_lonlat_outputs)
2177	else ! unstructured grid (e.g. dynamico)
2178	call write_output("area","Mesh area","m2",cell_area)
2179	endif
2180	endif
2181
2182	if (fast) then
2183	call write_output("fluxrad","fluxrad","W m-2",fluxrad)
2184	! call write_output("dplanck","dplanck","W.s m-2 K-1",dplanck)
2185	endif
2186
2187	if (callsoil) then
2188	! "soil" variables
2189	call write_output("capcal","Surface Heat Capacity","W.s m-2 K-1",capcal)
2190	call write_output("tsoil","tsoil","K",tsoil)
2191	call write_output("therm_inertia","therm_inertia","S.I.",therm_inertia)
2192	endif
2193
2194	! Total energy balance diagnostics
2195	if(callrad)then
2196	call write_output("ALB","Surface albedo"," ",albedo_equivalent)
2197	call write_output("ASR","absorbed stellar rad.","W m-2",fluxabs_sw)
2198	call write_output("ISR","incoming stellar rad.","W m-2",fluxtop_dn)
2199	call write_output("OLR","outgoing longwave rad.","W m-2",fluxtop_lw)
2200	call write_output("GND","heat flux from ground","W m-2",fluxgrd)
2201	if (.not.fast) then
2202	call write_output("DYN","dynamical heat input","W m-2",fluxdyn)
2203	endif
2204	endif ! end of 'callrad'
2205
2206	!-------- General 3D variables
2207
2208	if (.not.fast) then
2209	if (check_physics_outputs) then
2210	! Check the validity of updated fields at the end of the physics step
2211	call check_physics_fields("physiq:", zt, zu, zv, pplev, zq)
2212	endif
2213
2214	call write_output("zzlay","Midlayer altitude", "m",zzlay(:,:))
2215	call write_output("zzlev","Interlayer altitude", "m",zzlev(:,1:nlayer))
2216	call write_output('pphi','Geopotential','m2s-2',pphi)
2217
2218	call write_output("temperature","temperature","K",zt)
2219	call write_output("teta","potential temperature","K",zh)
2220	call write_output("u","Zonal wind","m.s-1",zu)
2221	call write_output("v","Meridional wind","m.s-1",zv)
2222	call write_output("w","Vertical wind","m.s-1",pw)
2223	call write_output("p","Pressure","Pa",pplay)
2224	call write_output("omega","omega","Pa/s",omega)
2225	endif
2226
2227	if(enertest) then
2228	if (calldifv) then
2229	call write_output("q2","turbulent kinetic energy","J.kg^-1",q2)
2230	call write_output("sensibFlux","sensible heat flux","w.m^-2",sensibFlux)
2231	endif
2232	if (corrk) then
2233	call write_output("dEzradsw","radiative heating","w.m^-2",dEzRadsw)
2234	call write_output("dEzradlw","radiative heating","w.m^-2",dEzRadlw)
2235	endif
2236	endif ! end of 'enertest'
2237
2238	! Temporary inclusions for heating diagnostics.
2239	if (.not.fast) then
2240	call write_output("zdtsw","SW heating","K s-1",zdtsw)
2241	call write_output("zdtlw","LW heating","K s-1",zdtlw)
2242	call write_output("dtrad","radiative heating","K s-1",dtrad)
2243	call write_output("zdtdyn","Dyn. heating","K s-1",zdtdyn)
2244	call write_output("zdudyn","Dyn. U","m s-2",zdudyn)
2245	call write_output("zdtconduc","tendancy conduc","K s-1",zdtconduc)
2246	call write_output("zdumolvis","tendancy molvis","m s-1",zdumolvis)
2247	call write_output("zdvmolvis","tendancy molvis","m s-1",zdvmolvis)
2248	call write_output("zdtdif","tendancy T diff","K s-1",zdtdif)
2249	call write_output("zdtsdif","tendancy Ts diff","K s-1",zdtsdif)
2250	call write_output("zdtadj","tendancy T adj","K s-1",zdtadj)
2251	endif
2252
2253	! Optical diagnostics.
2254	call write_output("tau_col",&
2255	"Aerosol surface opacity at reference visible wavelength","",tau_col)
2256	! Diagnostics of optical thickness (dtau = dtau_gas + dtau_rayaer + dtau_cont).
2257	! Warning this is exp(-tau), I let you postproc with -log to have tau itself
2258	call write_output('dtauv_01','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,1)) ! 5.398 um (17x27)
2259	call write_output('dtauv_23','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,23)) ! 0.941 um (17x27)
2260	call write_output('dtauv_24','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,24)) ! 0.700 um (17x27)
2261	call write_output('dtauv_27','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,27)) ! 0.119 um (17x27)
2262	call write_output('dtaui_01','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,1)) ! 550.0 um (17x27)
2263	call write_output('dtaui_17','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,17)) ! 3.531 um (17x27)
2264
2265	if (calllott) then
2266	call write_output("zdugw","tendancy U lott","m s-2",zdugw)
2267	call write_output("zdvgw","tendancy V lott","m s-2",zdvgw)
2268	call write_output("zdtgw","tendancy T lott","K s-1",zdtgw)
2269	endif
2270
2271	! Output tracers.
2272	if (tracer) then
2273
2274	do iq=1,nq
2275	if (.not.fast) then
2276	call write_output(noms(iq),noms(iq),'kg/kg',zq(:,:,iq))
2277	endif
2278	call write_output(trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2279	'kg m^-2',qcol(:,iq) )
2280	call write_output(trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2281	'kg m^-2',qsurf(:,iq) )
2282	enddo ! end of 'nq' loop
2283
2284	! N2 cycle
2285	call write_output('n2_iceflux','n2_iceflux',"kg m^-2 s^-1",flusurf(:,igcm_n2) )
2286	if (.not.fast) then
2287	call write_output("zdtc","tendancy T cond N2","K s-1",zdtc)
2288	call write_output("zdtsurfc","tendancy Ts cond N2","K s-1",zdtsurfc)
2289	call write_output("zduc","tendancy U cond N2","m s-1",zduc)
2290	call write_output("zdvc","tendancy V cond N2","m s-1",zdvc)
2291	call write_output("zdqc_n2","tendancy tracer cond N2","kg kg-1 s-1",zdqc(:,:,1))
2292	call write_output("zdqsc_n2","tendancy tracer surf cond N2","kg kg-1 s-1",zdqsc(:,1))
2293	call write_output("zdqdif_n2","tendancy tracer diff","kg kg-1 s-1",zdqdif(:,:,1))
2294	call write_output("zdqsdif_n2","tendancy tracer surf diff","kg kg-1 s-1",zdqsdif(:,1))
2295	call write_output("zdqadj_n2","tendancy tracer adj","K s-1",zdqadj(:,:,1))
2296	endif
2297
2298	! CH4 cycle
2299	if (methane) then
2300
2301	call write_output('ch4_iceflux','ch4_iceflux',&
2302	"kg m^-2 s^-1",flusurf(:,igcm_ch4_ice) )
2303	call write_output("vmr_ch4","vmr_ch4","%",vmr_ch4)
2304
2305	if (.not.fast) then
2306	call write_output("zrho_ch4","zrho_ch4","kg.m-3",zrho_ch4(:,:))
2307	!call write_output("rice_ch4","ch4 ice mass mean radius","m",rice_ch4)
2308	!call write_output("zq1temp_ch4"," "," ",zq1temp_ch4)
2309	!call write_output("qsat_ch4"," "," ",qsat_ch4)
2310	!call write_output("qsat_ch4_l1"," "," ",qsat_ch4_l1)
2311
2312	! 3D Tendancies
2313	call write_output("zdqcn2_ch4","zdq condn2 ch4","",&
2314	zdqc(:,:,igcm_ch4_gas))
2315	call write_output("zdqdif_ch4","zdqdif ch4","",&
2316	zdqdif(:,:,igcm_ch4_gas))
2317	call write_output("zdqsdif_ch4_ice","zdqsdif ch4","",&
2318	zdqsdif(:,igcm_ch4_ice))
2319	call write_output("zdqadj_ch4","zdqadj ch4","",&
2320	zdqadj(:,:,igcm_ch4_gas))
2321	endif
2322
2323	if (sedimentation) then
2324	call write_output("zdqsed_ch4","zdqsed ch4","",&
2325	zdqsed(:,:,igcm_ch4_gas))
2326	call write_output("zdqssed_ch4","zdqssed ch4","",&
2327	zdqssed(:,igcm_ch4_gas))
2328	endif
2329
2330	if (metcloud.and.(.not.fast)) then
2331	call write_output("zdtch4cloud","ch4 cloud","K s-1",&
2332	zdtch4cloud)
2333	call write_output("zdqch4cloud_gas","ch4 cloud","kg kg-1 s-1",&
2334	zdqch4cloud(:,:,igcm_ch4_gas))
2335	call write_output("zdqch4cloud_ice","ch4 cloud","kg kg-1 s-1",&
2336	zdqch4cloud(:,:,igcm_ch4_ice))
2337	endif
2338
2339	endif
2340
2341	! CO cycle
2342	if (carbox) then
2343	! call write_output("zdtcocloud","tendancy T cocloud","K",zdtcocloud)
2344	call write_output('co_iceflux','co_iceflux',&
2345	"kg m^-2 s^-1",flusurf(:,igcm_co_ice) )
2346	call write_output("vmr_co","vmr_co","%",vmr_co)
2347	if (.not.fast) THEN
2348	call write_output("zrho_co","zrho_co","kg.m-3",zrho_co(:,:))
2349	endif
2350	endif
2351
2352	! Haze or callmufi
2353	if (haze.or.callmufi) then
2354
2355	if (haze_radproffix)then
2356	call write_output('haze_reff','haze_reff','m',reffrad(:,:,1))
2357	end if
2358	!call write_output("zrho_haze","zrho_haze","kg.m-3",zrho_haze(:,:))
2359	!call write_output("zdqhaze_col","zdqhaze col","kg/m2/s",&
2360	! zdqhaze_col(:))
2361
2362	! 3D Tendencies
2363	call write_output("zdqrho_photprec","Photolysis rate",&
2364	"kg.m-3.s-1",zdqrho_photprec(:,:))
2365
2366	!call write_output("zdqhaze_ch4","Photolysis of CH4 tendancy","kg.kg-1.s-1",&
2367	! zdqhaze(:,:,igcm_ch4_gas))
2368	!call write_output("zdqhaze_prec","Precursor haze tendancy","kg.kg-1.s-1",&
2369	! zdqhaze(:,:,igcm_prec_haze))
2370
2371	call write_output("zdqphot_ch4","Photolysis of CH4 tendancy","kg.kg-1.s-1",&
2372	zdqphot_ch4(:,:))
2373	call write_output("zdqphot_prec","Precursor haze tendancy","kg.kg-1.s-1",&
2374	zdqphot_prec(:,:))
2375	call write_output("zdqconv_prec","Conversion of prec_haze into haze tendancy","kg.kg-1.s-1",&
2376	zdqconv_prec(:,:))
2377
2378	if (igcm_haze.ne.0) then
2379	call write_output("zdqhaze_haze","zdqhaze_haze","",&
2380	zdqhaze(:,:,igcm_haze))
2381	if (sedimentation) then
2382	call write_output("zdqssed_haze","zdqssed haze",&
2383	"kg/m2/s",zdqssed(:,igcm_haze))
2384	endif
2385	endif
2386	endif ! end Haze or callmufi
2387
2388	if (callmufi) then
2389	! Tracers:
2390	call write_output("m0as","Density number of spherical aerosols","m-3",zq(:,:,micro_indx(1))*int2ext(:,:))
2391	call write_output("m3as","Volume of spherical aerosols","m3.m-3",zq(:,:,micro_indx(2))*int2ext(:,:))
2392	call write_output("m0af","Density number of fractal aerosols","m-3",zq(:,:,micro_indx(3))*int2ext(:,:))
2393	call write_output("m3af","Volume of fractal aerosols","m3.m-3",zq(:,:,micro_indx(4))*int2ext(:,:))
2394
2395	! Diagnostics:
2396	call write_output("aers_prec","Spherical aerosols precipitations","kg.m-2.s-1",mp2m_aer_s_prec(:))
2397	call write_output("aerf_prec","Fractal aerosols precipitations","kg.m-2.s-1",mp2m_aer_f_prec(:))
2398	call write_output("aers_w","Spherical aerosol settling velocity","m.s-1",mp2m_aer_s_w(:,:))
2399	call write_output("aerf_w","Fractal aerosol settling velocity","m.s-1",mp2m_aer_f_w(:,:))
2400	call write_output("rcs","Characteristic radius of spherical aerosols","m",mp2m_rc_sph(:,:))
2401	call write_output("rcf","Characteristic radius of fractal aerosols","m",mp2m_rc_fra(:,:))
2402	endif ! end callmufi
2403
2404	endif ! end of 'tracer'
2405
2406	! Output spectrum.
2407	if(specOLR.and.corrk)then
2408	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
2409	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
2410	call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W/m^2/cm^-1",3,GSR_nu)
2411	endif
2412
2413	! XIOS outputs
2414	#ifdef CPP_XIOS
2415	! Send fields to XIOS: (NB these fields must also be defined as
2416	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
2417	CALL send_xios_field("controle",tab_cntrl_mod,1)
2418
2419	CALL send_xios_field("ap",ap,1)
2420	CALL send_xios_field("bp",bp,1)
2421	CALL send_xios_field("aps",aps,1)
2422	CALL send_xios_field("bps",bps,1)
2423
2424	if (lastcall.and.is_omp_master) then
2425	if (is_master) write(,) "physiq: call xios_context_finalize"
2426	call xios_context_finalize
2427	endif
2428	#endif
2429
2430	if (check_physics_outputs) then
2431	! Check the validity of updated fields at the end of the physics step
2432	call check_physics_fields("end of physiq:", zt, zu, zv, pplev, zq)
2433	endif
2434
2435	icount=icount+1
2436
2437	end subroutine physiq
2438
2439	end module physiq_mod

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

Download in other formats:

Original Format