Context Navigation

physiq_mod.F90 @ 3411

Last change on this file since 3411 was 3411, checked in by afalco, 3 months ago
Pluto PCM: Fixed some issue with no gcm. AF
Property svn:executable set to ``*
File size: 101.3 KB

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

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

Download in other formats:

Original Format