Context Navigation

physiq_mod.F90 @ 3482

Last change on this file since 3482 was 3455, checked in by afalco, 15 months ago
Pluto PCM: added conduction & molvis. AF
Property svn:executable set to ``*
File size: 106.1 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,lw_ch4,lw_co,&
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	callconduct,callmolvis,callmoldiff, &
51	aerohaze, corrk, diagdtau, &
52	diurnal, enertest, fat1au, &
53	icetstep, intheat, iradia, kastprof, &
54	lwrite, mass_redistrib, meanOLR, &
55	fast,fasthaze,haze,metcloud,monoxcloud, &
56	n2cond,nearn2cond,noseason_day,conservn2, &
57	convergeps,kbo,triton,paleo,paleoyears,glaflow, &
58	carbox, methane,condmetsurf,condcosurf, &
59	oldplutovdifc,oldplutocorrk,oldplutosedim, &
60	aerohaze,haze_proffix,haze_radproffix, &
61	source_haze, tsurfmax, albmin_ch4, &
62	season, sedimentation,generic_condensation, &
63	specOLR, &
64	startphy_file, testradtimes, &
65	tracer, UseTurbDiff, &
66	global1d, szangle
67	use generic_tracer_index_mod, only: generic_tracer_index
68	use check_fields_mod, only: check_physics_fields
69	use conc_mod, only: rnew, cpnew, ini_conc_mod
70	use phys_state_var_mod
71	use callcorrk_mod, only: callcorrk
72	use callcorrk_pluto_mod, only: callcorrk_pluto
73	use vdifc_mod, only: vdifc
74	use vdifc_pluto_mod, only: vdifc_pluto
75	use turbdiff_mod, only: turbdiff
76	use turb_mod, only : q2,sensibFlux,turb_resolved
77	use mass_redistribution_mod, only: mass_redistribution
78	use condensation_generic_mod, only: condensation_generic
79	use datafile_mod, only: datadir
80	#ifndef MESOSCALE
81	use vertical_layers_mod, only: presnivs, pseudoalt
82	use mod_phys_lmdz_omp_data, ONLY: is_omp_master
83	#else
84	use comm_wrf, only : comm_HR_SW, comm_HR_LW, &
85	comm_CLOUDFRAC,comm_TOTCLOUDFRAC, comm_RH, &
86	comm_HR_DYN, &
87	comm_DQICE,comm_DQVAP,comm_ALBEQ, &
88	comm_FLUXTOP_DN,comm_FLUXABS_SW, &
89	comm_FLUXTOP_LW,comm_FLUXSURF_SW, &
90	comm_FLUXSURF_LW,comm_FLXGRD, &
91	comm_DTRAIN,comm_DTLSC, &
92	comm_LATENT_HF
93
94	#endif
95
96	#ifdef CPP_XIOS
97	use xios_output_mod, only: initialize_xios_output, &
98	update_xios_timestep, &
99	send_xios_field
100	use wxios, only: wxios_context_init, xios_context_finalize
101	#endif
102
103	implicit none
104
105
106	!==================================================================
107	!
108	! Purpose
109	! -------
110	! Central subroutine for all the physics parameterisations in the
111	! universal model. Originally adapted from the Mars LMDZ model.
112	!
113	! The model can be run without or with tracer transport
114	! depending on the value of "tracer" in file "callphys.def".
115	!
116	!
117	! It includes:
118	!
119	! I. Initialization :
120	! I.1 Firstcall initializations.
121	! I.2 Initialization for every call to physiq.
122	!
123	! II. Compute radiative transfer tendencies (longwave and shortwave) :
124	! II.a Option 1 : Call correlated-k radiative transfer scheme.
125	! II.b Option 2 : Call Newtonian cooling scheme.
126	! II.! Option 3 : Atmosphere has no radiative effect.
127	!
128	! III. Vertical diffusion (turbulent mixing) :
129	!
130	! IV. Convection :
131	! IV.a Thermal plume model
132	! IV.b Dry convective adjusment
133	!
134	! V. Condensation and sublimation of gases (currently just N2).
135	!
136	! VI. Tracers
137	! VI.1. Water and water ice.
138	! VI.2 Photochemistry
139	! VI.3. Aerosols and particles.
140	! VI.4. Updates (pressure variations, surface budget).
141	! VI.5. Slab Ocean.
142	! VI.6. Surface Tracer Update.
143	!
144	! VII. Surface and sub-surface soil temperature.
145	!
146	! VIII. Perform diagnostics and write output files.
147	!
148	!
149	! arguments
150	! ---------
151	!
152	! INPUT
153	! -----
154	!
155	! ngrid Size of the horizontal grid.
156	! nlayer Number of vertical layers.
157	! nq Number of advected fields.
158	!
159	! firstcall True at the first call.
160	! lastcall True at the last call.
161	!
162	! pday Number of days counted from the North. Spring equinoxe.
163	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
164	! ptimestep timestep (s).
165	!
166	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa).
167	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
168	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2.s-2).
169	!
170	! pu(ngrid,nlayer) u, zonal component of the wind (ms-1).
171	! pv(ngrid,nlayer) v, meridional component of the wind (ms-1).
172	!
173	! pt(ngrid,nlayer) Temperature (K).
174	!
175	! pq(ngrid,nlayer,nq) Advected fields.
176	!
177	! pudyn(ngrid,nlayer) \
178	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
179	! ptdyn(ngrid,nlayer) / corresponding variables.
180	! pqdyn(ngrid,nlayer,nq) /
181	! flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
182	!
183	! OUTPUT
184	! ------
185	!
186	! pdu(ngrid,nlayer) \
187	! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding
188	! pdt(ngrid,nlayer) / variables due to physical processes.
189	! pdq(ngrid,nlayer) /
190	! pdpsrf(ngrid) /
191	!
192	!
193	! Authors
194	! -------
195	! Frederic Hourdin 15/10/93
196	! Francois Forget 1994
197	! Christophe Hourdin 02/1997
198	! Subroutine completely rewritten by F. Forget (01/2000)
199	! Water ice clouds: Franck Montmessin (update 06/2003)
200	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
201	! New correlated-k radiative scheme: R. Wordsworth (2009)
202	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
203	! Improved water cycle: R. Wordsworth / B. Charnay (2010)
204	! To F90: R. Wordsworth (2010)
205	! New turbulent diffusion scheme: J. Leconte (2012)
206	! Loops converted to F90 matrix format: J. Leconte (2012)
207	! No more ngrid/nq, F90 commons and adaptation to parallel: A. Spiga (2012)
208	! Purge of the code : M. Turbet (2015)
209	! Photochemical core developped by F. Lefevre: B. Charnay (2017)
210	! Purge for Pluto model : A. Falco (2024)
211	!==================================================================
212
213
214	! 0. Declarations :
215	! ------------------
216
217	include "netcdf.inc"
218
219	! Arguments :
220	! -----------
221
222	! INPUTS:
223	! -------
224
225	integer,intent(in) :: ngrid ! Number of atmospheric columns.
226	integer,intent(in) :: nlayer ! Number of atmospheric layers.
227	integer,intent(in) :: nq ! Number of tracers.
228
229	logical,intent(in) :: firstcall ! Signals first call to physics.
230	logical,intent(in) :: lastcall ! Signals last call to physics.
231
232	real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0.
233	real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
234	real,intent(in) :: ptimestep ! Physics timestep (s).
235	real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
236	real,intent(in) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
237	real,intent(in) :: pphi(ngrid,nlayer) ! Geopotential at mid-layer (m2s-2).
238	real,intent(in) :: pu(ngrid,nlayer) ! Zonal wind component (m/s).
239	real,intent(in) :: pv(ngrid,nlayer) ! Meridional wind component (m/s).
240	real,intent(in) :: pt(ngrid,nlayer) ! Temperature (K).
241	real,intent(in) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
242	real,intent(in) :: flxw(ngrid,nlayer) ! Vertical mass flux (ks/s) at lower boundary of layer
243
244	! OUTPUTS:
245	! --------
246
247	! Physical tendencies :
248
249	real,intent(out) :: pdu(ngrid,nlayer) ! Zonal wind tendencies (m/s/s).
250	real,intent(out) :: pdv(ngrid,nlayer) ! Meridional wind tendencies (m/s/s).
251	real,intent(out) :: pdt(ngrid,nlayer) ! Temperature tendencies (K/s).
252	real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
253	real,intent(out) :: pdpsrf(ngrid) ! Surface pressure tendency (Pa/s).
254
255	! Local saved variables:
256	! ----------------------
257	integer,save :: day_ini ! Initial date of the run (sol since Ls=0).
258	integer,save :: icount ! Counter of calls to physiq during the run.
259	!$OMP THREADPRIVATE(day_ini,icount)
260
261	!Pluto specific
262	REAL,save :: acond,bcond
263	REAL,save :: tcond1p4Pa
264	DATA tcond1p4Pa/38/
265	real :: tidat(ngrid,nsoilmx) ! thermal inertia soil
266
267
268
269	! Local variables :
270	! -----------------
271	! Tendencies for the paleoclimate mode
272	REAL qsurfyear(ngrid,nq) ! kg.m-2 averaged mass of ice lost/gained in the last Pluto year of the run
273	REAL phisfinew(ngrid) ! geopotential of the bedrock (= phisfi-qsurf/1000*g)
274	REAL qsurfpal(ngrid,nq) ! qsurf after a paleoclimate step : for physdem1 and restartfi
275	REAL phisfipal(ngrid) ! geopotential after a paleoclimate step : for physdem1 and restartfi
276	REAL oblipal ! change of obliquity
277	REAL peri_daypal ! new periday
278	REAL eccpal ! change of eccentricity
279	REAL tpalnew ! change of time
280	REAL adjustnew ! change in N2 ice albedo
281	REAL pdaypal ! new pday = day_ini + step
282	REAL zdt_tot ! time range corresponding to the flux of qsurfyear
283	REAL massacc(nq) ! accumulated mass
284	REAL masslost(nq) ! accumulated mass
285
286	REAL globave ! globalaverage 2D ps
287	REAL globaveice(nq) ! globalaverage 2D ice
288	REAL globavenewice(nq) ! globalaverage 2D ice
289	INTEGER lecttsoil ! lecture of tsoil from proftsoil
290	REAL qsurf1(ngrid,nq) ! saving qsurf to calculate flux over long timescales kg.m-2
291	REAL flusurf(ngrid,nq) ! flux cond/sub kg.m-2.s-1
292	REAL flusurfold(ngrid,nq) ! old flux cond/sub kg.m-2.s-1
293	REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
294
295	REAL,SAVE :: ptime0 ! store the first time
296	REAL dstep
297	REAL,SAVE :: glastep=20 ! step in pluto day to spread glacier
298
299
300	! Aerosol (dust or ice) extinction optical depth at reference wavelength
301	! for the "naerkind" optically active aerosols:
302
303	real,save,allocatable :: aerosol(:,:,:) ! Aerosols
304	!$OMP THREADPRIVATE(aerosol)
305	real zh(ngrid,nlayer) ! Potential temperature (K).
306	real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
307	real omega(ngrid,nlayer) ! omega velocity (Pa/s, >0 when downward)
308
309	integer i,l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_gas, igcm_generic_ice
310	logical call_ice_gas_generic ! to call only one time the ice/vap pair of a tracer
311
312	real zls ! Solar longitude (radians).
313	real zlss ! Sub solar point longitude (radians).
314	real zday ! Date (time since Ls=0, calculated in sols).
315	REAL,save :: saveday ! saved date
316	REAL,save :: savedeclin ! saved declin
317	real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers.
318	real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries.
319
320
321	! VARIABLES for the thermal plume model (AF24: deleted)
322
323	! TENDENCIES due to various processes :
324
325	! For Surface Temperature : (K/s)
326	real zdtsurf(ngrid) ! Cumulated tendencies.
327	real zdtsurfmr(ngrid) ! Mass_redistribution routine.
328	real zdtsurfc(ngrid) ! Condense_n2 routine.
329	real zdtsdif(ngrid) ! Turbdiff/vdifc routines.
330	! real zdtsurf_hyd(ngrid) ! Hydrol routine.
331
332	! For Atmospheric Temperatures : (K/s)
333	real dtlscale(ngrid,nlayer) ! Largescale routine.
334	real dt_generic_condensation(ngrid,nlayer) ! condensation_generic routine.
335	real zdtc(ngrid,nlayer) ! Condense_n2 routine.
336	real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
337	real zdtmr(ngrid,nlayer) ! Mass_redistribution routine.
338	real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
339	real zdtchim(ngrid,nlayer) ! Calchim routine.
340
341	! For Surface Tracers : (kg/m2/s)
342	real dqsurf(ngrid,nq) ! Cumulated tendencies.
343	!real zdqsurfc(ngrid) ! Condense_n2 routine.
344	REAL zdqsc(ngrid,nq) ! Condense_n2 routine.
345	real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines.
346	real zdqssed(ngrid,nq) ! Callsedim routine.
347	real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine.
348
349	! For Tracers : (kg/kg_of_air/s)
350	real zdqc(ngrid,nlayer,nq) ! Condense_n2 routine.
351	real zdqadj(ngrid,nlayer,nq) ! Convadj routine.
352	real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines.
353	real zdqevap(ngrid,nlayer) ! Turbdiff routine.
354	real zdqsed(ngrid,nlayer,nq) ! Callsedim routine.
355	real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine.
356	real dqvaplscale_generic(ngrid,nlayer,nq) ! condensation_generic routine.
357	real dqcldlscale_generic(ngrid,nlayer,nq) ! condensation_generic routine.
358	REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine
359	REAL,allocatable,save :: zdqschim(:,:) ! Calchim_asis routine
360	!$OMP THREADPRIVATE(zdqchim,zdqschim)
361
362
363	!! PLUTO variables
364	REAL zdqch4cloud(ngrid,nlayer,nq)
365	REAL zdqsch4cloud(ngrid,nq)
366	REAL zdtch4cloud(ngrid,nlayer)
367	REAL zdqcocloud(ngrid,nlayer,nq)
368	REAL zdqscocloud(ngrid,nq)
369	REAL zdtcocloud(ngrid,nlayer)
370	REAL rice_ch4(ngrid,nlayer) ! Methane ice geometric mean radius (m)
371	REAL rice_co(ngrid,nlayer) ! CO ice geometric mean radius (m)
372
373	REAL zdqsch4fast(ngrid) ! used only for fast model nogcm
374	REAL zdqch4fast(ngrid) ! used only for fast model nogcm
375	REAL zdqscofast(ngrid) ! used only for fast model nogcm
376	REAL zdqcofast(ngrid) ! used only for fast model nogcm
377	REAL zdqflow(ngrid,nq)
378
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
871	! --------------------------------------------------------
872	! 1.3 thermosphere
873	! --------------------------------------------------------
874
875	! ajout de la conduction depuis la thermosphere
876	IF (callconduct) THEN
877
878	call conduction (ngrid,nlayer,ptimestep, &
879	pplay,pt,zzlay,zzlev,zdtconduc,tsurf)
880	DO l=1,nlayer
881	DO ig=1,ngrid
882	pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
883	ENDDO
884	ENDDO
885
886	ENDIF
887
888	! ajout de la viscosite moleculaire
889	IF (callmolvis) THEN
890	call molvis(ngrid,nlayer,ptimestep, &
891	pplay,pt,zzlay,zzlev, &
892	zdtconduc,pu,tsurf,zdumolvis)
893	call molvis(ngrid,nlayer,ptimestep, &
894	pplay,pt,zzlay,zzlev, &
895	zdtconduc,pv,tsurf,zdvmolvis)
896
897	DO l=1,nlayer
898	DO ig=1,ngrid
899	! pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
900	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
901	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
902	ENDDO
903	ENDDO
904	ENDIF
905
906	IF (callmoldiff) THEN
907	call moldiff_red(ngrid,nlayer,nq, &
908	pplay,pplev,pt,pdt,pq,pdq,ptimestep, &
909	zzlay,zdtconduc,zdqmoldiff)
910
911	DO l=1,nlayer
912	DO ig=1,ngrid
913	DO iq=1, nq
914	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
915	ENDDO
916	ENDDO
917	ENDDO
918	ENDIF
919
920	if (conservn2) then
921	write(,) 'conservn2 thermo'
922	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
923	endif
924
925
926	!---------------------------------
927	! II. Compute radiative tendencies
928	!---------------------------------
929	! Saving qsurf to compute paleo flux condensation/sublimation
930	DO iq=1, nq
931	DO ig=1,ngrid
932	IF (qsurf(ig,iq).lt.0.) then
933	qsurf(ig,iq)=0.
934	ENDIF
935	qsurf1(ig,iq)=qsurf(ig,iq)
936	ENDDO
937	ENDDO
938
939
940	! Compute local stellar zenith angles
941	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
942	fract = 0
943	if (diurnal) then
944	ztim1=SIN(declin)
945	ztim2=COS(declin)COS(2.pi*(zday-.5))
946	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
947
948	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
949	ztim1,ztim2,ztim3,mu0,fract)
950	else if(diurnal .eqv. .false.) then
951
952	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad)
953	! WARNING: this function appears not to work in 1D
954
955	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
956	mu0 = cos(pi*szangle/180.0)
957	fract= 1/(4*mu0) ! AF24: from pluto.old
958	endif
959
960	endif
961
962
963	! Pluto albedo /IT changes depending on surface ices (only in 3D)
964	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
965	if (ngrid.ne.1) then
966
967	!! Specific to change albedo of N2 so that Psurf
968	!! converges toward 1.4 Pa in "1989" seasons for Triton
969	!! converges toward 1.1 Pa in "2015" seasons for Pluto
970	if (convergeps) then
971	if (triton) then
972	! 1989 declination
973	if (declin180./pi.gt.-46..and.declin180./pi.lt.-45. &
974	.and.zday.gt.saveday+1000. &
975	.and.declin.lt.savedeclin) then
976	call globalaverage2d(ngrid,pplev(:,1),globave)
977	if (globave.gt.1.5) then
978	adjust=adjust+0.005
979	else if (globave.lt.1.3) then
980	adjust=adjust-0.005
981	endif
982	saveday=zday
983	endif
984	else
985	! Pluto : 2015 declination current epoch
986	if (declin180./pi.gt.50.and.declin180./pi.lt.51. &
987	.and.zday.gt.saveday+10000. &
988	.and.declin.gt.savedeclin) then
989	call globalaverage2d(ngrid,pplev(:,1),globave)
990	if (globave.gt.1.2) then
991	adjust=adjust+0.005
992	else if (globave.lt.1.) then
993	adjust=adjust-0.005
994	endif
995	saveday=zday
996	endif
997	endif
998	end if
999	end if ! if ngrid ne 1
1000
1001	call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat, &
1002	capcal,adjust,dist_star,albedo,emis,flusurfold,ptimestep, &
1003	zls)
1004	! do i=2,ngrid
1005	! albedo(i,:) = albedo(1,:)
1006	! enddo
1007
1008	if (callrad) then
1009	if( mod(icount-1,iradia).eq.0.or.lastcall) then
1010
1011	! Eclipse incoming sunlight !AF24: removed
1012
1013	!! call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1014	!! call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1015
1016
1017	if (corrk) then
1018
1019	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1020	! II.a Call correlated-k radiative transfer scheme
1021	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1022	if(kastprof)then
1023	print*,'kastprof should not = true here'
1024	call abort
1025	endif
1026	! if(water) then !AF24: removed
1027
1028	if(generic_condensation) then
1029	do iq=1,nq
1030
1031	call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
1032
1033	if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
1034
1035	epsi_generic=constants_epsi_generic(iq)
1036
1037	muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_gas))
1038	muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_gas))
1039
1040	endif
1041	end do ! do iq=1,nq loop on tracers
1042	! take into account generic condensable specie (GCS) effect on mean molecular weight
1043
1044	else
1045	muvar(1:ngrid,1:nlayer+1)=mugaz
1046	endif
1047
1048	! if(ok_slab_ocean) then !AF24: removed
1049
1050	! standard callcorrk
1051	if (oldplutocorrk) then
1052	call callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf, &
1053	albedo(:,1),emis,mu0,pplev,pplay,pt, &
1054	zzlay,tsurf,fract,dist_star,aerosol, &
1055	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, &
1056	fluxabs_sw,fluxtop_dn,reffrad,tau_col,ptime,pday, &
1057	cloudfrac,totcloudfrac,.false., &
1058	firstcall,lastcall)
1059	albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
1060	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+ &
1061	fluxsurf_sw(1:ngrid)*(1.-albedo(1:ngrid,1))
1062	else
1063	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
1064	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
1065	zzlay,tsurf,fract,dist_star,aerosol,muvar, &
1066	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, &
1067	fluxsurfabs_sw,fluxtop_lw, &
1068	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,GSR_nu, &
1069	int_dtaui,int_dtauv, &
1070	tau_col,cloudfrac,totcloudfrac, &
1071	.false.,firstcall,lastcall)
1072	! Radiative flux from the sky absorbed by the surface (W.m-2).
1073	GSR=0.0
1074	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+ &
1075	fluxsurfabs_sw(1:ngrid)
1076	endif ! oldplutocorrk
1077	!GG (feb2021): Option to "artificially" decrease the raditive time scale in
1078	!the deep atmosphere press > 0.1 bar. Suggested by MT.
1079	!! COEFF_RAD to be "tuned" to facilitate convergence of tendency
1080
1081	!coeff_rad=0. ! 0 values, it doesn't accelerate the convergence
1082	!coeff_rad=0.5
1083	!coeff_rad=1.
1084	!do l=1, nlayer
1085	! do ig=1,ngrid
1086	! if(pplay(ig,l).ge.1.d4) then
1087	! zdtsw(ig,l)=zdtsw(ig,l)(pplay(ig,l)/1.d4)*coeff_rad
1088	! zdtlw(ig,l)=zdtlw(ig,l)(pplay(ig,l)/1.d4)*coeff_rad
1089	! endif
1090	! enddo
1091	!enddo
1092
1093	! AF24: removed CLFvarying Option
1094
1095
1096	!if(noradsurf)then ! no lower surface; SW flux just disappears
1097	! GSR = SUM(fluxsurf_sw(1:ngrid)*cell_area(1:ngrid))/totarea
1098	! fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
1099	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
1100	!endif
1101
1102	! Net atmospheric radiative heating rate (K.s-1)
1103	dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
1104
1105	! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that !
1106	if (firstcall .and. albedo_spectral_mode) then
1107	call spectral_albedo_calc(albedo_snow_SPECTV)
1108	endif
1109
1110	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1111	! II.b Call Newtonian cooling scheme !AF24: removed
1112	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1113	else
1114	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1115	! II.! Atmosphere has no radiative effect
1116	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1117	fluxtop_dn(1:ngrid) = fract(1:ngrid)mu0(1:ngrid)Fat1AU/dist_star**2
1118	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
1119	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
1120	endif
1121	fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
1122	print*,'------------WARNING---WARNING------------' ! by MT2015.
1123	print*,'You are in corrk=false mode, '
1124	print*,'and the surface albedo is taken equal to the first visible spectral value'
1125
1126	albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
1127	fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
1128	fluxabs_sw(1:ngrid)=fluxsurfabs_sw(1:ngrid)
1129	fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid)
1130	fluxtop_lw(1:ngrid) = emis(1:ngrid)sigmatsurf(1:ngrid)**4
1131
1132	dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating
1133
1134	endif ! end of corrk
1135
1136	endif ! of if(mod(icount-1,iradia).eq.0)
1137
1138
1139	! Transformation of the radiative tendencies
1140	! ------------------------------------------
1141	zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
1142	zplanck(1:ngrid)=emis(1:ngrid)sigmazplanck(1:ngrid)*zplanck(1:ngrid)
1143	fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
1144	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
1145
1146	! Test of energy conservation
1147	!----------------------------
1148	if(enertest)then
1149	call planetwide_sumval(cppmassarea(:,:)zdtsw(:,:)/totarea_planet,dEtotSW)
1150	call planetwide_sumval(cppmassarea(:,:)zdtlw(:,:)/totarea_planet,dEtotLW)
1151	!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
1152	call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
1153	call planetwide_sumval((fluxsurf_lw(:)emis(:)-zplanck(:))cell_area(:)/totarea_planet,dEtotsLW)
1154	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
1155	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
1156	if (is_master) then
1157	print*,'---------------------------------------------------------------'
1158	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
1159	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
1160	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
1161	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
1162	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
1163	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
1164	endif
1165	endif ! end of 'enertest'
1166
1167	endif ! of if (callrad)
1168
1169	!! call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1170	!! call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1171
1172	if (conservn2) then
1173	write(,) 'conservn2 radiat'
1174	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
1175	endif
1176
1177	! --------------------------------------------
1178	! III. Vertical diffusion (turbulent mixing) :
1179	! --------------------------------------------
1180
1181	if (calldifv) then
1182
1183	zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
1184
1185	if (oldplutovdifc) then
1186	zdum1(:,:) = 0
1187	zdum2(:,:) = 0
1188	zdh(:,:)=pdt(:,:)/zpopsk(:,:)
1189
1190	! Calling vdif (Martian version WITH N2 condensation)
1191	CALL vdifc_pluto(ngrid,nlayer,nq,zpopsk, &
1192	ptimestep,capcal,lwrite, &
1193	pplay,pplev,zzlay,zzlev,z0, &
1194	pu,pv,zh,pq,pt,tsurf,emis,qsurf, &
1195	zdum1,zdum2,zdh,pdq,pdt,zflubid, &
1196	zdudif,zdvdif,zdhdif,zdtsdif,q2, &
1197	zdqdif,zdqsdif,qsat_ch4,qsat_ch4_l1) !,zq1temp_ch4,qsat_ch4)
1198
1199	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1200
1201	bcond=1./tcond1p4Pa
1202	acond=r/lw_n2
1203
1204	!------------------------------------------------------------------
1205	! test methane conservation
1206	! if(methane)then
1207	! call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
1208	! & ptimestep,pplev,zdqdif,zdqsdif,'CH4',' vdifc ')
1209	! endif ! methane
1210	!------------------------------------------------------------------
1211	! test CO conservation
1212	! if(carbox)then
1213	! call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
1214	! & ptimestep,pplev,zdqdif,zdqsdif,'CO ',' vdifc ')
1215	! endif ! carbox
1216	!------------------------------------------------------------------
1217
1218	! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
1219	else if (UseTurbDiff) then
1220
1221	call turbdiff(ngrid,nlayer,nq, &
1222	ptimestep,capcal, &
1223	pplay,pplev,zzlay,zzlev,z0, &
1224	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
1225	pdt,pdq,zflubid, &
1226	zdudif,zdvdif,zdtdif,zdtsdif, &
1227	sensibFlux,q2,zdqdif,zdqevap,zdqsdif, &
1228	taux,tauy)
1229
1230	else ! if (oldplutovdifc) .and. (UseTurbDiff)
1231
1232	zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1233
1234	call vdifc(ngrid,nlayer,nq,zpopsk, &
1235	ptimestep,capcal,lwrite, &
1236	pplay,pplev,zzlay,zzlev,z0, &
1237	pu,pv,zh,pq,tsurf,emis,qsurf, &
1238	zdh,pdq,zflubid, &
1239	zdudif,zdvdif,zdhdif,zdtsdif, &
1240	sensibFlux,q2,zdqdif,zdqsdif)
1241
1242	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1243	zdqevap(1:ngrid,1:nlayer)=0.
1244
1245	end if !end of 'UseTurbDiff'
1246
1247	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
1248
1249	!!! this is always done, except for turbulence-resolving simulations
1250	if (.not. turb_resolved) then
1251	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
1252	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
1253	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
1254	endif
1255
1256	! if(ok_slab_ocean)then !AF24: removed
1257	! flux_sens_lat(1:ngrid)=(zdtsdif(1:ngrid)*capcal(1:ngrid)-fluxrad(1:ngrid))
1258	! endif
1259
1260	!! call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_gas))
1261
1262	if (tracer) then
1263	pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
1264	dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
1265	end if ! of if (tracer)
1266
1267	!! call writediagfi(ngrid,"vdifc_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1268	!! call writediagfi(ngrid,"vdifc_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1269
1270	! test energy conservation
1271	!-------------------------
1272	if(enertest)then
1273
1274	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
1275	do ig = 1, ngrid
1276	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
1277	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
1278	enddo
1279
1280	call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
1281	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
1282	call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
1283	call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
1284
1285	if (is_master) then
1286
1287	if (UseTurbDiff) then
1288	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1289	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
1290	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1291	else
1292	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1293	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
1294	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1295	end if
1296	endif ! end of 'is_master'
1297
1298	! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
1299	endif ! end of 'enertest'
1300
1301
1302	! ! Test water conservation. !AF24: removed
1303
1304	else ! calldifv
1305
1306	! if(.not.newtonian)then
1307	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
1308
1309	! ------------------------------------------------------------------
1310	! Methane surface sublimation and condensation in fast model (nogcm)
1311	! ------------------------------------------------------------------
1312	if ((methane).and.(fast).and.condmetsurf) THEN
1313
1314	call ch4surf(ngrid,nlayer,nq,ptimestep, &
1315	tsurf,zdtsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf, &
1316	zdqch4fast,zdqsch4fast)
1317
1318	dqsurf(1:ngrid,igcm_ch4_ice)= dqsurf(1:ngrid,igcm_ch4_ice) + &
1319	zdqsch4fast(1:ngrid)
1320	pdq(1:ngrid,1,igcm_ch4_gas)= pdq(1:ngrid,1,igcm_ch4_gas) + &
1321	zdqch4fast(1:ngrid)
1322	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+lw_ch4*zdqsch4fast(1:ngrid)/capcal(1:ngrid)
1323	end if
1324	! ------------------------------------------------------------------
1325	! CO surface sublimation and condensation in fast model (nogcm)
1326	! ------------------------------------------------------------------
1327	if ((carbox).and.(fast).and.condcosurf) THEN
1328
1329	call cosurf(ngrid,nlayer,nq,ptimestep, &
1330	tsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf, &
1331	zdqcofast,zdqscofast)
1332
1333	dqsurf(1:ngrid,igcm_co_ice)= dqsurf(1:ngrid,igcm_co_ice) + &
1334	zdqscofast(1:ngrid)
1335	pdq(1:ngrid,1,igcm_co_gas)= pdq(1:ngrid,1,igcm_co_gas) + &
1336	zdqcofast(1:ngrid)
1337	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+lw_co*zdqscofast(1:ngrid)/capcal(1:ngrid)
1338	end if
1339
1340
1341	endif ! end of 'calldifv'
1342
1343	if (conservn2) then
1344	write(,) 'conservn2 diff'
1345	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)+ &
1346	dqsurf(:,1)*ptimestep)
1347	endif
1348
1349	!-------------------
1350	! IV. Convection :
1351	!-------------------
1352
1353	! ~~~~~~~~~~~~~~~~~~~~~~~~~~
1354	! IV.a Thermal plume model : AF24: removed
1355	! ~~~~~~~~~~~~~~~~~~~~~~~~~~
1356
1357	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1358	! IV.b Dry convective adjustment :
1359	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1360
1361	if(calladj) then
1362
1363	zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1364	zduadj(1:ngrid,1:nlayer)=0.0
1365	zdvadj(1:ngrid,1:nlayer)=0.0
1366	zdhadj(1:ngrid,1:nlayer)=0.0
1367
1368
1369	call convadj(ngrid,nlayer,nq,ptimestep, &
1370	pplay,pplev,zpopsk, &
1371	pu,pv,zh,pq, &
1372	pdu,pdv,zdh,pdq, &
1373	zduadj,zdvadj,zdhadj, &
1374	zdqadj)
1375
1376	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
1377	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer)
1378	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
1379	zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1380
1381	if(tracer) then
1382	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
1383	end if
1384
1385	! Test energy conservation
1386	if(enertest)then
1387	call planetwide_sumval(cppmassarea(:,:)zdtadj(:,:)/totarea_planet,dEtot)
1388	if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
1389	endif
1390
1391	! ! Test water conservation !AF24: removed
1392
1393	endif ! end of 'calladj'
1394	!----------------------------------------------
1395	! Non orographic Gravity Waves: AF24: removed
1396	!---------------------------------------------
1397
1398	!-----------------------------------------------
1399	! V. Nitrogen condensation-sublimation :
1400	!-----------------------------------------------
1401
1402	if (n2cond) then
1403
1404	if (.not.tracer) then
1405	print*,'We need a N2 ice tracer to condense N2'
1406	call abort
1407	endif
1408
1409	zdqc(:,:,:)=0.
1410	zdqsc(:,:)=0.
1411	call condense_n2(ngrid,nlayer,nq,ptimestep, &
1412	capcal,pplay,pplev,tsurf,pt, &
1413	pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, &
1414	qsurf(1,igcm_n2),albedo,emis, &
1415	zdtc,zdtsurfc,pdpsrf,zduc,zdvc, &
1416	zdqc,zdqsc(1,igcm_n2))
1417
1418	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
1419	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer)+zdvc(1:ngrid,1:nlayer)
1420	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer)+zduc(1:ngrid,1:nlayer)
1421	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
1422
1423	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq)
1424	dqsurf(1:ngrid,igcm_n2) = dqsurf(1:ngrid,igcm_n2) + zdqsc(1:ngrid,igcm_n2)
1425
1426	!! call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1427	!! call writediagfi(ngrid,"condense_n2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1428
1429	! test energy conservation
1430	if(enertest)then
1431	call planetwide_sumval(cppmassarea(:,:)zdtc(:,:)/totarea_planet,dEtot)
1432	call planetwide_sumval(capcal(:)zdtsurfc(:)cell_area(:)/totarea_planet,dEtots)
1433	if (is_master) then
1434	print*,'In n2cloud atmospheric energy change =',dEtot,' W m-2'
1435	print*,'In n2cloud surface energy change =',dEtots,' W m-2'
1436	endif
1437	endif
1438
1439	endif ! end of 'n2cond'
1440
1441	if (conservn2) then
1442	write(,) 'conservn2 n2cond'
1443	call testconservmass(ngrid,nlayer,pplev(:,1)+ &
1444	pdpsrf(:)ptimestep,qsurf(:,1)+dqsurf(:,1)ptimestep)
1445	endif
1446
1447	!---------------------------------------------
1448	! VI. Specific parameterizations for tracers
1449	!---------------------------------------------
1450
1451	if (tracer) then
1452
1453
1454
1455	! 7a. Methane, CO, and ice
1456	! ---------------------------------------
1457	! Methane ice condensation in the atmosphere
1458	! ----------------------------------------
1459	rice_ch4(:,:)=0 ! initialization needed for callsedim
1460	zdqch4cloud(:,:,:)=0.
1461	if ((methane).and.(metcloud).and.(.not.fast)) THEN
1462	call ch4cloud(ngrid,nlayer,naerkind,ptimestep, &
1463	pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, &
1464	pq,pdq,zdqch4cloud,zdqsch4cloud,zdtch4cloud, &
1465	nq,rice_ch4)
1466
1467	DO l=1,nlayer
1468	DO ig=1,ngrid
1469	pdq(ig,l,igcm_ch4_gas)=pdq(ig,l,igcm_ch4_gas)+ &
1470	zdqch4cloud(ig,l,igcm_ch4_gas)
1471	pdq(ig,l,igcm_ch4_ice)=pdq(ig,l,igcm_ch4_ice)+ &
1472	zdqch4cloud(ig,l,igcm_ch4_ice)
1473	ENDDO
1474	ENDDO
1475
1476	! Increment methane ice surface tracer tendency
1477	DO ig=1,ngrid
1478	dqsurf(ig,igcm_ch4_ice)=dqsurf(ig,igcm_ch4_ice)+ &
1479	zdqsch4cloud(ig,igcm_ch4_ice)
1480	ENDDO
1481
1482	! update temperature tendancy
1483	DO ig=1,ngrid
1484	DO l=1,nlayer
1485	pdt(ig,l)=pdt(ig,l)+zdtch4cloud(ig,l)
1486	ENDDO
1487	ENDDO
1488	end if
1489
1490	! ---------------------------------------
1491	! CO ice condensation in the atmosphere
1492	! ----------------------------------------
1493	zdqcocloud(:,:,:)=0.
1494	IF ((carbox).and.(monoxcloud).and.(.not.fast)) THEN
1495	call cocloud(ngrid,nlayer,naerkind,ptimestep, &
1496	pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, &
1497	pq,pdq,zdqcocloud,zdqscocloud,zdtcocloud, &
1498	nq,rice_co,qsurf(1,igcm_n2),dqsurf(1,igcm_n2))
1499
1500	DO l=1,nlayer
1501	DO ig=1,ngrid
1502	pdq(ig,l,igcm_co_gas)=pdq(ig,l,igcm_co_gas)+ &
1503	zdqcocloud(ig,l,igcm_co_gas)
1504	pdq(ig,l,igcm_co_ice)=pdq(ig,l,igcm_co_ice)+ &
1505	zdqcocloud(ig,l,igcm_co_ice)
1506	ENDDO
1507	ENDDO
1508
1509	! Increment CO ice surface tracer tendency
1510	DO ig=1,ngrid
1511	dqsurf(ig,igcm_co_ice)=dqsurf(ig,igcm_co_ice)+ &
1512	zdqscocloud(ig,igcm_co_ice)
1513	ENDDO
1514
1515	! update temperature tendancy
1516	DO ig=1,ngrid
1517	DO l=1,nlayer
1518	pdt(ig,l)=pdt(ig,l)+zdtcocloud(ig,l)
1519	ENDDO
1520	ENDDO
1521	ELSE
1522	rice_co(:,:)=0 ! initialization needed for callsedim
1523	END IF ! of IF (carbox)
1524
1525	!------------------------------------------------------------------
1526	! test methane conservation
1527	! if(methane)then
1528	! call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
1529	! & ptimestep,pplev,zdqch4cloud,zdqsch4cloud,'CH4','ch4clou')
1530	! endif ! methane
1531	!------------------------------------------------------------------
1532	! test CO conservation
1533	! if(carbox)then
1534	! call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
1535	! & ptimestep,pplev,zdqcocloud,zdqscocloud,'CO ','cocloud')
1536	! endif ! carbox
1537	!------------------------------------------------------------------
1538
1539	! 7b. Haze particle production
1540	! -------------------
1541	IF (haze) THEN
1542
1543	zdqphot_prec(:,:)=0.
1544	zdqphot_ch4(:,:)=0.
1545	zdqhaze(:,:,:)=0
1546	! Forcing to a fixed haze profile if haze_proffix
1547	if (haze_proffix.and.i_haze.gt.0.) then
1548	call haze_prof(ngrid,nlayer,zzlay,pplay,pt, &
1549	reffrad,profmmr)
1550	zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze)) &
1551	/ptimestep
1552	else
1553	call hazecloud(ngrid,nlayer,nq,ptimestep, &
1554	pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,zdqhaze, &
1555	zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
1556	endif
1557
1558	DO iq=1, nq ! should be updated
1559	DO l=1,nlayer
1560	DO ig=1,ngrid
1561	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqhaze(ig,l,iq)
1562	ENDDO
1563	ENDDO
1564	ENDDO
1565
1566	ENDIF
1567
1568	IF (fast.and.fasthaze) THEN
1569	call prodhaze(ngrid,nlayer,nq,ptimestep,pplev,pq,pdq,dist_star, &
1570	mu0,declin,zdqprodhaze,zdqsprodhaze,gradflux,fluxbot, &
1571	fluxlym_sol_bot,fluxlym_ipm_bot,flym_sol,flym_ipm)
1572
1573	DO ig=1,ngrid
1574	pdq(ig,1,igcm_ch4_gas)=pdq(ig,1,igcm_ch4_gas)+ &
1575	zdqprodhaze(ig,igcm_ch4_gas)
1576	pdq(ig,1,igcm_prec_haze)=pdq(ig,1,igcm_prec_haze)+ &
1577	zdqprodhaze(ig,igcm_prec_haze)
1578	pdq(ig,1,igcm_haze)=abs(pdq(ig,1,igcm_haze)+ &
1579	zdqprodhaze(ig,igcm_haze))
1580	qsurf(ig,igcm_haze)= qsurf(ig,igcm_haze)+ &
1581	zdqsprodhaze(ig)*ptimestep
1582	ENDDO
1583
1584	ENDIF
1585
1586
1587
1588	! -------------------------
1589	! VI.3. Aerosol particles
1590	! -------------------------
1591
1592	!Generic Condensation
1593	if (generic_condensation) then
1594	call condensation_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay, &
1595	pt,pq,pdt,pdq,dt_generic_condensation, &
1596	dqvaplscale_generic,dqcldlscale_generic,rneb_generic)
1597	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dt_generic_condensation(1:ngrid,1:nlayer)
1598	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+dqvaplscale_generic(1:ngrid,1:nlayer,1:nq)
1599	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+dqcldlscale_generic(1:ngrid,1:nlayer,1:nq)
1600
1601	if(enertest)then
1602	do ig=1,ngrid
1603	genericconddE(ig) = cppSUM(mass(:,:)dt_generic_condensation(:,:))
1604	enddo
1605
1606	call planetwide_sumval(cppmassarea(:,:)dt_generic_condensation(:,:)/totarea_planet,dEtot)
1607
1608	if (is_master) print*,'In generic condensation atmospheric energy change =',dEtot,' W m-2'
1609	end if
1610
1611	! Let's loop on tracers
1612	cloudfrac(:,:)=0.0
1613	do iq=1,nq
1614	call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
1615	if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
1616	do l = 1, nlayer
1617	do ig=1,ngrid
1618	cloudfrac(ig,l)=rneb_generic(ig,l,iq)
1619	enddo
1620	enddo
1621	endif
1622	end do ! do iq=1,nq loop on tracers
1623	! endif ! .not. water
1624
1625	endif !generic_condensation
1626
1627	! Sedimentation.
1628	if (sedimentation) then
1629
1630	zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
1631	zdqssed(1:ngrid,1:nq) = 0.0
1632
1633	if (oldplutosedim)then
1634	call callsedim_pluto(ngrid,nlayer,ptimestep, &
1635	pplev,zzlev,pt,pdt,rice_ch4,rice_co, &
1636	pq,pdq,zdqsed,zdqssed,nq,pphi)
1637
1638	else
1639	call callsedim(ngrid,nlayer,ptimestep, &
1640	pplev,zzlev,pt,pdt,pq,pdq, &
1641	zdqsed,zdqssed,nq,pphi)
1642	endif
1643
1644	! Whether it falls as rain or snow depends only on the surface temperature
1645	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
1646	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
1647
1648	end if ! end of 'sedimentation'
1649
1650	! ---------------
1651	! VI.4. Updates
1652	! ---------------
1653
1654	! Updating Atmospheric Mass and Tracers budgets.
1655	if(mass_redistrib) then
1656
1657	zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * 0
1658	! ( zdqevap(1:ngrid,1:nlayer) &
1659	! ! + zdqrain(1:ngrid,1:nlayer,igcm_h2o_gas) &
1660	! ! + dqmoist(1:ngrid,1:nlayer,igcm_h2o_gas) &
1661	! + dqvaplscale(1:ngrid,1:nlayer) )
1662
1663	do ig = 1, ngrid
1664	zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
1665	enddo
1666
1667	! call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
1668	! call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
1669	call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
1670
1671	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
1672	capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
1673	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
1674	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
1675
1676	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
1677	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
1678	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
1679	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
1680	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
1681	pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
1682	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
1683
1684	endif
1685
1686	! call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1687
1688	!! call writediagfi(ngrid,"slab_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1689	!! call writediagfi(ngrid,"slab_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1690
1691
1692	! ------------------
1693	! VI.5. Slab Ocean !AF24: removed
1694	! ------------------
1695
1696	! -----------------------------
1697	! VI.6. Surface Tracer Update
1698	! -----------------------------
1699
1700	! AF24: deleted hydrology
1701
1702	qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
1703
1704	! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
1705	! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
1706	qsurf_hist(:,:) = qsurf(:,:)
1707
1708	! if(ice_update)then
1709	! ice_min(1:ngrid)=min(ice_min(1:ngrid),qsurf(1:ngrid,igcm_h2o_ice))
1710	! endif
1711
1712	endif! end of if 'tracer'
1713
1714	if (conservn2) then
1715	write(,) 'conservn2 tracer'
1716	call testconservmass(ngrid,nlayer,pplev(:,1)+ &
1717	pdpsrf(:)*ptimestep,qsurf(:,1))
1718	endif
1719
1720	DO ig=1,ngrid
1721	flusurf(ig,igcm_n2)=(qsurf(ig,igcm_n2)- &
1722	qsurf1(ig,igcm_n2))/ptimestep
1723	flusurfold(ig,igcm_n2)=flusurf(ig,igcm_n2)
1724	if (methane) then
1725	flusurf(ig,igcm_ch4_ice)=(qsurf(ig,igcm_ch4_ice)- &
1726	qsurf1(ig,igcm_ch4_ice))/ptimestep
1727	flusurfold(ig,igcm_ch4_ice)=flusurf(ig,igcm_ch4_ice)
1728	endif
1729	if (carbox) then
1730	flusurf(ig,igcm_co_ice)=(qsurf(ig,igcm_co_ice)- &
1731	qsurf1(ig,igcm_co_ice))/ptimestep
1732	!flusurfold(ig,igcm_co_ice)=flusurf(ig,igcm_co_ice)
1733	endif
1734	ENDDO
1735
1736	!! Special source of haze particle !
1737	! todo: should be placed in haze and use tendency of n2 instead of flusurf
1738	IF (source_haze) THEN
1739	write(,) "Source haze not supported yet."
1740	stop
1741	! call hazesource(ngrid,nlayer,nq,ptimestep, &
1742	! pplev,flusurf,mu0,zdq_source)
1743
1744	DO iq=1, nq
1745	DO l=1,nlayer
1746	DO ig=1,ngrid
1747	pdq(ig,l,iq)=pdq(ig,l,iq)+zdq_source(ig,l,iq)
1748	ENDDO
1749	ENDDO
1750	ENDDO
1751	ENDIF
1752
1753
1754	!------------------------------------------------
1755	! VII. Surface and sub-surface soil temperature
1756	!------------------------------------------------
1757
1758	! For diagnostic
1759	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1760	DO ig=1,ngrid
1761	rho(ig,1) = pplay(ig,1)/(r*pt(ig,1))
1762	sensiblehf1(ig)=rho(ig,1)cpp(0.4/log(zzlay(ig,1)/z0))*2 &
1763	(pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1))*0.5 &
1764	(tsurf(ig)-pt(ig,1))
1765	if (calldifv) then
1766	sensiblehf2(ig)=zflubid(ig)-capcal(ig)*zdtsdif(ig)
1767	end if
1768
1769	ENDDO
1770
1771
1772	! VII.1 Increment surface temperature
1773	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1774	tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
1775
1776	! Prevent surface (.e.g. non volatile ch4) to exceed max temperature
1777	! Lellouch et al., 2000,2011
1778	IF (tsurfmax) THEN
1779	DO ig=1,ngrid
1780	if (albedo_equivalent(ig).gt.albmin_ch4.and. &
1781	qsurf(ig,igcm_n2).eq.0.) then
1782	tsurf(ig)=min(tsurf(ig),54.)
1783	endif
1784	ENDDO
1785	ENDIF
1786
1787	! VII.2 Compute soil temperatures and subsurface heat flux.
1788	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1789	if (callsoil) then
1790	call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, &
1791	ptimestep,tsurf,tsoil,capcal,fluxgrd)
1792	endif
1793
1794	! ! For output :
1795	! tidat_out(:,:)=0.
1796	! DO l=1,min(nlayermx,nsoilmx)
1797	! tidat_out(:,l)=tidat(:,l)
1798	! ENDDO
1799
1800	! Test energy conservation
1801	if(enertest)then
1802	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf(:)/totarea_planet,dEtots)
1803	if (is_master) print*,'Surface energy change =',dEtots,' W m-2'
1804	endif
1805
1806
1807
1808	! VII.3 multiply tendencies of cond/subli for paleo loop only in the
1809	! last Pluto year of the simulation
1810	! Year day must be adapted in the startfi for each object
1811	! Paleo uses year_day to calculate the annual mean tendancies
1812	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1813	IF (paleo) then
1814	if (zday.gt.day_ini+ptime0+nday-year_day) then
1815	DO iq=1,nq
1816	DO ig=1,ngrid
1817	qsurfyear(ig,iq)=qsurfyear(ig,iq)+ &
1818	(qsurf(ig,iq)-qsurf1(ig,iq)) !kg m-2 !ptimestep
1819	ENDDO
1820	ENDDO
1821	endif
1822	endif
1823
1824	! VII.4 Glacial flow at each timestep glastep or at lastcall
1825	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1826	IF (fast.and.glaflow) THEN
1827	if ((mod(zday-day_ini-ptime0,glastep)).lt.1. &
1828	.or.lastcall) then
1829	IF (lastcall) then
1830	dstep=mod(zday-day_ini-ptime0,glastep)*daysec
1831	else
1832	dstep=glastep*daysec
1833	endif
1834	zdqflow(:,:)=qsurf(:,:)
1835	IF (paleo) then
1836	call spreadglacier_paleo(ngrid,nq,qsurf, &
1837	phisfinew,dstep,tsurf)
1838	else
1839	call spreadglacier_simple(ngrid,nq,qsurf,dstep)
1840	endif
1841	zdqflow(:,:)=(zdqflow(:,:)-qsurf(:,:))/dstep
1842
1843	if (conservn2) then
1844	write(,) 'conservn2 glaflow'
1845	call testconservmass(ngrid,nlayer,pplev(:,1)+ &
1846	pdpsrf(:)*ptimestep,qsurf(:,1))
1847	endif
1848
1849	endif
1850	ENDIF
1851
1852	!---------------------------------------------------
1853	! VIII. Perform diagnostics and write output files
1854	!---------------------------------------------------
1855
1856	! Note : For output only: the actual model integration is performed in the dynamics.
1857
1858
1859	! Temperature, zonal and meridional winds.
1860	zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
1861	zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
1862	zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
1863
1864	!! Recast thermal plume vertical velocity array for outputs
1865	!! AF24: removed
1866
1867	! Diagnostic.
1868	zdtdyn(1:ngrid,1:nlayer) = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep
1869	ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
1870
1871	zdudyn(1:ngrid,1:nlayer) = (pu(1:ngrid,1:nlayer)-zuprevious(1:ngrid,1:nlayer)) / ptimestep
1872	zuprevious(1:ngrid,1:nlayer) = zu(1:ngrid,1:nlayer)
1873
1874	if(firstcall)then
1875	zdtdyn(1:ngrid,1:nlayer)=0.0
1876	zdudyn(1:ngrid,1:nlayer)=0.0
1877	endif
1878
1879	! Dynamical heating diagnostic.
1880	do ig=1,ngrid
1881	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp
1882	enddo
1883
1884	! Tracers.
1885	zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
1886
1887	! Surface pressure.
1888	ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
1889	call globalaverage2d(ngrid,ps,globave)
1890
1891	! pressure density !pluto specific
1892	IF (.not.fast) then !
1893	do ig=1,ngrid
1894	do l=1,nlayer
1895	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1896	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1897	rho(ig,l) = zplay(ig,l)/(r*zt(ig,l))
1898	enddo
1899	zplev(ig,nlayer+1)=pplev(ig,nlayer+1)/pplev(ig,1)*ps(ig)
1900	enddo
1901	ENDIF
1902
1903
1904	! Surface and soil temperature information
1905	call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
1906	call planetwide_minval(tsurf(:),Ts2)
1907	call planetwide_maxval(tsurf(:),Ts3)
1908	if(callsoil)then
1909	TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
1910	if (is_master) then
1911	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
1912	print*,Ts1,Ts2,Ts3,TsS
1913	end if
1914	else
1915	if (is_master) then
1916	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
1917	print*,Ts1,Ts2,Ts3
1918	endif
1919	end if
1920
1921
1922	! Check the energy balance of the simulation during the run
1923	if(corrk)then
1924
1925	call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
1926	call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
1927	call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
1928	call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
1929	call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
1930	do ig=1,ngrid
1931	if(fluxtop_dn(ig).lt.0.0)then
1932	print*,'fluxtop_dn has gone crazy'
1933	print*,'fluxtop_dn=',fluxtop_dn(ig)
1934	print*,'tau_col=',tau_col(ig)
1935	print*,'aerosol=',aerosol(ig,:,:)
1936	print*,'temp= ',pt(ig,:)
1937	print*,'pplay= ',pplay(ig,:)
1938	call abort
1939	endif
1940	end do
1941
1942	if(ngrid.eq.1)then
1943	DYN=0.0
1944	endif
1945
1946	if (is_master) then
1947	print*,' ISR ASR OLR GND DYN [W m^-2]'
1948	print*, ISR,ASR,OLR,GND,DYN
1949	endif
1950
1951	if(enertest .and. is_master)then
1952	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
1953	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
1954	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
1955	endif
1956
1957	if(meanOLR .and. is_master)then
1958	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
1959	! to record global radiative balance
1960	open(92,file="rad_bal.out",form='formatted',position='append')
1961	write(92,*) zday,ISR,ASR,OLR
1962	close(92)
1963	open(93,file="tem_bal.out",form='formatted',position='append')
1964	if(callsoil)then
1965	write(93,*) zday,Ts1,Ts2,Ts3,TsS
1966	else
1967	write(93,*) zday,Ts1,Ts2,Ts3
1968	endif
1969	close(93)
1970	endif
1971	endif
1972
1973	endif ! end of 'corrk'
1974
1975
1976	! Diagnostic to test radiative-convective timescales in code.
1977	if(testradtimes)then
1978	open(38,file="tau_phys.out",form='formatted',position='append')
1979	ig=1
1980	do l=1,nlayer
1981	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1982	enddo
1983	close(38)
1984	print*,'As testradtimes enabled,'
1985	print*,'exiting physics on first call'
1986	call abort
1987	endif
1988
1989
1990	! Compute column amounts (kg m-2) if tracers are enabled.
1991	if(tracer)then
1992	qcol(1:ngrid,1:nq)=0.0
1993	do iq=1,nq
1994	do ig=1,ngrid
1995	qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
1996	enddo
1997	enddo
1998
1999	endif ! end of 'tracer'
2000
2001
2002
2003	! Calculate RH_generic (Generic Relative Humidity) for diagnostic.
2004	if(generic_condensation)then
2005	RH_generic(:,:,:)=0.0
2006	do iq=1,nq
2007
2008	call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
2009
2010	if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
2011
2012	do l = 1, nlayer
2013	do ig=1,ngrid
2014	call Psat_generic(zt(ig,l),pplay(ig,l),metallicity,psat_tmp_generic,qsat_generic(ig,l,iq))
2015	RH_generic(ig,l,iq) = zq(ig,l,igcm_generic_gas) / qsat_generic(ig,l,iq)
2016	enddo
2017	enddo
2018
2019	end if
2020
2021	end do ! iq=1,nq
2022
2023	endif ! end of 'generic_condensation'
2024
2025
2026	if (methane) then
2027	IF (fast) then ! zq is the mixing ratio supposingly mixed in all atmosphere
2028	DO ig=1,ngrid
2029	vmr_ch4(ig)=zq(ig,1,igcm_ch4_gas)* &
2030	mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.
2031	ENDDO
2032	ELSE
2033	DO ig=1,ngrid
2034	! compute vmr methane
2035	vmr_ch4(ig)=qcol(ig,igcm_ch4_gas)* &
2036	g/ps(ig)mmol(igcm_n2)/mmol(igcm_ch4_gas)100.
2037	! compute density methane
2038	DO l=1,nlayer
2039	zrho_ch4(ig,l)=zq(ig,l,igcm_ch4_gas)*rho(ig,l)
2040	ENDDO
2041	ENDDO
2042	ENDIF
2043	endif
2044
2045	if (carbox) then
2046	IF (fast) then
2047	DO ig=1,ngrid
2048	vmr_co(ig)=zq(ig,1,igcm_co_gas)* &
2049	mmol(igcm_n2)/mmol(igcm_co_gas)*100.
2050	ENDDO
2051	ELSE
2052	DO ig=1,ngrid
2053	! compute vmr CO
2054	vmr_co(ig)=qcol(ig,igcm_co_gas)* &
2055	g/ps(ig)mmol(igcm_n2)/mmol(igcm_co_gas)100.
2056	! compute density CO
2057	DO l=1,nlayer
2058	zrho_co(ig,l)=zq(ig,l,igcm_co_gas)*rho(ig,l)
2059	ENDDO
2060	ENDDO
2061	ENDIF
2062	endif
2063
2064	zrho_haze(:,:)=0.
2065	zdqrho_photprec(:,:)=0.
2066	IF (haze.and.aerohaze) then
2067	DO ig=1,ngrid
2068	DO l=1,nlayer
2069	zrho_haze(ig,l)=zq(ig,l,igcm_haze)*rho(ig,l)
2070	zdqrho_photprec(ig,l)=zdqphot_prec(ig,l)*rho(ig,l)
2071	ENDDO
2072	ENDDO
2073	ENDIF
2074
2075	IF (fasthaze) then
2076	DO ig=1,ngrid
2077	qcol(ig,igcm_haze)=zq(ig,1,igcm_haze)*pplev(ig,1)/g
2078	qcol(ig,igcm_prec_haze)=zq(ig,1,igcm_prec_haze)*pplev(ig,1)/g
2079	ENDDO
2080	ENDIF
2081
2082	! Info about Ls, declin...
2083	IF (fast) THEN
2084	if (is_master) write(,),'Ls=',zls180./pi,' dec=',declin180./pi
2085	if (is_master) write(,),'zday=',zday,' ps=',globave
2086	IF (lastcall) then
2087	if (is_master) write(,),'lastcall'
2088	ENDIF
2089	ELSE
2090	if (is_master) write(,),'Ls=',zls180./pi,'decli=',declin180./pi,'zday=',zday
2091	ENDIF
2092
2093	lecttsoil=0 ! default value for lecttsoil
2094	call getin_p("lecttsoil",lecttsoil)
2095	IF (lastcall.and.(ngrid.EQ.1).and.(lecttsoil.eq.1)) THEN
2096	! save tsoil temperature profile for 1D profile
2097	OPEN(13,file='proftsoil.out',form='formatted')
2098	DO i=1,nsoilmx
2099	write(13,*) tsoil(1,i)
2100	ENDDO
2101	CLOSE(13)
2102	ENDIF
2103
2104
2105	if (is_master) print,'--> Ls =',zls180./pi
2106
2107
2108	!----------------------------------------------------------------------
2109	! Writing NetCDF file "RESTARTFI" at the end of the run
2110	!----------------------------------------------------------------------
2111
2112	! Note: 'restartfi' is stored just before dynamics are stored
2113	! in 'restart'. Between now and the writting of 'restart',
2114	! there will have been the itau=itau+1 instruction and
2115	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
2116	! thus we store for time=time+dtvr
2117
2118
2119
2120	if(lastcall) then
2121	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
2122
2123	!! Update surface ice distribution to iterate to steady state if requested
2124	!! AF24: removed
2125
2126	! endif
2127	if (paleo) then
2128	! time range for tendencies of ice flux qsurfyear
2129	zdt_tot=year_day ! Last year of simulation
2130
2131	masslost(:)=0.
2132	massacc(:)=0.
2133
2134	DO ig=1,ngrid
2135	! update new reservoir of ice on the surface
2136	DO iq=1,nq
2137	! kg/m2 to be sublimed or condensed during paleoyears
2138	qsurfyear(ig,iq)=qsurfyear(ig,iq)* &
2139	paleoyears365.25/(zdt_totdaysec/86400.)
2140
2141	! special case if we sublime the entire reservoir
2142	!! AF: TODO : fix following lines (real_area), using line below:
2143	! call planetwide_sumval((-qsurfyear(:,iq)-qsurf(:,iq))*cell_area(:),masslost)
2144
2145	! IF (-qsurfyear(ig,iq).gt.qsurf(ig,iq)) THEN
2146	! masslost(iq)=masslost(iq)+real_area(ig)* &
2147	! (-qsurfyear(ig,iq)-qsurf(ig,iq))
2148	! qsurfyear(ig,iq)=-qsurf(ig,iq)
2149	! ENDIF
2150
2151	! IF (qsurfyear(ig,iq).gt.0.) THEN
2152	! massacc(iq)=massacc(iq)+real_area(ig)*qsurfyear(ig,iq)
2153	! ENDIF
2154
2155
2156	ENDDO
2157	ENDDO
2158
2159	DO ig=1,ngrid
2160	DO iq=1,nq
2161	qsurfpal(ig,iq)=qsurf(ig,iq)+qsurfyear(ig,iq)
2162	IF (qsurfyear(ig,iq).gt.0.) THEN
2163	qsurfpal(ig,iq)=qsurfpal(ig,iq)- &
2164	qsurfyear(ig,iq)*masslost(iq)/massacc(iq)
2165	ENDIF
2166	ENDDO
2167	ENDDO
2168	! Finally ensure conservation of qsurf
2169	DO iq=1,nq
2170	call globalaverage2d(ngrid,qsurf(:,iq),globaveice(iq))
2171	call globalaverage2d(ngrid,qsurfpal(:,iq), &
2172	globavenewice(iq))
2173	IF (globavenewice(iq).gt.0.) THEN
2174	qsurfpal(:,iq)=qsurfpal(:,iq)* &
2175	globaveice(iq)/globavenewice(iq)
2176	ENDIF
2177	ENDDO
2178
2179	! update new geopotential depending on the ice reservoir
2180	phisfipal(:)=phisfinew(:)+qsurfpal(:,igcm_n2)*g/1000.
2181	!phisfipal(ig)=phisfi(ig)
2182
2183
2184	if (kbo.or.triton) then ! case of Triton : we do not change the orbital parameters
2185
2186	pdaypal=pday ! no increment of pdaypal to keep same evolution of the subsolar point
2187	eccpal=1.-periastr/((periastr+apoastr)/2.) !no change of ecc
2188	peri_daypal=peri_day ! no change
2189	oblipal=obliquit ! no change
2190	tpalnew=tpal
2191	adjustnew=adjust
2192
2193	else ! Pluto
2194	! update new pday and tpal (Myr) to be set in startfi controle
2195	pdaypal=int(day_ini+paleoyears*365.25/6.3872)
2196	tpalnew=tpal+paleoyears*1.e-6 ! Myrs
2197
2198	! update new N2 ice adjustment (not tested yet on Pluto)
2199	adjustnew=adjust
2200
2201	! update milankovitch parameters : obliquity,Lsp,ecc
2202	call calcmilank(tpalnew,oblipal,peri_daypal,eccpal)
2203	!peri_daypal=peri_day
2204	!eccpal=0.009
2205
2206	endif
2207
2208	if (is_master) write(,) "Paleo peri=",peri_daypal," tpal=",tpalnew
2209	if (is_master) write(,) "Paleo eccpal=",eccpal," tpal=",tpalnew
2210
2211
2212	#ifndef MESOSCALE
2213	! create restartfi
2214	if (ngrid.ne.1) then
2215	!TODO: import this routine from pluto.old
2216	! call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq, &
2217	! ptimestep,pdaypal, &
2218	! ztime_fin,tsurf,tsoil,emis,q2,qsurfpal, &
2219	! cell_area,albedodat,tidat,zmea,zstd,zsig, &
2220	! zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal, &
2221	! peri_daypal)
2222	endif
2223	else ! 'paleo'
2224
2225	if (ngrid.ne.1) then
2226	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
2227
2228	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
2229	ptimestep,ztime_fin, &
2230	tsurf,tsoil,emis,q2,qsurf_hist)
2231	endif
2232	#endif
2233	endif ! end of 'paleo'
2234	! if(ok_slab_ocean) then
2235	! call ocean_slab_final!(tslab, seaice)
2236	! end if
2237
2238	endif ! end of 'lastcall'
2239
2240
2241	! -----------------------------------------------------------------
2242	! WSTATS: Saving statistics
2243	! -----------------------------------------------------------------
2244	! ("stats" stores and accumulates key variables in file "stats.nc"
2245	! which can later be used to make the statistic files of the run:
2246	! if flag "callstats" from callphys.def is .true.)
2247
2248
2249	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
2250	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2251	call wstats(ngrid,"fluxsurf_lw", &
2252	"Thermal IR radiative flux to surface","W.m-2",2, &
2253	fluxsurf_lw)
2254	call wstats(ngrid,"fluxtop_lw", &
2255	"Thermal IR radiative flux to space","W.m-2",2, &
2256	fluxtop_lw)
2257
2258	! call wstats(ngrid,"fluxsurf_sw", &
2259	! "Solar radiative flux to surface","W.m-2",2, &
2260	! fluxsurf_sw_tot)
2261	! call wstats(ngrid,"fluxtop_sw", &
2262	! "Solar radiative flux to space","W.m-2",2, &
2263	! fluxtop_sw_tot)
2264
2265
2266	call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
2267	call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
2268	call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
2269	!call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
2270	!call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
2271	call wstats(ngrid,"p","Pressure","Pa",3,pplay)
2272	call wstats(ngrid,"emis","Emissivity","",2,emis)
2273	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
2274	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
2275	call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
2276	call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
2277	call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
2278
2279	if (tracer) then
2280	do iq=1,nq
2281	call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
2282	call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2283	'kg m^-2',2,qsurf(1,iq) )
2284	call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2285	'kg m^-2',2,qcol(1,iq) )
2286
2287	! call wstats(ngrid,trim(noms(iq))//'_reff', &
2288	! trim(noms(iq))//'_reff', &
2289	! 'm',3,reffrad(1,1,iq))
2290
2291	end do
2292
2293	endif ! end of 'tracer'
2294
2295	!AF24: deleted slab ocean and water
2296
2297	if(lastcall.and.callstats) then
2298	write (,) "Writing stats..."
2299	call mkstats(ierr)
2300	endif
2301
2302
2303	#ifndef MESOSCALE
2304
2305	!-----------------------------------------------------------------------------------------------------
2306	! OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
2307	!
2308	! Note 1 : output with period "ecritphy", set in "run.def"
2309	!
2310	! Note 2 : writediagfi can also be called from any other subroutine for any variable,
2311	! but its preferable to keep all the calls in one place ...
2312	!-----------------------------------------------------------------------------------------------------
2313
2314	call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
2315	call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
2316	call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
2317	call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
2318	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2319	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
2320	call writediagfi(ngrid,"emis","Emissivity","",2,emis)
2321
2322	!! Pluto outputs
2323	call writediagfi(ngrid,"dist_star","dist_star","AU",0,dist_star)
2324
2325	if (fast) then
2326	call writediagfi(ngrid,"globave","surf press","Pa",0,globave)
2327	call writediagfi(ngrid,"fluxrad","fluxrad","W m-2",2,fluxrad)
2328	call writediagfi(ngrid,"fluxgrd","fluxgrd","W m-2",2,fluxgrd)
2329	call writediagfi(ngrid,"capcal","capcal","W.s m-2 K-1",2,capcal)
2330	! call writediagfi(ngrid,"dplanck","dplanck","W.s m-2 K-1",2,dplanck)
2331	call writediagfi(ngrid,"tsoil","tsoil","K",3,tsoil)
2332	else
2333	if (check_physics_outputs) then
2334	! Check the validity of updated fields at the end of the physics step
2335	call check_physics_fields("HERE physiq:", zt, zu, zv, pplev, zq)
2336	endif
2337
2338	call writediagfi(ngrid,"temp","temperature","K",3,zt)
2339	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
2340	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
2341	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
2342	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
2343	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
2344	endif
2345
2346	! Subsurface temperatures
2347	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2348	! call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
2349
2350	! Total energy balance diagnostics
2351	if(callrad)then
2352
2353	call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
2354	call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
2355	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
2356	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
2357	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
2358	call writediagfi(ngrid,"shad","rings"," ", 2, fract)
2359
2360	! call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
2361	! call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
2362
2363	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
2364	! endif
2365
2366	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
2367
2368	endif ! end of 'callrad'
2369
2370	if(enertest) then
2371	if (calldifv) then
2372
2373	call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
2374	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
2375
2376	! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
2377	! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
2378	! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
2379
2380	endif
2381
2382	if (corrk) then
2383	call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
2384	call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
2385	endif
2386
2387	if (generic_condensation) then
2388	call writediagfi(ngrid,"genericconddE","heat from generic condensation","W m-2",2,genericconddE)
2389	call writediagfi(ngrid,"dt_generic_condensation","heating from generic condensation","K s-1",3,dt_generic_condensation)
2390	endif
2391
2392	endif ! end of 'enertest'
2393
2394	! Diagnostics of optical thickness
2395	! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
2396	if (diagdtau) then
2397	do nw=1,L_NSPECTV
2398	write(str2,'(i2.2)') nw
2399	call writediagfi(ngrid,'dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',1,int_dtauv(:,nlayer:1:-1,nw))
2400	enddo
2401	do nw=1,L_NSPECTI
2402	write(str2,'(i2.2)') nw
2403	call writediagfi(ngrid,'dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',1,int_dtaui(:,nlayer:1:-1,nw))
2404	enddo
2405	endif
2406
2407	! Temporary inclusions for heating diagnostics.
2408	call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
2409	call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
2410	call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
2411	call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
2412
2413	! For Debugging.
2414	!call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
2415	!call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
2416
2417	! Output tracers.
2418	if (tracer) then
2419	! call writediagfi(ngrid,"zdtc","tendancy T cond N2","K",3,zdtc)
2420
2421	do iq=1,nq
2422	call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
2423	! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2424	! 'kg m^-2',2,qsurf_hist(1,iq) )
2425	call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2426	'kg m^-2',2,qcol(1,iq) )
2427	call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2428	'kg m^-2',2,qsurf(1,iq) )
2429
2430	if(generic_condensation)then
2431	call writediagfi(ngrid,"rneb_generic","GCS cloud fraction (generic condensation)"," ",3,rneb_generic)
2432	call writediagfi(ngrid,"CLF","GCS cloud fraction"," ",3,cloudfrac)
2433	call writediagfi(ngrid,"RH_generic","GCS relative humidity"," ",3,RH_generic)
2434	endif
2435	! call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
2436	enddo ! end of 'nq' loop
2437
2438	!Pluto specific
2439	call writediagfi(ngrid,'n2_iceflux','n2_iceflux',"kg m^-2 s^-1",2,flusurf(1,igcm_n2) )
2440	if (haze_radproffix)then
2441	call writediagfi(ngrid,'haze_reff','haze_reff','m',3,reffrad(1,1,1))
2442	end if
2443	if (methane) then
2444	! call writediagfi(ngrid,"rice_ch4","ch4 ice mass mean radius","m",3,rice_ch4)
2445	! call writediagfi(ngrid,"zq1temp_ch4"," "," ",2,zq1temp_ch4)
2446	! call writediagfi(ngrid,"qsat_ch4"," "," ",2,qsat_ch4)
2447	! call writediagfi(ngrid,"qsat_ch4_l1"," "," ",2,qsat_ch4_l1)
2448
2449	call writediagfi(ngrid,'ch4_iceflux','ch4_iceflux',&
2450	"kg m^-2 s^-1",2,flusurf(1,igcm_ch4_ice) )
2451	call writediagfi(ngrid,"vmr_ch4","vmr_ch4","%",2,vmr_ch4)
2452	if (.not.fast) then
2453	call writediagfi(ngrid,"zrho_ch4","zrho_ch4","kg.m-3",3,zrho_ch4(:,:))
2454	endif
2455
2456	! Tendancies
2457	call writediagfi(ngrid,"zdqcn2_ch4","zdq condn2 ch4","",&
2458	3,zdqc(:,:,igcm_ch4_gas))
2459	call writediagfi(ngrid,"zdqdif_ch4","zdqdif ch4","",&
2460	3,zdqdif(:,:,igcm_ch4_gas))
2461	call writediagfi(ngrid,"zdqsdif_ch4_ice","zdqsdif ch4","",&
2462	2,zdqsdif(:,igcm_ch4_ice))
2463	call writediagfi(ngrid,"zdqadj_ch4","zdqadj ch4","",&
2464	3,zdqadj(:,:,igcm_ch4_gas))
2465	if (sedimentation) then
2466	call writediagfi(ngrid,"zdqsed_ch4","zdqsed ch4","",&
2467	3,zdqsed(:,:,igcm_ch4_gas))
2468	call writediagfi(ngrid,"zdqssed_ch4","zdqssed ch4","",&
2469	2,zdqssed(:,igcm_ch4_gas))
2470	endif
2471	if (metcloud.and.(.not.fast)) then
2472	call writediagfi(ngrid,"zdtch4cloud","ch4 cloud","T s-1",&
2473	3,zdtch4cloud)
2474	call writediagfi(ngrid,"zdqch4cloud","ch4 cloud","T s-1",&
2475	3,zdqch4cloud(1,1,igcm_ch4_gas))
2476	endif
2477
2478	endif
2479
2480	if (carbox) then
2481	! call writediagfi(ngrid,"zdtcocloud","tendancy T cocloud","K",3,zdtcocloud)
2482	call writediagfi(ngrid,'co_iceflux','co_iceflux',&
2483	"kg m^-2 s^-1",2,flusurf(1,igcm_co_ice) )
2484	call writediagfi(ngrid,"vmr_co","vmr_co","%",2,vmr_co)
2485	if (.not.fast) THEN
2486	call writediagfi(ngrid,"zrho_co","zrho_co","kg.m-3",3,zrho_co(:,:))
2487	endif
2488	endif
2489
2490	if (haze) then
2491	! call writediagfi(ngrid,"zrho_haze","zrho_haze","kg.m-3",3,zrho_haze(:,:))
2492	call writediagfi(ngrid,"zdqrho_photprec","zdqrho_photprec",&
2493	"kg.m-3.s-1",3,zdqrho_photprec(:,:))
2494	call writediagfi(ngrid,"zdqphot_prec","zdqphot_prec","",&
2495	3,zdqphot_prec(:,:))
2496	call writediagfi(ngrid,"zdqhaze_ch4","zdqhaze_ch4","",&
2497	3,zdqhaze(:,:,igcm_ch4_gas))
2498	call writediagfi(ngrid,"zdqhaze_prec","zdqhaze_prec","",&
2499	3,zdqhaze(:,:,igcm_prec_haze))
2500	if (igcm_haze.ne.0) then
2501	call writediagfi(ngrid,"zdqhaze_haze","zdqhaze_haze","",&
2502	3,zdqhaze(:,:,igcm_haze))
2503	if (sedimentation) then
2504	call writediagfi(ngrid,"zdqssed_haze","zdqssed haze",&
2505	"kg/m2/s",2,zdqssed(:,igcm_haze))
2506	endif
2507	endif
2508	call writediagfi(ngrid,"zdqphot_ch4","zdqphot_ch4","",&
2509	3,zdqphot_ch4(:,:))
2510	call writediagfi(ngrid,"zdqconv_prec","zdqconv_prec","",&
2511	3,zdqconv_prec(:,:))
2512	! call writediagfi(ngrid,"zdqhaze_col","zdqhaze col","kg/m2/s",
2513	! & 2,zdqhaze_col(:))
2514	endif
2515
2516	if (aerohaze) then
2517	call writediagfi(ngrid,"tau_col",&
2518	"Total aerosol optical depth","opacity",2,tau_col)
2519	endif
2520
2521	endif ! end of 'tracer'
2522
2523
2524	! Output spectrum.
2525	if(specOLR.and.corrk)then
2526	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
2527	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
2528	call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W/m^2/cm^-1",3,GSR_nu)
2529	endif
2530
2531	#else
2532	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
2533	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
2534	comm_ALBEQ(1:ngrid)=albedo_equivalent(1:ngrid)
2535	if (.not.calldifv) comm_LATENT_HF(:)=0.0
2536
2537	if ((tracer).and.(generic_condensation)) then
2538
2539	! If you have set generic_condensation (and not water) and you have set several GCS
2540	! then the outputs given to WRF will be only the ones for the last generic tracer
2541	! (Because it is rewritten every tracer in the loop)
2542	! WRF can take only one moist tracer
2543
2544	do iq=1,nq
2545	call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
2546
2547	if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
2548
2549	reffrad_generic_zeros_for_wrf(:,:) = 1.
2550
2551	comm_CLOUDFRAC(1:ngrid,1:nlayer) = cloudfrac(1:ngrid,1:nlayer)
2552	comm_TOTCLOUDFRAC(1:ngrid) = totcloudfrac(1:ngrid) !??????
2553	comm_SURFRAIN(1:ngrid) = zdqsrain_generic(1:ngrid,iq)
2554	comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_gas)
2555	comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_generic_ice)
2556	! comm_H2OICE_REFF(1:ngrid,1:nlayer) = reffrad_generic_zeros_for_wrf(1:ngrid,1:nlayer) ! for now zeros !
2557	!comm_H2OICE_REFF(1:ngrid,1:nlayer) = 0*zdtrain_generic(1:ngrid,1:nlayer) ! for now zeros !
2558	comm_REEVAP(1:ngrid) = reevap_precip_generic(1:ngrid,iq)
2559	comm_DTRAIN(1:ngrid,1:nlayer) = zdtrain_generic(1:ngrid,1:nlayer)
2560	comm_DTLSC(1:ngrid,1:nlayer) = dt_generic_condensation(1:ngrid,1:nlayer)
2561	comm_RH(1:ngrid,1:nlayer) = RH_generic(1:ngrid,1:nlayer,iq)
2562
2563	endif
2564	end do ! do iq=1,nq loop on tracers
2565
2566	else
2567	comm_CLOUDFRAC(1:ngrid,1:nlayer)=0.
2568	comm_TOTCLOUDFRAC(1:ngrid)=0.
2569	comm_SURFRAIN(1:ngrid)=0.
2570	comm_DQVAP(1:ngrid,1:nlayer)=0.
2571	comm_DQICE(1:ngrid,1:nlayer)=0.
2572	! comm_H2OICE_REFF(1:ngrid,1:nlayer)=0.
2573	comm_REEVAP(1:ngrid)=0.
2574	comm_DTRAIN(1:ngrid,1:nlayer)=0.
2575	comm_DTLSC(1:ngrid,1:nlayer)=0.
2576	comm_RH(1:ngrid,1:nlayer)=0.
2577
2578	endif ! if water, if generic_condensation, else
2579
2580	comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid)
2581	comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid)
2582	comm_FLUXTOP_LW(1:ngrid)=fluxtop_lw(1:ngrid)
2583	comm_FLUXSURF_SW(1:ngrid)=fluxsurf_sw(1:ngrid)
2584	comm_FLUXSURF_LW(1:ngrid)=fluxsurf_lw(1:ngrid)
2585	comm_FLXGRD(1:ngrid)=fluxgrd(1:ngrid)
2586	sensibFlux(1:ngrid) = zflubid(1:ngrid) - capcal(1:ngrid)*zdtsdif(1:ngrid) !!! ????
2587	comm_HR_DYN(1:ngrid,1:nlayer) = zdtdyn(1:ngrid,1:nlayer)
2588	#endif
2589
2590	! XIOS outputs
2591	#ifdef CPP_XIOS
2592	! Send fields to XIOS: (NB these fields must also be defined as
2593	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
2594	CALL send_xios_field("ls",zls)
2595
2596	CALL send_xios_field("ps",ps)
2597	CALL send_xios_field("area",cell_area)
2598	CALL send_xios_field("p",pplay)
2599	CALL send_xios_field("temperature",zt)
2600	CALL send_xios_field("u",zu)
2601	CALL send_xios_field("v",zv)
2602	call send_xios_field("w",pw)
2603
2604	CALL send_xios_field("ISR",fluxtop_dn)
2605	CALL send_xios_field("OLR",fluxtop_lw)
2606	CALL send_xios_field("ASR",fluxabs_sw)
2607
2608	if (specOLR .and. corrk) then
2609	call send_xios_field("OLR3D",OLR_nu)
2610	call send_xios_field("IR_Bandwidth",DWNI)
2611	call send_xios_field("VI_Bandwidth",DWNV)
2612	call send_xios_field("OSR3D",OSR_nu)
2613	call send_xios_field("GSR3D",GSR_nu)
2614	endif
2615
2616	if (lastcall.and.is_omp_master) then
2617	write(,) "physiq: call xios_context_finalize"
2618	call xios_context_finalize
2619	endif
2620	#endif
2621
2622	if (check_physics_outputs) then
2623	! Check the validity of updated fields at the end of the physics step
2624	call check_physics_fields("end of physiq:", zt, zu, zv, pplev, zq)
2625	endif
2626
2627	icount=icount+1
2628
2629	end subroutine physiq
2630
2631	end module physiq_mod

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

Download in other formats:

Original Format