Context Navigation

physiq_mod.F90 @ 3241

Last change on this file since 3241 was 3241, checked in by afalco, 9 months ago
Pluto PCM: outputs for CH4, CO, etc AF
Property svn:executable set to ``*
File size: 98.4 KB

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

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

Download in other formats:

Original Format