Context Navigation

physiq_mod.F90 @ 3936

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

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

Download in other formats:

Original Format