Context Navigation

physiq_mod.F90 @ 3329

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

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

Download in other formats:

Original Format