Context Navigation

physiq_mod.F90 @ 3412

Last change on this file since 3412 was 3412, checked in by afalco, 3 months ago
Pluto PCM: Imported glaciers & conserv mass routines from pluto.old. AF
Property svn:executable set to ``*
File size: 104.4 KB

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

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

Download in other formats:

Original Format