Context Navigation

physiq_mod.F90 @ 3374

Last change on this file since 3374 was 3374, checked in by afalco, 5 months ago
Pluto PCM: rice initialized to 0. AF
Property svn:executable set to ``*
File size: 101.1 KB

Line
1	module physiq_mod
2
3	implicit none
4
5	contains
6
7	subroutine physiq(ngrid,nlayer,nq, &
8	firstcall,lastcall, &
9	pday,ptime,ptimestep, &
10	pplev,pplay,pphi, &
11	pu,pv,pt,pq, &
12	flxw, &
13	pdu,pdv,pdt,pdq,pdpsrf)
14
15	!!
16	use write_field_phy, only: Writefield_phy
17	!!
18	use ioipsl_getin_p_mod, only: getin_p
19	use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind, corrkdir, banddir
20	use generic_cloud_common_h, only : epsi_generic, Psat_generic
21	use gases_h, only: gnom, gfrac
22	use radcommon_h, only: sigma, glat, grav, BWNV, WNOI, DWNI, DWNV, WNOV
23	use suaer_corrk_mod, only: suaer_corrk
24	use radii_mod, only: su_aer_radii,haze_reffrad_fix
25	use aerosol_mod, only: iaero_haze, i_haze, haze_prof
26	use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, &
27	dryness
28	use comdiurn_h, only: coslat, sinlat, coslon, sinlon
29	use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
30	use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
31	use geometry_mod, only: latitude, longitude, cell_area
32	USE comgeomfi_h, only: totarea, totarea_planet
33	USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
34	igcm_n2,igcm_ch4_gas,igcm_ch4_ice,igcm_haze,&
35	igcm_co_gas,igcm_co_ice,igcm_prec_haze,lw_n2,&
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	if (diurnal) then
867	ztim1=SIN(declin)
868	ztim2=COS(declin)COS(2.pi*(zday-.5))
869	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
870
871	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
872	ztim1,ztim2,ztim3,mu0,fract)
873	else if(diurnal .eqv. .false.) then
874
875	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad)
876	! WARNING: this function appears not to work in 1D
877
878	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
879	mu0 = cos(pi*szangle/180.0)
880	fract= 1/(4*mu0) ! AF24: from pluto.old
881	endif
882
883	endif
884
885
886	! Pluto albedo /IT changes depending on surface ices (only in 3D)
887	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
888	if (ngrid.ne.1) then
889
890	!! Specific to change albedo of N2 so that Psurf
891	!! converges toward 1.4 Pa in "1989" seasons for Triton
892	!! converges toward 1.1 Pa in "2015" seasons for Pluto
893	if (convergeps) then
894	if (triton) then
895	! 1989 declination
896	if (declin180./pi.gt.-46..and.declin180./pi.lt.-45. &
897	.and.zday.gt.saveday+1000. &
898	.and.declin.lt.savedeclin) then
899	call globalaverage2d(ngrid,pplev(:,1),globave)
900	if (globave.gt.1.5) then
901	adjust=adjust+0.005
902	else if (globave.lt.1.3) then
903	adjust=adjust-0.005
904	endif
905	saveday=zday
906	endif
907	else
908	! Pluto : 2015 declination current epoch
909	if (declin180./pi.gt.50.and.declin180./pi.lt.51. &
910	.and.zday.gt.saveday+10000. &
911	.and.declin.gt.savedeclin) then
912	call globalaverage2d(ngrid,pplev(:,1),globave)
913	if (globave.gt.1.2) then
914	adjust=adjust+0.005
915	else if (globave.lt.1.) then
916	adjust=adjust-0.005
917	endif
918	saveday=zday
919	endif
920	endif
921	end if
922	end if ! if ngrid ne 1
923
924	call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat, &
925	capcal,adjust,dist_star,albedo,emis,flusurfold,ptimestep, &
926	zls)
927	! do i=2,ngrid
928	! albedo(i,:) = albedo(1,:)
929	! enddo
930
931	if (callrad) then
932	if( mod(icount-1,iradia).eq.0.or.lastcall) then
933
934	! Eclipse incoming sunlight !AF24: removed
935
936	!! call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
937	!! call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
938
939
940	if (corrk) then
941
942	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
943	! II.a Call correlated-k radiative transfer scheme
944	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
945	if(kastprof)then
946	print*,'kastprof should not = true here'
947	call abort
948	endif
949	! if(water) then !AF24: removed
950
951	if(generic_condensation) then
952	do iq=1,nq
953
954	call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
955
956	if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
957
958	epsi_generic=constants_epsi_generic(iq)
959
960	muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_gas))
961	muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_gas))
962
963	endif
964	end do ! do iq=1,nq loop on tracers
965	! take into account generic condensable specie (GCS) effect on mean molecular weight
966
967	else
968	muvar(1:ngrid,1:nlayer+1)=mugaz
969	endif
970
971	! if(ok_slab_ocean) then !AF24: removed
972
973	! standard callcorrk
974	if (oldplutocorrk) then
975	call callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf, &
976	albedo(:,1),emis,mu0,pplev,pplay,pt, &
977	zzlay,tsurf,fract,dist_star,aerosol, &
978	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, &
979	fluxabs_sw,fluxtop_dn,reffrad,tau_col,ptime,pday, &
980	cloudfrac,totcloudfrac,.false., &
981	firstcall,lastcall)
982	albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
983	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+ &
984	fluxsurf_sw(1:ngrid)*(1.-albedo(1:ngrid,1))
985	else
986	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
987	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
988	zzlay,tsurf,fract,dist_star,aerosol,muvar, &
989	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, &
990	fluxsurfabs_sw,fluxtop_lw, &
991	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,GSR_nu, &
992	int_dtaui,int_dtauv, &
993	tau_col,cloudfrac,totcloudfrac, &
994	.false.,firstcall,lastcall)
995	! Radiative flux from the sky absorbed by the surface (W.m-2).
996	GSR=0.0
997	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+ &
998	fluxsurfabs_sw(1:ngrid)
999	endif ! oldplutocorrk
1000	!GG (feb2021): Option to "artificially" decrease the raditive time scale in
1001	!the deep atmosphere press > 0.1 bar. Suggested by MT.
1002	!! COEFF_RAD to be "tuned" to facilitate convergence of tendency
1003
1004	!coeff_rad=0. ! 0 values, it doesn't accelerate the convergence
1005	!coeff_rad=0.5
1006	!coeff_rad=1.
1007	!do l=1, nlayer
1008	! do ig=1,ngrid
1009	! if(pplay(ig,l).ge.1.d4) then
1010	! zdtsw(ig,l)=zdtsw(ig,l)(pplay(ig,l)/1.d4)*coeff_rad
1011	! zdtlw(ig,l)=zdtlw(ig,l)(pplay(ig,l)/1.d4)*coeff_rad
1012	! endif
1013	! enddo
1014	!enddo
1015
1016	! AF24: removed CLFvarying Option
1017
1018
1019	!if(noradsurf)then ! no lower surface; SW flux just disappears
1020	! GSR = SUM(fluxsurf_sw(1:ngrid)*cell_area(1:ngrid))/totarea
1021	! fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
1022	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
1023	!endif
1024
1025	! Net atmospheric radiative heating rate (K.s-1)
1026	dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
1027
1028	! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that !
1029	if (firstcall .and. albedo_spectral_mode) then
1030	call spectral_albedo_calc(albedo_snow_SPECTV)
1031	endif
1032
1033	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1034	! II.b Call Newtonian cooling scheme !AF24: removed
1035	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1036	else
1037	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1038	! II.! Atmosphere has no radiative effect
1039	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1040	fluxtop_dn(1:ngrid) = fract(1:ngrid)mu0(1:ngrid)Fat1AU/dist_star**2
1041	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
1042	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
1043	endif
1044	fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
1045	print*,'------------WARNING---WARNING------------' ! by MT2015.
1046	print*,'You are in corrk=false mode, '
1047	print*,'and the surface albedo is taken equal to the first visible spectral value'
1048
1049	albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
1050	fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
1051	fluxabs_sw(1:ngrid)=fluxsurfabs_sw(1:ngrid)
1052	fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid)
1053	fluxtop_lw(1:ngrid) = emis(1:ngrid)sigmatsurf(1:ngrid)**4
1054
1055	dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating
1056
1057	endif ! end of corrk
1058
1059	endif ! of if(mod(icount-1,iradia).eq.0)
1060
1061
1062	! Transformation of the radiative tendencies
1063	! ------------------------------------------
1064	zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
1065	zplanck(1:ngrid)=emis(1:ngrid)sigmazplanck(1:ngrid)*zplanck(1:ngrid)
1066	fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
1067	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
1068
1069	! Test of energy conservation
1070	!----------------------------
1071	if(enertest)then
1072	call planetwide_sumval(cppmassarea(:,:)zdtsw(:,:)/totarea_planet,dEtotSW)
1073	call planetwide_sumval(cppmassarea(:,:)zdtlw(:,:)/totarea_planet,dEtotLW)
1074	!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
1075	call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
1076	call planetwide_sumval((fluxsurf_lw(:)emis(:)-zplanck(:))cell_area(:)/totarea_planet,dEtotsLW)
1077	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
1078	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
1079	if (is_master) then
1080	print*,'---------------------------------------------------------------'
1081	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
1082	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
1083	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
1084	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
1085	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
1086	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
1087	endif
1088	endif ! end of 'enertest'
1089
1090	endif ! of if (callrad)
1091
1092	!! call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1093	!! call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1094
1095
1096	! --------------------------------------------
1097	! III. Vertical diffusion (turbulent mixing) :
1098	! --------------------------------------------
1099
1100	if (calldifv) then
1101
1102	zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
1103
1104	if (oldplutovdifc) then
1105	zdum1(:,:) = 0
1106	zdum2(:,:) = 0
1107	zdh(:,:)=pdt(:,:)/zpopsk(:,:)
1108
1109	! Calling vdif (Martian version WITH N2 condensation)
1110	CALL vdifc_pluto(ngrid,nlayer,nq,zpopsk, &
1111	ptimestep,capcal,lwrite, &
1112	pplay,pplev,zzlay,zzlev,z0, &
1113	pu,pv,zh,pq,pt,tsurf,emis,qsurf, &
1114	zdum1,zdum2,zdh,pdq,pdt,zflubid, &
1115	zdudif,zdvdif,zdhdif,zdtsdif,q2, &
1116	zdqdif,zdqsdif,qsat_ch4,qsat_ch4_l1) !,zq1temp_ch4,qsat_ch4)
1117
1118	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1119
1120	bcond=1./tcond1p4Pa
1121	acond=r/lw_n2
1122
1123	!------------------------------------------------------------------
1124	! test methane conservation
1125	! if(methane)then
1126	! call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
1127	! & ptimestep,pplev,zdqdif,zdqsdif,'CH4',' vdifc ')
1128	! endif ! methane
1129	!------------------------------------------------------------------
1130	! test CO conservation
1131	! if(carbox)then
1132	! call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
1133	! & ptimestep,pplev,zdqdif,zdqsdif,'CO ',' vdifc ')
1134	! endif ! carbox
1135	!------------------------------------------------------------------
1136
1137	! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
1138	else if (UseTurbDiff) then
1139
1140	call turbdiff(ngrid,nlayer,nq, &
1141	ptimestep,capcal, &
1142	pplay,pplev,zzlay,zzlev,z0, &
1143	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
1144	pdt,pdq,zflubid, &
1145	zdudif,zdvdif,zdtdif,zdtsdif, &
1146	sensibFlux,q2,zdqdif,zdqevap,zdqsdif, &
1147	taux,tauy)
1148
1149	else ! if (oldplutovdifc) .and. (UseTurbDiff)
1150
1151	zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1152
1153	call vdifc(ngrid,nlayer,nq,zpopsk, &
1154	ptimestep,capcal,lwrite, &
1155	pplay,pplev,zzlay,zzlev,z0, &
1156	pu,pv,zh,pq,tsurf,emis,qsurf, &
1157	zdh,pdq,zflubid, &
1158	zdudif,zdvdif,zdhdif,zdtsdif, &
1159	sensibFlux,q2,zdqdif,zdqsdif)
1160
1161	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1162	zdqevap(1:ngrid,1:nlayer)=0.
1163
1164	end if !end of 'UseTurbDiff'
1165
1166	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
1167
1168	!!! this is always done, except for turbulence-resolving simulations
1169	if (.not. turb_resolved) then
1170	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
1171	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
1172	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
1173	endif
1174
1175	! if(ok_slab_ocean)then !AF24: removed
1176	! flux_sens_lat(1:ngrid)=(zdtsdif(1:ngrid)*capcal(1:ngrid)-fluxrad(1:ngrid))
1177	! endif
1178
1179	!! call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_gas))
1180
1181	if (tracer) then
1182	pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
1183	dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
1184	end if ! of if (tracer)
1185
1186	!! call writediagfi(ngrid,"vdifc_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1187	!! call writediagfi(ngrid,"vdifc_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1188
1189	! test energy conservation
1190	!-------------------------
1191	if(enertest)then
1192
1193	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
1194	do ig = 1, ngrid
1195	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
1196	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
1197	enddo
1198
1199	call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
1200	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
1201	call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
1202	call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
1203
1204	if (is_master) then
1205
1206	if (UseTurbDiff) then
1207	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1208	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
1209	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1210	else
1211	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1212	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
1213	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1214	end if
1215	endif ! end of 'is_master'
1216
1217	! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
1218	endif ! end of 'enertest'
1219
1220
1221	! ! Test water conservation. !AF24: removed
1222
1223	else ! calldifv
1224
1225	! if(.not.newtonian)then
1226	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
1227
1228	endif ! end of 'calldifv'
1229
1230
1231	!-------------------
1232	! IV. Convection :
1233	!-------------------
1234
1235	! ~~~~~~~~~~~~~~~~~~~~~~~~~~
1236	! IV.a Thermal plume model : AF24: removed
1237	! ~~~~~~~~~~~~~~~~~~~~~~~~~~
1238
1239	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1240	! IV.b Dry convective adjustment :
1241	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1242
1243	if(calladj) then
1244
1245	zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1246	zduadj(1:ngrid,1:nlayer)=0.0
1247	zdvadj(1:ngrid,1:nlayer)=0.0
1248	zdhadj(1:ngrid,1:nlayer)=0.0
1249
1250
1251	call convadj(ngrid,nlayer,nq,ptimestep, &
1252	pplay,pplev,zpopsk, &
1253	pu,pv,zh,pq, &
1254	pdu,pdv,zdh,pdq, &
1255	zduadj,zdvadj,zdhadj, &
1256	zdqadj)
1257
1258	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
1259	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer)
1260	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
1261	zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1262
1263	if(tracer) then
1264	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
1265	end if
1266
1267	! Test energy conservation
1268	if(enertest)then
1269	call planetwide_sumval(cppmassarea(:,:)zdtadj(:,:)/totarea_planet,dEtot)
1270	if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
1271	endif
1272
1273	! ! Test water conservation !AF24: removed
1274
1275	endif ! end of 'calladj'
1276	!----------------------------------------------
1277	! Non orographic Gravity Waves: AF24: removed
1278	!---------------------------------------------
1279
1280	!-----------------------------------------------
1281	! V. Nitrogen condensation-sublimation :
1282	!-----------------------------------------------
1283
1284	if (n2cond) then
1285
1286	if (.not.tracer) then
1287	print*,'We need a N2 ice tracer to condense N2'
1288	call abort
1289	endif
1290
1291	call condense_n2(ngrid,nlayer,nq,ptimestep, &
1292	capcal,pplay,pplev,tsurf,pt, &
1293	pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, &
1294	qsurf(1,igcm_n2),albedo,emis, &
1295	zdtc,zdtsurfc,pdpsrf,zduc,zdvc, &
1296	zdqc,zdqsc(1,igcm_n2))
1297
1298	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
1299	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer)+zdvc(1:ngrid,1:nlayer)
1300	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer)+zduc(1:ngrid,1:nlayer)
1301	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
1302
1303	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq)
1304	dqsurf(1:ngrid,igcm_n2) = dqsurf(1:ngrid,igcm_n2) + zdqsc(1:ngrid,igcm_n2)
1305
1306	!! call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1307	!! call writediagfi(ngrid,"condense_n2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1308
1309	! test energy conservation
1310	if(enertest)then
1311	call planetwide_sumval(cppmassarea(:,:)zdtc(:,:)/totarea_planet,dEtot)
1312	call planetwide_sumval(capcal(:)zdtsurfc(:)cell_area(:)/totarea_planet,dEtots)
1313	if (is_master) then
1314	print*,'In n2cloud atmospheric energy change =',dEtot,' W m-2'
1315	print*,'In n2cloud surface energy change =',dEtots,' W m-2'
1316	endif
1317	endif
1318
1319	endif ! end of 'n2cond'
1320
1321
1322	!---------------------------------------------
1323	! VI. Specific parameterizations for tracers
1324	!---------------------------------------------
1325
1326	if (tracer) then
1327
1328
1329
1330	! 7a. Methane, CO, and ice
1331	! ---------------------------------------
1332	! Methane ice condensation in the atmosphere
1333	! ----------------------------------------
1334	rice_ch4(:,:)=0 ! initialization needed for callsedim
1335	zdqch4cloud(:,:,:)=0.
1336	if ((methane).and.(metcloud).and.(.not.fast)) THEN
1337	call ch4cloud(ngrid,nlayer,naerkind,ptimestep, &
1338	pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, &
1339	pq,pdq,zdqch4cloud,zdqsch4cloud,zdtch4cloud, &
1340	nq,rice_ch4)
1341
1342	DO l=1,nlayer
1343	DO ig=1,ngrid
1344	pdq(ig,l,igcm_ch4_gas)=pdq(ig,l,igcm_ch4_gas)+ &
1345	zdqch4cloud(ig,l,igcm_ch4_gas)
1346	pdq(ig,l,igcm_ch4_ice)=pdq(ig,l,igcm_ch4_ice)+ &
1347	zdqch4cloud(ig,l,igcm_ch4_ice)
1348	ENDDO
1349	ENDDO
1350
1351	! Increment methane ice surface tracer tendency
1352	DO ig=1,ngrid
1353	dqsurf(ig,igcm_ch4_ice)=dqsurf(ig,igcm_ch4_ice)+ &
1354	zdqsch4cloud(ig,igcm_ch4_ice)
1355	ENDDO
1356
1357	! update temperature tendancy
1358	DO ig=1,ngrid
1359	DO l=1,nlayer
1360	pdt(ig,l)=pdt(ig,l)+zdtch4cloud(ig,l)
1361	ENDDO
1362	ENDDO
1363	end if
1364
1365	! ---------------------------------------
1366	! CO ice condensation in the atmosphere
1367	! ----------------------------------------
1368	zdqcocloud(:,:,:)=0.
1369	IF ((carbox).and.(monoxcloud).and.(.not.fast)) THEN
1370	call cocloud(ngrid,nlayer,naerkind,ptimestep, &
1371	pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, &
1372	pq,pdq,zdqcocloud,zdqscocloud,zdtcocloud, &
1373	nq,rice_co,qsurf(1,igcm_n2),dqsurf(1,igcm_n2))
1374
1375	DO l=1,nlayer
1376	DO ig=1,ngrid
1377	pdq(ig,l,igcm_co_gas)=pdq(ig,l,igcm_co_gas)+ &
1378	zdqcocloud(ig,l,igcm_co_gas)
1379	pdq(ig,l,igcm_co_ice)=pdq(ig,l,igcm_co_ice)+ &
1380	zdqcocloud(ig,l,igcm_co_ice)
1381	ENDDO
1382	ENDDO
1383
1384	! Increment CO ice surface tracer tendency
1385	DO ig=1,ngrid
1386	dqsurf(ig,igcm_co_ice)=dqsurf(ig,igcm_co_ice)+ &
1387	zdqscocloud(ig,igcm_co_ice)
1388	ENDDO
1389
1390	! update temperature tendancy
1391	DO ig=1,ngrid
1392	DO l=1,nlayer
1393	pdt(ig,l)=pdt(ig,l)+zdtcocloud(ig,l)
1394	ENDDO
1395	ENDDO
1396	ELSE
1397	rice_co(:,:)=0 ! initialization needed for callsedim
1398	END IF ! of IF (carbox)
1399
1400	!------------------------------------------------------------------
1401	! test methane conservation
1402	! if(methane)then
1403	! call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
1404	! & ptimestep,pplev,zdqch4cloud,zdqsch4cloud,'CH4','ch4clou')
1405	! endif ! methane
1406	!------------------------------------------------------------------
1407	! test CO conservation
1408	! if(carbox)then
1409	! call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
1410	! & ptimestep,pplev,zdqcocloud,zdqscocloud,'CO ','cocloud')
1411	! endif ! carbox
1412	!------------------------------------------------------------------
1413
1414	! 7b. Haze particle production
1415	! -------------------
1416	IF (haze) THEN
1417
1418	zdqphot_prec(:,:)=0.
1419	zdqphot_ch4(:,:)=0.
1420	zdqhaze(:,:,:)=0
1421	! Forcing to a fixed haze profile if haze_proffix
1422	if (haze_proffix.and.i_haze.gt.0.) then
1423	call haze_prof(ngrid,nlayer,zzlay,pplay,pt, &
1424	reffrad,profmmr)
1425	zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze)) &
1426	/ptimestep
1427	else
1428	call hazecloud(ngrid,nlayer,nq,ptimestep, &
1429	pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,zdqhaze, &
1430	zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
1431	endif
1432
1433	DO iq=1, nq ! should be updated
1434	DO l=1,nlayer
1435	DO ig=1,ngrid
1436	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqhaze(ig,l,iq)
1437	ENDDO
1438	ENDDO
1439	ENDDO
1440
1441	ENDIF
1442
1443	IF (fast.and.fasthaze) THEN
1444	call prodhaze(ngrid,nlayer,nq,ptimestep,pplev,pq,pdq,dist_star, &
1445	mu0,declin,zdqprodhaze,zdqsprodhaze,gradflux,fluxbot, &
1446	fluxlym_sol_bot,fluxlym_ipm_bot,flym_sol,flym_ipm)
1447
1448	DO ig=1,ngrid
1449	pdq(ig,1,igcm_ch4_gas)=pdq(ig,1,igcm_ch4_gas)+ &
1450	zdqprodhaze(ig,igcm_ch4_gas)
1451	pdq(ig,1,igcm_prec_haze)=pdq(ig,1,igcm_prec_haze)+ &
1452	zdqprodhaze(ig,igcm_prec_haze)
1453	pdq(ig,1,igcm_haze)=abs(pdq(ig,1,igcm_haze)+ &
1454	zdqprodhaze(ig,igcm_haze))
1455	qsurf(ig,igcm_haze)= qsurf(ig,igcm_haze)+ &
1456	zdqsprodhaze(ig)*ptimestep
1457	ENDDO
1458
1459	ENDIF
1460
1461
1462
1463	! -------------------------
1464	! VI.3. Aerosol particles
1465	! -------------------------
1466
1467	!Generic Condensation
1468	if (generic_condensation) then
1469	call condensation_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay, &
1470	pt,pq,pdt,pdq,dt_generic_condensation, &
1471	dqvaplscale_generic,dqcldlscale_generic,rneb_generic)
1472	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dt_generic_condensation(1:ngrid,1:nlayer)
1473	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+dqvaplscale_generic(1:ngrid,1:nlayer,1:nq)
1474	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+dqcldlscale_generic(1:ngrid,1:nlayer,1:nq)
1475
1476	if(enertest)then
1477	do ig=1,ngrid
1478	genericconddE(ig) = cppSUM(mass(:,:)dt_generic_condensation(:,:))
1479	enddo
1480
1481	call planetwide_sumval(cppmassarea(:,:)dt_generic_condensation(:,:)/totarea_planet,dEtot)
1482
1483	if (is_master) print*,'In generic condensation atmospheric energy change =',dEtot,' W m-2'
1484	end if
1485
1486	! Let's loop on tracers
1487	cloudfrac(:,:)=0.0
1488	do iq=1,nq
1489	call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
1490	if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
1491	do l = 1, nlayer
1492	do ig=1,ngrid
1493	cloudfrac(ig,l)=rneb_generic(ig,l,iq)
1494	enddo
1495	enddo
1496	endif
1497	end do ! do iq=1,nq loop on tracers
1498	! endif ! .not. water
1499
1500	endif !generic_condensation
1501
1502	! Sedimentation.
1503	if (sedimentation) then
1504
1505	zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
1506	zdqssed(1:ngrid,1:nq) = 0.0
1507
1508	if (oldplutosedim)then
1509	call callsedim_pluto(ngrid,nlayer,ptimestep, &
1510	pplev,zzlev,pt,pdt,rice_ch4,rice_co, &
1511	pq,pdq,zdqsed,zdqssed,nq,pphi)
1512
1513	else
1514	call callsedim(ngrid,nlayer,ptimestep, &
1515	pplev,zzlev,pt,pdt,pq,pdq, &
1516	zdqsed,zdqssed,nq,pphi)
1517	endif
1518
1519	! Whether it falls as rain or snow depends only on the surface temperature
1520	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
1521	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
1522
1523	end if ! end of 'sedimentation'
1524
1525	! ---------------
1526	! VI.4. Updates
1527	! ---------------
1528
1529	! Updating Atmospheric Mass and Tracers budgets.
1530	if(mass_redistrib) then
1531
1532	zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * 0
1533	! ( zdqevap(1:ngrid,1:nlayer) &
1534	! ! + zdqrain(1:ngrid,1:nlayer,igcm_h2o_gas) &
1535	! ! + dqmoist(1:ngrid,1:nlayer,igcm_h2o_gas) &
1536	! + dqvaplscale(1:ngrid,1:nlayer) )
1537
1538	do ig = 1, ngrid
1539	zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
1540	enddo
1541
1542	! call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
1543	! call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
1544	call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
1545
1546	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
1547	capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
1548	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
1549	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
1550
1551	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
1552	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
1553	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
1554	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
1555	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
1556	pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
1557	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
1558
1559	endif
1560
1561	! call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1562
1563	!! call writediagfi(ngrid,"slab_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1564	!! call writediagfi(ngrid,"slab_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1565
1566
1567	! ------------------
1568	! VI.5. Slab Ocean !AF24: removed
1569	! ------------------
1570
1571	! -----------------------------
1572	! VI.6. Surface Tracer Update
1573	! -----------------------------
1574
1575	! AF24: deleted hydrology
1576
1577	qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
1578
1579	! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
1580	! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
1581	qsurf_hist(:,:) = qsurf(:,:)
1582
1583	! if(ice_update)then
1584	! ice_min(1:ngrid)=min(ice_min(1:ngrid),qsurf(1:ngrid,igcm_h2o_ice))
1585	! endif
1586
1587	endif! end of if 'tracer'
1588
1589	if (conservn2) then
1590	write(,) 'conservn2 tracer'
1591	! call testconservmass(ngrid,nlayer,pplev(:,1)+ &
1592	! pdpsrf(:)*ptimestep,qsurf(:,1))
1593	endif
1594
1595	DO ig=1,ngrid
1596	flusurf(ig,igcm_n2)=(qsurf(ig,igcm_n2)- &
1597	qsurf1(ig,igcm_n2))/ptimestep
1598	flusurfold(ig,igcm_n2)=flusurf(ig,igcm_n2)
1599	if (methane) then
1600	flusurf(ig,igcm_ch4_ice)=(qsurf(ig,igcm_ch4_ice)- &
1601	qsurf1(ig,igcm_ch4_ice))/ptimestep
1602	flusurfold(ig,igcm_ch4_ice)=flusurf(ig,igcm_ch4_ice)
1603	endif
1604	if (carbox) then
1605	flusurf(ig,igcm_co_ice)=(qsurf(ig,igcm_co_ice)- &
1606	qsurf1(ig,igcm_co_ice))/ptimestep
1607	!flusurfold(ig,igcm_co_ice)=flusurf(ig,igcm_co_ice)
1608	endif
1609	ENDDO
1610
1611	!! Special source of haze particle !
1612	! todo: should be placed in haze and use tendency of n2 instead of flusurf
1613	IF (source_haze) THEN
1614	write(,) "Source haze not supported yet."
1615	stop
1616	! call hazesource(ngrid,nlayer,nq,ptimestep, &
1617	! pplev,flusurf,mu0,zdq_source)
1618
1619	DO iq=1, nq
1620	DO l=1,nlayer
1621	DO ig=1,ngrid
1622	pdq(ig,l,iq)=pdq(ig,l,iq)+zdq_source(ig,l,iq)
1623	ENDDO
1624	ENDDO
1625	ENDDO
1626	ENDIF
1627
1628
1629	!------------------------------------------------
1630	! VII. Surface and sub-surface soil temperature
1631	!------------------------------------------------
1632
1633	! For diagnostic
1634	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1635	DO ig=1,ngrid
1636	rho(ig,1) = pplay(ig,1)/(r*pt(ig,1))
1637	sensiblehf1(ig)=rho(ig,1)cpp(0.4/log(zzlay(ig,1)/z0))*2 &
1638	(pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1))*0.5 &
1639	(tsurf(ig)-pt(ig,1))
1640	sensiblehf2(ig)=zflubid(ig)-capcal(ig)*zdtsdif(ig)
1641
1642	ENDDO
1643
1644
1645	! ! Increment surface temperature
1646	! if(ok_slab_ocean)then !AF24: removed
1647
1648	tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
1649	! Compute soil temperatures and subsurface heat flux.
1650	if (callsoil) then
1651	call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, &
1652	ptimestep,tsurf,tsoil,capcal,fluxgrd)
1653	endif
1654
1655
1656	! if (ok_slab_ocean) then !AF24: removed
1657
1658	! Test energy conservation
1659	if(enertest)then
1660	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf(:)/totarea_planet,dEtots)
1661	if (is_master) print*,'Surface energy change =',dEtots,' W m-2'
1662	endif
1663
1664
1665	!---------------------------------------------------
1666	! VIII. Perform diagnostics and write output files
1667	!---------------------------------------------------
1668
1669	! Note : For output only: the actual model integration is performed in the dynamics.
1670
1671
1672	! Temperature, zonal and meridional winds.
1673	zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
1674	zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
1675	zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
1676
1677	!! Recast thermal plume vertical velocity array for outputs
1678	!! AF24: removed
1679
1680	! Diagnostic.
1681	zdtdyn(1:ngrid,1:nlayer) = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep
1682	ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
1683
1684	zdudyn(1:ngrid,1:nlayer) = (pu(1:ngrid,1:nlayer)-zuprevious(1:ngrid,1:nlayer)) / ptimestep
1685	zuprevious(1:ngrid,1:nlayer) = zu(1:ngrid,1:nlayer)
1686
1687	if(firstcall)then
1688	zdtdyn(1:ngrid,1:nlayer)=0.0
1689	zdudyn(1:ngrid,1:nlayer)=0.0
1690	endif
1691
1692	! Dynamical heating diagnostic.
1693	do ig=1,ngrid
1694	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp
1695	enddo
1696
1697	! Tracers.
1698	zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
1699
1700	! Surface pressure.
1701	ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
1702
1703	! pressure density !pluto specific
1704	IF (.not.fast) then !
1705	do ig=1,ngrid
1706	do l=1,nlayer
1707	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1708	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1709	rho(ig,l) = zplay(ig,l)/(r*zt(ig,l))
1710	enddo
1711	zplev(ig,nlayer+1)=pplev(ig,nlayer+1)/pplev(ig,1)*ps(ig)
1712	enddo
1713	ENDIF
1714
1715
1716	! Surface and soil temperature information
1717	call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
1718	call planetwide_minval(tsurf(:),Ts2)
1719	call planetwide_maxval(tsurf(:),Ts3)
1720	if(callsoil)then
1721	TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
1722	if (is_master) then
1723	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
1724	print*,Ts1,Ts2,Ts3,TsS
1725	end if
1726	else
1727	if (is_master) then
1728	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
1729	print*,Ts1,Ts2,Ts3
1730	endif
1731	end if
1732
1733
1734	! Check the energy balance of the simulation during the run
1735	if(corrk)then
1736
1737	call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
1738	call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
1739	call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
1740	call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
1741	call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
1742	do ig=1,ngrid
1743	if(fluxtop_dn(ig).lt.0.0)then
1744	print*,'fluxtop_dn has gone crazy'
1745	print*,'fluxtop_dn=',fluxtop_dn(ig)
1746	print*,'tau_col=',tau_col(ig)
1747	print*,'aerosol=',aerosol(ig,:,:)
1748	print*,'temp= ',pt(ig,:)
1749	print*,'pplay= ',pplay(ig,:)
1750	call abort
1751	endif
1752	end do
1753
1754	if(ngrid.eq.1)then
1755	DYN=0.0
1756	endif
1757
1758	if (is_master) then
1759	print*,' ISR ASR OLR GND DYN [W m^-2]'
1760	print*, ISR,ASR,OLR,GND,DYN
1761	endif
1762
1763	if(enertest .and. is_master)then
1764	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
1765	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
1766	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
1767	endif
1768
1769	if(meanOLR .and. is_master)then
1770	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
1771	! to record global radiative balance
1772	open(92,file="rad_bal.out",form='formatted',position='append')
1773	write(92,*) zday,ISR,ASR,OLR
1774	close(92)
1775	open(93,file="tem_bal.out",form='formatted',position='append')
1776	if(callsoil)then
1777	write(93,*) zday,Ts1,Ts2,Ts3,TsS
1778	else
1779	write(93,*) zday,Ts1,Ts2,Ts3
1780	endif
1781	close(93)
1782	endif
1783	endif
1784
1785	endif ! end of 'corrk'
1786
1787
1788	! Diagnostic to test radiative-convective timescales in code.
1789	if(testradtimes)then
1790	open(38,file="tau_phys.out",form='formatted',position='append')
1791	ig=1
1792	do l=1,nlayer
1793	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1794	enddo
1795	close(38)
1796	print*,'As testradtimes enabled,'
1797	print*,'exiting physics on first call'
1798	call abort
1799	endif
1800
1801
1802	! Compute column amounts (kg m-2) if tracers are enabled.
1803	if(tracer)then
1804	qcol(1:ngrid,1:nq)=0.0
1805	do iq=1,nq
1806	do ig=1,ngrid
1807	qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
1808	enddo
1809	enddo
1810
1811	endif ! end of 'tracer'
1812
1813
1814	! ! Test for water conservation. !AF24: removed
1815
1816	! Calculate RH_generic (Generic Relative Humidity) for diagnostic.
1817	if(generic_condensation)then
1818	RH_generic(:,:,:)=0.0
1819	do iq=1,nq
1820
1821	call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
1822
1823	if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
1824
1825	do l = 1, nlayer
1826	do ig=1,ngrid
1827	call Psat_generic(zt(ig,l),pplay(ig,l),metallicity,psat_tmp_generic,qsat_generic(ig,l,iq))
1828	RH_generic(ig,l,iq) = zq(ig,l,igcm_generic_gas) / qsat_generic(ig,l,iq)
1829	enddo
1830	enddo
1831
1832	end if
1833
1834	end do ! iq=1,nq
1835
1836	endif ! end of 'generic_condensation'
1837
1838
1839	if (methane) then
1840	IF (fast) then ! zq is the mixing ratio supposingly mixed in all atmosphere
1841	DO ig=1,ngrid
1842	vmr_ch4(ig)=zq(ig,1,igcm_ch4_gas)* &
1843	mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.
1844	ENDDO
1845	ELSE
1846	DO ig=1,ngrid
1847	! compute vmr methane
1848	vmr_ch4(ig)=qcol(ig,igcm_ch4_gas)* &
1849	g/ps(ig)mmol(igcm_n2)/mmol(igcm_ch4_gas)100.
1850	! compute density methane
1851	DO l=1,nlayer
1852	zrho_ch4(ig,l)=zq(ig,l,igcm_ch4_gas)*rho(ig,l)
1853	ENDDO
1854	ENDDO
1855	ENDIF
1856	endif
1857
1858	if (carbox) then
1859	IF (fast) then
1860	DO ig=1,ngrid
1861	vmr_co(ig)=zq(ig,1,igcm_co_gas)* &
1862	mmol(igcm_n2)/mmol(igcm_co_gas)*100.
1863	ENDDO
1864	ELSE
1865	DO ig=1,ngrid
1866	! compute vmr CO
1867	vmr_co(ig)=qcol(ig,igcm_co_gas)* &
1868	g/ps(ig)mmol(igcm_n2)/mmol(igcm_co_gas)100.
1869	! compute density CO
1870	DO l=1,nlayer
1871	zrho_co(ig,l)=zq(ig,l,igcm_co_gas)*rho(ig,l)
1872	ENDDO
1873	ENDDO
1874	ENDIF
1875	endif
1876
1877	zrho_haze(:,:)=0.
1878	zdqrho_photprec(:,:)=0.
1879	IF (haze.and.aerohaze) then
1880	DO ig=1,ngrid
1881	DO l=1,nlayer
1882	zrho_haze(ig,l)=zq(ig,l,igcm_haze)*rho(ig,l)
1883	zdqrho_photprec(ig,l)=zdqphot_prec(ig,l)*rho(ig,l)
1884	ENDDO
1885	ENDDO
1886	ENDIF
1887
1888	IF (fasthaze) then
1889	DO ig=1,ngrid
1890	qcol(ig,igcm_haze)=zq(ig,1,igcm_haze)*pplev(ig,1)/g
1891	qcol(ig,igcm_prec_haze)=zq(ig,1,igcm_prec_haze)*pplev(ig,1)/g
1892	ENDDO
1893	ENDIF
1894
1895	! Info about Ls, declin...
1896	IF (fast) THEN
1897	if (is_master) write(,),'Ls=',zls180./pi,' dec=',declin180./pi
1898	if (is_master) write(,),'zday=',zday,' ps=',globave
1899	IF (lastcall) then
1900	if (is_master) write(,),'lastcall'
1901	ENDIF
1902	ELSE
1903	if (is_master) write(,),'Ls=',zls180./pi,'decli=',declin180./pi,'zday=',zday
1904	ENDIF
1905
1906	lecttsoil=0 ! default value for lecttsoil
1907	call getin_p("lecttsoil",lecttsoil)
1908	IF (lastcall.and.(ngrid.EQ.1).and.(lecttsoil.eq.1)) THEN
1909	! save tsoil temperature profile for 1D profile
1910	OPEN(13,file='proftsoil.out',form='formatted')
1911	DO i=1,nsoilmx
1912	write(13,*) tsoil(1,i)
1913	ENDDO
1914	CLOSE(13)
1915	ENDIF
1916
1917
1918	if (is_master) print,'--> Ls =',zls180./pi
1919
1920
1921	!----------------------------------------------------------------------
1922	! Writing NetCDF file "RESTARTFI" at the end of the run
1923	!----------------------------------------------------------------------
1924
1925	! Note: 'restartfi' is stored just before dynamics are stored
1926	! in 'restart'. Between now and the writting of 'restart',
1927	! there will have been the itau=itau+1 instruction and
1928	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1929	! thus we store for time=time+dtvr
1930
1931
1932
1933	if(lastcall) then
1934	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1935
1936	!! Update surface ice distribution to iterate to steady state if requested
1937	!! AF24: removed
1938
1939	! endif
1940	if (paleo) then
1941	! time range for tendencies of ice flux qsurfyear
1942	zdt_tot=year_day ! Last year of simulation
1943
1944	masslost(:)=0.
1945	massacc(:)=0.
1946
1947	DO ig=1,ngrid
1948	! update new reservoir of ice on the surface
1949	DO iq=1,nq
1950	! kg/m2 to be sublimed or condensed during paleoyears
1951	qsurfyear(ig,iq)=qsurfyear(ig,iq)* &
1952	paleoyears365.25/(zdt_totdaysec/86400.)
1953
1954	! special case if we sublime the entire reservoir
1955	!! AF: TODO : fix following lines (real_area), using line below:
1956	! call planetwide_sumval((-qsurfyear(:,iq)-qsurf(:,iq))*cell_area(:),masslost)
1957
1958	! IF (-qsurfyear(ig,iq).gt.qsurf(ig,iq)) THEN
1959	! masslost(iq)=masslost(iq)+real_area(ig)* &
1960	! (-qsurfyear(ig,iq)-qsurf(ig,iq))
1961	! qsurfyear(ig,iq)=-qsurf(ig,iq)
1962	! ENDIF
1963
1964	! IF (qsurfyear(ig,iq).gt.0.) THEN
1965	! massacc(iq)=massacc(iq)+real_area(ig)*qsurfyear(ig,iq)
1966	! ENDIF
1967
1968
1969	ENDDO
1970	ENDDO
1971
1972	DO ig=1,ngrid
1973	DO iq=1,nq
1974	qsurfpal(ig,iq)=qsurf(ig,iq)+qsurfyear(ig,iq)
1975	IF (qsurfyear(ig,iq).gt.0.) THEN
1976	qsurfpal(ig,iq)=qsurfpal(ig,iq)- &
1977	qsurfyear(ig,iq)*masslost(iq)/massacc(iq)
1978	ENDIF
1979	ENDDO
1980	ENDDO
1981	! Finally ensure conservation of qsurf
1982	DO iq=1,nq
1983	call globalaverage2d(ngrid,qsurf(:,iq),globaveice(iq))
1984	call globalaverage2d(ngrid,qsurfpal(:,iq), &
1985	globavenewice(iq))
1986	IF (globavenewice(iq).gt.0.) THEN
1987	qsurfpal(:,iq)=qsurfpal(:,iq)* &
1988	globaveice(iq)/globavenewice(iq)
1989	ENDIF
1990	ENDDO
1991
1992	! update new geopotential depending on the ice reservoir
1993	phisfipal(:)=phisfinew(:)+qsurfpal(:,igcm_n2)*g/1000.
1994	!phisfipal(ig)=phisfi(ig)
1995
1996
1997	if (kbo.or.triton) then ! case of Triton : we do not change the orbital parameters
1998
1999	pdaypal=pday ! no increment of pdaypal to keep same evolution of the subsolar point
2000	eccpal=1.-periastr/((periastr+apoastr)/2.) !no change of ecc
2001	peri_daypal=peri_day ! no change
2002	oblipal=obliquit ! no change
2003	tpalnew=tpal
2004	adjustnew=adjust
2005
2006	else ! Pluto
2007	! update new pday and tpal (Myr) to be set in startfi controle
2008	pdaypal=int(day_ini+paleoyears*365.25/6.3872)
2009	tpalnew=tpal+paleoyears*1.e-6 ! Myrs
2010
2011	! update new N2 ice adjustment (not tested yet on Pluto)
2012	adjustnew=adjust
2013
2014	! update milankovitch parameters : obliquity,Lsp,ecc
2015	call calcmilank(tpalnew,oblipal,peri_daypal,eccpal)
2016	!peri_daypal=peri_day
2017	!eccpal=0.009
2018
2019	endif
2020
2021	if (is_master) write(,) "Paleo peri=",peri_daypal," tpal=",tpalnew
2022	if (is_master) write(,) "Paleo eccpal=",eccpal," tpal=",tpalnew
2023
2024
2025	#ifndef MESOSCALE
2026	! create restartfi
2027	if (ngrid.ne.1) then
2028	!TODO: import this routine from pluto.old
2029	! call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq, &
2030	! ptimestep,pdaypal, &
2031	! ztime_fin,tsurf,tsoil,emis,q2,qsurfpal, &
2032	! cell_area,albedodat,tidat,zmea,zstd,zsig, &
2033	! zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal, &
2034	! peri_daypal)
2035	endif
2036	else ! 'paleo'
2037
2038	if (ngrid.ne.1) then
2039	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
2040
2041	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
2042	ptimestep,ztime_fin, &
2043	tsurf,tsoil,emis,q2,qsurf_hist)
2044	endif
2045	#endif
2046	endif ! end of 'paleo'
2047	! if(ok_slab_ocean) then
2048	! call ocean_slab_final!(tslab, seaice)
2049	! end if
2050
2051	endif ! end of 'lastcall'
2052
2053
2054	! -----------------------------------------------------------------
2055	! WSTATS: Saving statistics
2056	! -----------------------------------------------------------------
2057	! ("stats" stores and accumulates key variables in file "stats.nc"
2058	! which can later be used to make the statistic files of the run:
2059	! if flag "callstats" from callphys.def is .true.)
2060
2061
2062	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
2063	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2064	call wstats(ngrid,"fluxsurf_lw", &
2065	"Thermal IR radiative flux to surface","W.m-2",2, &
2066	fluxsurf_lw)
2067	call wstats(ngrid,"fluxtop_lw", &
2068	"Thermal IR radiative flux to space","W.m-2",2, &
2069	fluxtop_lw)
2070
2071	! call wstats(ngrid,"fluxsurf_sw", &
2072	! "Solar radiative flux to surface","W.m-2",2, &
2073	! fluxsurf_sw_tot)
2074	! call wstats(ngrid,"fluxtop_sw", &
2075	! "Solar radiative flux to space","W.m-2",2, &
2076	! fluxtop_sw_tot)
2077
2078
2079	call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
2080	call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
2081	call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
2082	!call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
2083	!call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
2084	call wstats(ngrid,"p","Pressure","Pa",3,pplay)
2085	call wstats(ngrid,"emis","Emissivity","",2,emis)
2086	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
2087	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
2088	call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
2089	call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
2090	call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
2091
2092	if (tracer) then
2093	do iq=1,nq
2094	call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
2095	call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2096	'kg m^-2',2,qsurf(1,iq) )
2097	call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2098	'kg m^-2',2,qcol(1,iq) )
2099
2100	! call wstats(ngrid,trim(noms(iq))//'_reff', &
2101	! trim(noms(iq))//'_reff', &
2102	! 'm',3,reffrad(1,1,iq))
2103
2104	end do
2105
2106	endif ! end of 'tracer'
2107
2108	!AF24: deleted slab ocean and water
2109
2110	if(lastcall.and.callstats) then
2111	write (,) "Writing stats..."
2112	call mkstats(ierr)
2113	endif
2114
2115
2116	#ifndef MESOSCALE
2117
2118	!-----------------------------------------------------------------------------------------------------
2119	! OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
2120	!
2121	! Note 1 : output with period "ecritphy", set in "run.def"
2122	!
2123	! Note 2 : writediagfi can also be called from any other subroutine for any variable,
2124	! but its preferable to keep all the calls in one place ...
2125	!-----------------------------------------------------------------------------------------------------
2126
2127	call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
2128	call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
2129	call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
2130	call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
2131	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2132	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
2133
2134	!! Pluto outputs
2135	! call writediagfi(ngrid,"rice_ch4","ch4 ice mass mean radius","m",3,rice_ch4)
2136	call writediagfi(ngrid,"dist_star","dist_star","AU",0,dist_star)
2137
2138	if (fast) then
2139	call writediagfi(ngrid,"globave","surf press","Pa",0,globave)
2140	!AF: TODO which outputs?
2141	else
2142	if (check_physics_outputs) then
2143	! Check the validity of updated fields at the end of the physics step
2144	call check_physics_fields("HERE physiq:", zt, zu, zv, pplev, zq)
2145	endif
2146
2147	call writediagfi(ngrid,"emis","Emissivity","",2,emis)
2148	call writediagfi(ngrid,"temp","temperature","K",3,zt)
2149	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
2150	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
2151	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
2152	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
2153	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
2154	endif
2155
2156	! Subsurface temperatures
2157	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2158	! call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
2159
2160	! ! Oceanic layers !AF24: removed
2161
2162	! ! Thermal plume model !AF24: removed
2163
2164	! GW non-oro outputs !AF24: removed
2165
2166	! Total energy balance diagnostics
2167	if(callrad)then
2168
2169	call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
2170	call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
2171	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
2172	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
2173	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
2174	call writediagfi(ngrid,"shad","rings"," ", 2, fract)
2175
2176	! call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
2177	! call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
2178	! call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
2179	! call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
2180	! call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
2181	! call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
2182
2183	! if(ok_slab_ocean) then
2184	! call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrdocean)
2185	! else
2186	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
2187	! endif
2188
2189	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
2190
2191	endif ! end of 'callrad'
2192
2193	if(enertest) then
2194
2195	if (calldifv) then
2196
2197	call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
2198	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
2199
2200	! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
2201	! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
2202	! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
2203
2204	endif
2205
2206	if (corrk) then
2207	call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
2208	call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
2209	endif
2210
2211	! if(watercond) then !AF24: removed
2212
2213	if (generic_condensation) then
2214
2215	call writediagfi(ngrid,"genericconddE","heat from generic condensation","W m-2",2,genericconddE)
2216	call writediagfi(ngrid,"dt_generic_condensation","heating from generic condensation","K s-1",3,dt_generic_condensation)
2217
2218	endif
2219
2220	endif ! end of 'enertest'
2221
2222	! Diagnostics of optical thickness
2223	! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
2224	if (diagdtau) then
2225	do nw=1,L_NSPECTV
2226	write(str2,'(i2.2)') nw
2227	call writediagfi(ngrid,'dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',1,int_dtauv(:,nlayer:1:-1,nw))
2228	enddo
2229	do nw=1,L_NSPECTI
2230	write(str2,'(i2.2)') nw
2231	call writediagfi(ngrid,'dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',1,int_dtaui(:,nlayer:1:-1,nw))
2232	enddo
2233	endif
2234
2235
2236	! Temporary inclusions for heating diagnostics.
2237	call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
2238	call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
2239	call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
2240	call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
2241
2242	!Pluto specific
2243	call writediagfi(ngrid,"zdtc","tendancy T cond N2","K",3,zdtc)
2244	! call writediagfi(ngrid,"zdtch4cloud","tendancy T ch4cloud","K",3,zdtch4cloud)
2245	! call writediagfi(ngrid,"zdtcocloud","tendancy T cocloud","K",3,zdtcocloud)
2246	! call writediagfi(ngrid,"zq1temp_ch4"," "," ",2,zq1temp_ch4)
2247	! call writediagfi(ngrid,"qsat_ch4"," "," ",2,qsat_ch4)
2248	! call writediagfi(ngrid,"qsat_ch4_l1"," "," ",2,qsat_ch4_l1)
2249	! call writediagfi(ngrid,"senshf1","senshf1"," ",2,sensiblehf1)
2250	! call writediagfi(ngrid,"senshf2","senshf2"," ",2,sensiblehf2)
2251
2252
2253	! For Debugging.
2254	!call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
2255	!call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
2256
2257
2258	! Output aerosols.!AF: TODO: write haze aerosols
2259	! if (igcm_n2_ice.ne.0.and.iaero_haze.ne.0) &
2260	! call writediagfi(ngrid,'N2ice_reff','N2ice_reff','m',3,reffrad(1,1,iaero_haze))
2261	! if (igcm_n2_ice.ne.0.and.iaero_haze.ne.0) &
2262	! call writediagfi(ngrid,'N2ice_reffcol','N2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_haze))
2263	! if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) & !AF24: removed
2264
2265	! Output tracers.
2266	if (tracer) then
2267
2268	do iq=1,nq
2269	call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
2270	call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2271	'kg m^-2',2,qsurf_hist(1,iq) )
2272	call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2273	'kg m^-2',2,qcol(1,iq) )
2274	! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2275	! 'kg m^-2',2,qsurf(1,iq) )
2276
2277	! if(watercond.or.CLFvarying)then !AF24: removed
2278
2279	if(generic_condensation)then
2280	call writediagfi(ngrid,"rneb_generic","GCS cloud fraction (generic condensation)"," ",3,rneb_generic)
2281	call writediagfi(ngrid,"CLF","GCS cloud fraction"," ",3,cloudfrac)
2282	call writediagfi(ngrid,"RH_generic","GCS relative humidity"," ",3,RH_generic)
2283	endif
2284
2285	! if(generic_rain)then !AF24: removed
2286	! if((hydrology).and.(.not.ok_slab_ocean))then !AF24: removed
2287
2288	call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
2289
2290	enddo ! end of 'nq' loop
2291
2292	!Pluto specific
2293	call writediagfi(ngrid,'n2_iceflux','n2_iceflux',"kg m^-2 s^-1",2,flusurf(1,igcm_n2) )
2294	! call writediagfi(ngrid,'haze_reff','haze_reff','m',3,reffrad(1,1,1))
2295	if (methane) then
2296	call writediagfi(ngrid,'ch4_iceflux','ch4_iceflux',&
2297	"kg m^-2 s^-1",2,flusurf(1,igcm_ch4_ice) )
2298	call writediagfi(ngrid,"vmr_ch4","vmr_ch4","%",2,vmr_ch4)
2299	if (.not.fast) then
2300	call writediagfi(ngrid,"zrho_ch4","zrho_ch4","kg.m-3",3,zrho_ch4(:,:))
2301	endif
2302
2303	! Tendancies
2304	call writediagfi(ngrid,"zdqch4cloud","ch4 cloud","T s-1",&
2305	3,zdqch4cloud(1,1,igcm_ch4_gas))
2306	call writediagfi(ngrid,"zdqcn2_ch4","zdq condn2 ch4","",&
2307	3,zdqc(:,:,igcm_ch4_gas))
2308	call writediagfi(ngrid,"zdqdif_ch4","zdqdif ch4","",&
2309	3,zdqdif(:,:,igcm_ch4_gas))
2310	call writediagfi(ngrid,"zdqsdif_ch4_ice","zdqsdif ch4","",&
2311	2,zdqsdif(:,igcm_ch4_ice))
2312	call writediagfi(ngrid,"zdqadj_ch4","zdqadj ch4","",&
2313	3,zdqadj(:,:,igcm_ch4_gas))
2314	if (sedimentation) then
2315	call writediagfi(ngrid,"zdqsed_ch4","zdqsed ch4","",&
2316	3,zdqsed(:,:,igcm_ch4_gas))
2317	call writediagfi(ngrid,"zdqssed_ch4","zdqssed ch4","",&
2318	2,zdqssed(:,igcm_ch4_gas))
2319	endif
2320	if (metcloud) then
2321	call writediagfi(ngrid,"zdtch4cloud","ch4 cloud","T s-1",&
2322	3,zdtch4cloud)
2323	endif
2324
2325	endif
2326
2327	if (carbox) then
2328	call writediagfi(ngrid,'co_iceflux','co_iceflux',&
2329	"kg m^-2 s^-1",2,flusurf(1,igcm_co_ice) )
2330	call writediagfi(ngrid,"vmr_co","vmr_co","%",2,vmr_co)
2331	if (.not.fast) THEN
2332	call writediagfi(ngrid,"zrho_co","zrho_co","kg.m-3",3,zrho_co(:,:))
2333	endif
2334	endif
2335
2336	if (haze) then
2337	! call writediagfi(ngrid,"zrho_haze","zrho_haze","kg.m-3",3,zrho_haze(:,:))
2338	call writediagfi(ngrid,"zdqrho_photprec","zdqrho_photprec",&
2339	"kg.m-3.s-1",3,zdqrho_photprec(:,:))
2340	call writediagfi(ngrid,"zdqphot_prec","zdqphot_prec","",&
2341	3,zdqphot_prec(:,:))
2342	call writediagfi(ngrid,"zdqhaze_ch4","zdqhaze_ch4","",&
2343	3,zdqhaze(:,:,igcm_ch4_gas))
2344	call writediagfi(ngrid,"zdqhaze_prec","zdqhaze_prec","",&
2345	3,zdqhaze(:,:,igcm_prec_haze))
2346	if (igcm_haze.ne.0) then
2347	call writediagfi(ngrid,"zdqhaze_haze","zdqhaze_haze","",&
2348	3,zdqhaze(:,:,igcm_haze))
2349	if (sedimentation) then
2350	call writediagfi(ngrid,"zdqssed_haze","zdqssed haze",&
2351	"kg/m2/s",2,zdqssed(:,igcm_haze))
2352	endif
2353	endif
2354	call writediagfi(ngrid,"zdqphot_ch4","zdqphot_ch4","",&
2355	3,zdqphot_ch4(:,:))
2356	call writediagfi(ngrid,"zdqconv_prec","zdqconv_prec","",&
2357	3,zdqconv_prec(:,:))
2358	! call writediagfi(ngrid,"zdqhaze_col","zdqhaze col","kg/m2/s",
2359	! & 2,zdqhaze_col(:))
2360	endif
2361
2362	if (aerohaze) then
2363	call writediagfi(ngrid,"tau_col",&
2364	"Total aerosol optical depth","opacity",2,tau_col)
2365	endif
2366
2367	endif ! end of 'tracer'
2368
2369
2370	! Output spectrum.
2371	if(specOLR.and.corrk)then
2372	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
2373	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
2374	call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W/m^2/cm^-1",3,GSR_nu)
2375	endif
2376
2377	#else
2378	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
2379	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
2380	comm_ALBEQ(1:ngrid)=albedo_equivalent(1:ngrid)
2381	if (.not.calldifv) comm_LATENT_HF(:)=0.0
2382	! if ((tracer).and.(water)) then !AF24: removed
2383
2384	if ((tracer).and.(generic_condensation)) then
2385	! .and.(.not. water)
2386
2387	! If you have set generic_condensation (and not water) and you have set several GCS
2388	! then the outputs given to WRF will be only the ones for the last generic tracer
2389	! (Because it is rewritten every tracer in the loop)
2390	! WRF can take only one moist tracer
2391
2392	do iq=1,nq
2393	call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
2394
2395	if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
2396
2397	reffrad_generic_zeros_for_wrf(:,:) = 1.
2398
2399	comm_CLOUDFRAC(1:ngrid,1:nlayer) = cloudfrac(1:ngrid,1:nlayer)
2400	comm_TOTCLOUDFRAC(1:ngrid) = totcloudfrac(1:ngrid) !??????
2401	comm_SURFRAIN(1:ngrid) = zdqsrain_generic(1:ngrid,iq)
2402	comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_gas)
2403	comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_generic_ice)
2404	! comm_H2OICE_REFF(1:ngrid,1:nlayer) = reffrad_generic_zeros_for_wrf(1:ngrid,1:nlayer) ! for now zeros !
2405	!comm_H2OICE_REFF(1:ngrid,1:nlayer) = 0*zdtrain_generic(1:ngrid,1:nlayer) ! for now zeros !
2406	comm_REEVAP(1:ngrid) = reevap_precip_generic(1:ngrid,iq)
2407	comm_DTRAIN(1:ngrid,1:nlayer) = zdtrain_generic(1:ngrid,1:nlayer)
2408	comm_DTLSC(1:ngrid,1:nlayer) = dt_generic_condensation(1:ngrid,1:nlayer)
2409	comm_RH(1:ngrid,1:nlayer) = RH_generic(1:ngrid,1:nlayer,iq)
2410
2411	endif
2412	end do ! do iq=1,nq loop on tracers
2413
2414	else
2415	comm_CLOUDFRAC(1:ngrid,1:nlayer)=0.
2416	comm_TOTCLOUDFRAC(1:ngrid)=0.
2417	comm_SURFRAIN(1:ngrid)=0.
2418	comm_DQVAP(1:ngrid,1:nlayer)=0.
2419	comm_DQICE(1:ngrid,1:nlayer)=0.
2420	! comm_H2OICE_REFF(1:ngrid,1:nlayer)=0.
2421	comm_REEVAP(1:ngrid)=0.
2422	comm_DTRAIN(1:ngrid,1:nlayer)=0.
2423	comm_DTLSC(1:ngrid,1:nlayer)=0.
2424	comm_RH(1:ngrid,1:nlayer)=0.
2425
2426	endif ! if water, if generic_condensation, else
2427
2428	comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid)
2429	comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid)
2430	comm_FLUXTOP_LW(1:ngrid)=fluxtop_lw(1:ngrid)
2431	comm_FLUXSURF_SW(1:ngrid)=fluxsurf_sw(1:ngrid)
2432	comm_FLUXSURF_LW(1:ngrid)=fluxsurf_lw(1:ngrid)
2433	comm_FLXGRD(1:ngrid)=fluxgrd(1:ngrid)
2434	sensibFlux(1:ngrid) = zflubid(1:ngrid) - capcal(1:ngrid)*zdtsdif(1:ngrid) !!! ????
2435	comm_HR_DYN(1:ngrid,1:nlayer) = zdtdyn(1:ngrid,1:nlayer)
2436
2437	! if (turb_resolved) then
2438	! open(17,file='lsf.txt',form='formatted',status='old')
2439	! rewind(17)
2440	! DO l=1,nlayer
2441	! read(17,*) lsf_dt(l),lsf_dq(l)
2442	! ENDDO
2443	! close(17)
2444	! do ig=1,ngrid
2445	! if ((tracer).and.(water)) then
2446	! pdq(ig,:,igcm_h2o_gas) = pdq(ig,:,igcm_h2o_gas) + lsf_dq(:)
2447	! endif
2448	! pdt(ig,:) = pdt(ig,:) + lsf_dt(:)
2449	! comm_HR_DYN(ig,:) = lsf_dt(:)
2450	! enddo
2451	! endif
2452	#endif
2453
2454	! XIOS outputs
2455	#ifdef CPP_XIOS
2456	! Send fields to XIOS: (NB these fields must also be defined as
2457	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
2458	CALL send_xios_field("ls",zls)
2459
2460	CALL send_xios_field("ps",ps)
2461	CALL send_xios_field("area",cell_area)
2462	CALL send_xios_field("p",pplay)
2463	CALL send_xios_field("temperature",zt)
2464	CALL send_xios_field("u",zu)
2465	CALL send_xios_field("v",zv)
2466	CALL send_xios_field("omega",omega)
2467
2468	! IF (calltherm) THEN !AF24: removed
2469	! IF (water) THEN !AF24: removed
2470
2471	CALL send_xios_field("ISR",fluxtop_dn)
2472	CALL send_xios_field("OLR",fluxtop_lw)
2473	CALL send_xios_field("ASR",fluxabs_sw)
2474
2475	if (specOLR .and. corrk) then
2476	call send_xios_field("OLR3D",OLR_nu)
2477	call send_xios_field("IR_Bandwidth",DWNI)
2478	call send_xios_field("VI_Bandwidth",DWNV)
2479	call send_xios_field("OSR3D",OSR_nu)
2480	call send_xios_field("GSR3D",GSR_nu)
2481	endif
2482
2483	if (lastcall.and.is_omp_master) then
2484	write(,) "physiq: call xios_context_finalize"
2485	call xios_context_finalize
2486	endif
2487	#endif
2488
2489	if (check_physics_outputs) then
2490	! Check the validity of updated fields at the end of the physics step
2491	call check_physics_fields("end of physiq:", zt, zu, zv, pplev, zq)
2492	endif
2493
2494	icount=icount+1
2495
2496	end subroutine physiq
2497
2498	end module physiq_mod

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

Download in other formats:

Original Format