Context Navigation

physiq_mod.F90 @ 3572

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

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

Download in other formats:

Original Format