Context Navigation

physiq_mod.F90 @ 3952

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

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

Download in other formats:

Original Format