Context Navigation

physiq_mod.F90 @ 3506

Last change on this file since 3506 was 3506, checked in by afalco, 13 days ago
Pluto: import write_output function from Mars. xios specific outputs. AF
Property svn:executable set to ``*
File size: 108.4 KB

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

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

Download in other formats:

Original Format