Context Navigation

physiq_mod.F90 @ 3585

Last change on this file since 3585 was 3585, checked in by debatzbr, 12 days ago
Connecting microphysics to radiative transfer + miscellaneous cleans
Property svn:executable set to ``*
File size: 96.8 KB

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

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

Download in other formats:

Original Format