Context Navigation

physiq_mod.F90 @ 3971

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

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

Download in other formats:

Original Format