Context Navigation

physiq_mod.F90 @ 3508

Last change on this file since 3508 was 3508, checked in by afalco, 13 days ago
Pluto: physiq_mod to call write_output, which writes both via XIOS and diagfi routines. Included xml example files for XIOS. AF
Property svn:executable set to ``*
File size: 106.8 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	use write_output_mod, only: write_output
103
104	implicit none
105
106
107	!==================================================================
108	!
109	! Purpose
110	! -------
111	! Central subroutine for all the physics parameterisations in the
112	! universal model. Originally adapted from the Mars LMDZ model.
113	!
114	! The model can be run without or with tracer transport
115	! depending on the value of "tracer" in file "callphys.def".
116	!
117	!
118	! It includes:
119	!
120	! I. Initialization :
121	! I.1 Firstcall initializations.
122	! I.2 Initialization for every call to physiq.
123	!
124	! II. Compute radiative transfer tendencies (longwave and shortwave) :
125	! II.a Option 1 : Call correlated-k radiative transfer scheme.
126	! II.b Option 2 : Call Newtonian cooling scheme.
127	! II.! Option 3 : Atmosphere has no radiative effect.
128	!
129	! III. Vertical diffusion (turbulent mixing) :
130	!
131	! IV. Convection :
132	! IV.a Thermal plume model
133	! IV.b Dry convective adjusment
134	!
135	! V. Condensation and sublimation of gases (currently just N2).
136	!
137	! VI. Tracers
138	! VI.1. Water and water ice.
139	! VI.2 Photochemistry
140	! VI.3. Aerosols and particles.
141	! VI.4. Updates (pressure variations, surface budget).
142	! VI.5. Slab Ocean.
143	! VI.6. Surface Tracer Update.
144	!
145	! VII. Surface and sub-surface soil temperature.
146	!
147	! VIII. Perform diagnostics and write output files.
148	!
149	!
150	! arguments
151	! ---------
152	!
153	! INPUT
154	! -----
155	!
156	! ngrid Size of the horizontal grid.
157	! nlayer Number of vertical layers.
158	! nq Number of advected fields.
159	!
160	! firstcall True at the first call.
161	! lastcall True at the last call.
162	!
163	! pday Number of days counted from the North. Spring equinoxe.
164	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
165	! ptimestep timestep (s).
166	!
167	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa).
168	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
169	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2.s-2).
170	!
171	! pu(ngrid,nlayer) u, zonal component of the wind (ms-1).
172	! pv(ngrid,nlayer) v, meridional component of the wind (ms-1).
173	!
174	! pt(ngrid,nlayer) Temperature (K).
175	!
176	! pq(ngrid,nlayer,nq) Advected fields.
177	!
178	! pudyn(ngrid,nlayer) \
179	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
180	! ptdyn(ngrid,nlayer) / corresponding variables.
181	! pqdyn(ngrid,nlayer,nq) /
182	! flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
183	!
184	! OUTPUT
185	! ------
186	!
187	! pdu(ngrid,nlayer) \
188	! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding
189	! pdt(ngrid,nlayer) / variables due to physical processes.
190	! pdq(ngrid,nlayer) /
191	! pdpsrf(ngrid) /
192	!
193	!
194	! Authors
195	! -------
196	! Frederic Hourdin 15/10/93
197	! Francois Forget 1994
198	! Christophe Hourdin 02/1997
199	! Subroutine completely rewritten by F. Forget (01/2000)
200	! Water ice clouds: Franck Montmessin (update 06/2003)
201	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
202	! New correlated-k radiative scheme: R. Wordsworth (2009)
203	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
204	! Improved water cycle: R. Wordsworth / B. Charnay (2010)
205	! To F90: R. Wordsworth (2010)
206	! New turbulent diffusion scheme: J. Leconte (2012)
207	! Loops converted to F90 matrix format: J. Leconte (2012)
208	! No more ngrid/nq, F90 commons and adaptation to parallel: A. Spiga (2012)
209	! Purge of the code : M. Turbet (2015)
210	! Photochemical core developped by F. Lefevre: B. Charnay (2017)
211	! Purge for Pluto model : A. Falco (2024)
212	!==================================================================
213
214
215	! 0. Declarations :
216	! ------------------
217
218	include "netcdf.inc"
219
220	! Arguments :
221	! -----------
222
223	! INPUTS:
224	! -------
225
226	integer,intent(in) :: ngrid ! Number of atmospheric columns.
227	integer,intent(in) :: nlayer ! Number of atmospheric layers.
228	integer,intent(in) :: nq ! Number of tracers.
229
230	logical,intent(in) :: firstcall ! Signals first call to physics.
231	logical,intent(in) :: lastcall ! Signals last call to physics.
232
233	real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0.
234	real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
235	real,intent(in) :: ptimestep ! Physics timestep (s).
236	real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
237	real,intent(in) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
238	real,intent(in) :: pphi(ngrid,nlayer) ! Geopotential at mid-layer (m2s-2).
239	real,intent(in) :: pu(ngrid,nlayer) ! Zonal wind component (m/s).
240	real,intent(in) :: pv(ngrid,nlayer) ! Meridional wind component (m/s).
241	real,intent(in) :: pt(ngrid,nlayer) ! Temperature (K).
242	real,intent(in) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
243	real,intent(in) :: flxw(ngrid,nlayer) ! Vertical mass flux (ks/s) at lower boundary of layer
244
245	! OUTPUTS:
246	! --------
247
248	! Physical tendencies :
249
250	real,intent(out) :: pdu(ngrid,nlayer) ! Zonal wind tendencies (m/s/s).
251	real,intent(out) :: pdv(ngrid,nlayer) ! Meridional wind tendencies (m/s/s).
252	real,intent(out) :: pdt(ngrid,nlayer) ! Temperature tendencies (K/s).
253	real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
254	real,intent(out) :: pdpsrf(ngrid) ! Surface pressure tendency (Pa/s).
255
256	! Local saved variables:
257	! ----------------------
258	integer,save :: day_ini ! Initial date of the run (sol since Ls=0).
259	integer,save :: icount ! Counter of calls to physiq during the run.
260	!$OMP THREADPRIVATE(day_ini,icount)
261
262	!Pluto specific
263	REAL,save :: acond,bcond
264	REAL,save :: tcond1p4Pa
265	DATA tcond1p4Pa/38/
266
267
268
269	! Local variables :
270	! -----------------
271	! Tendencies for the paleoclimate mode
272	REAL qsurfyear(ngrid,nq) ! kg.m-2 averaged mass of ice lost/gained in the last Pluto year of the run
273	REAL phisfinew(ngrid) ! geopotential of the bedrock (= phisfi-qsurf/1000*g)
274	REAL qsurfpal(ngrid,nq) ! qsurf after a paleoclimate step : for physdem1 and restartfi
275	REAL phisfipal(ngrid) ! geopotential after a paleoclimate step : for physdem1 and restartfi
276	REAL oblipal ! change of obliquity
277	REAL peri_daypal ! new periday
278	REAL eccpal ! change of eccentricity
279	REAL tpalnew ! change of time
280	REAL adjustnew ! change in N2 ice albedo
281	REAL pdaypal ! new pday = day_ini + step
282	REAL zdt_tot ! time range corresponding to the flux of qsurfyear
283	REAL massacc(nq) ! accumulated mass
284	REAL masslost(nq) ! accumulated mass
285
286	REAL globave ! globalaverage 2D ps
287	REAL globaveice(nq) ! globalaverage 2D ice
288	REAL globavenewice(nq) ! globalaverage 2D ice
289	INTEGER lecttsoil ! lecture of tsoil from proftsoil
290	REAL qsurf1(ngrid,nq) ! saving qsurf to calculate flux over long timescales kg.m-2
291	REAL flusurf(ngrid,nq) ! flux cond/sub kg.m-2.s-1
292	REAL flusurfold(ngrid,nq) ! old flux cond/sub kg.m-2.s-1
293	REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
294
295	REAL,SAVE :: ptime0 ! store the first time
296	REAL dstep
297	REAL,SAVE :: glastep=20 ! step in pluto day to spread glacier
298
299
300	! Aerosol (dust or ice) extinction optical depth at reference wavelength
301	! for the "naerkind" optically active aerosols:
302
303	real,save,allocatable :: aerosol(:,:,:) ! Aerosols
304	!$OMP THREADPRIVATE(aerosol)
305	real zh(ngrid,nlayer) ! Potential temperature (K).
306	real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
307	real omega(ngrid,nlayer) ! omega velocity (Pa/s, >0 when downward)
308
309	integer i,l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_gas, igcm_generic_ice
310	logical call_ice_gas_generic ! to call only one time the ice/vap pair of a tracer
311
312	real zls ! Solar longitude (radians).
313	real zlss ! Sub solar point longitude (radians).
314	real zday ! Date (time since Ls=0, calculated in sols).
315	REAL,save :: saveday ! saved date
316	REAL,save :: savedeclin ! saved declin
317	real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers.
318	real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries.
319
320
321	! VARIABLES for the thermal plume model (AF24: deleted)
322
323	! TENDENCIES due to various processes :
324
325	! For Surface Temperature : (K/s)
326	real zdtsurf(ngrid) ! Cumulated tendencies.
327	real zdtsurfmr(ngrid) ! Mass_redistribution routine.
328	real zdtsurfc(ngrid) ! Condense_n2 routine.
329	real zdtsdif(ngrid) ! Turbdiff/vdifc routines.
330	! real zdtsurf_hyd(ngrid) ! Hydrol routine.
331
332	! For Atmospheric Temperatures : (K/s)
333	real dtlscale(ngrid,nlayer) ! Largescale routine.
334	real dt_generic_condensation(ngrid,nlayer) ! condensation_generic routine.
335	real zdtc(ngrid,nlayer) ! Condense_n2 routine.
336	real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
337	real zdtmr(ngrid,nlayer) ! Mass_redistribution routine.
338	real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
339	real zdtchim(ngrid,nlayer) ! Calchim routine.
340
341	! For Surface Tracers : (kg/m2/s)
342	real dqsurf(ngrid,nq) ! Cumulated tendencies.
343	!real zdqsurfc(ngrid) ! Condense_n2 routine.
344	REAL zdqsc(ngrid,nq) ! Condense_n2 routine.
345	real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines.
346	real zdqssed(ngrid,nq) ! Callsedim routine.
347	real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine.
348
349	! For Tracers : (kg/kg_of_air/s)
350	real zdqc(ngrid,nlayer,nq) ! Condense_n2 routine.
351	real zdqadj(ngrid,nlayer,nq) ! Convadj routine.
352	real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines.
353	real zdqevap(ngrid,nlayer) ! Turbdiff routine.
354	real zdqsed(ngrid,nlayer,nq) ! Callsedim routine.
355	real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine.
356	real dqvaplscale_generic(ngrid,nlayer,nq) ! condensation_generic routine.
357	real dqcldlscale_generic(ngrid,nlayer,nq) ! condensation_generic routine.
358	REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine
359	REAL,allocatable,save :: zdqschim(:,:) ! Calchim_asis routine
360	!$OMP THREADPRIVATE(zdqchim,zdqschim)
361
362
363	!! PLUTO variables
364	REAL zdqch4cloud(ngrid,nlayer,nq)
365	REAL zdqsch4cloud(ngrid,nq)
366	REAL zdtch4cloud(ngrid,nlayer)
367	REAL zdqcocloud(ngrid,nlayer,nq)
368	REAL zdqscocloud(ngrid,nq)
369	REAL zdtcocloud(ngrid,nlayer)
370	REAL rice_ch4(ngrid,nlayer) ! Methane ice geometric mean radius (m)
371	REAL rice_co(ngrid,nlayer) ! CO ice geometric mean radius (m)
372
373	REAL zdqsch4fast(ngrid) ! used only for fast model nogcm
374	REAL zdqch4fast(ngrid) ! used only for fast model nogcm
375	REAL zdqscofast(ngrid) ! used only for fast model nogcm
376	REAL zdqcofast(ngrid) ! used only for fast model nogcm
377	REAL zdqflow(ngrid,nq)
378
379	REAL zdtconduc(ngrid,nlayer) ! (K/s)
380	REAL zdumolvis(ngrid,nlayer)
381	REAL zdvmolvis(ngrid,nlayer)
382	real zdqmoldiff(ngrid,nlayer,nq)
383
384	! Haze relatated tendancies
385	REAL zdqhaze(ngrid,nlayer,nq)
386	REAL zdqprodhaze(ngrid,nq)
387	REAL zdqsprodhaze(ngrid)
388	REAL zdqhaze_col(ngrid)
389	REAL zdqphot_prec(ngrid,nlayer)
390	REAL zdqphot_ch4(ngrid,nlayer)
391	REAL zdqconv_prec(ngrid,nlayer)
392	REAL zdq_source(ngrid,nlayer,nq)
393	! Fast Haze relatated tendancies
394	REAL fluxbot(ngrid)
395	REAL gradflux(ngrid)
396	REAL fluxlym_sol_bot(ngrid) ! Solar flux Lyman alpha ph.m-2.s-1 reaching the surface
397	REAL fluxlym_ipm_bot(ngrid) ! IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 reaching the surface
398	REAL flym_sol(ngrid) ! Incident Solar flux Lyman alpha ph.m-2.s-1
399	REAL flym_ipm(ngrid) ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1
400	REAL zfluxuv ! Lyman alpha flux at 1AU
401
402
403
404	REAL array_zero1(ngrid)
405	REAL array_zero2(ngrid,nlayer)
406
407	! For Winds : (m/s/s)
408	real zdvadj(ngrid,nlayer), zduadj(ngrid,nlayer) ! Convadj routine.
409	real zdumr(ngrid,nlayer), zdvmr(ngrid,nlayer) ! Mass_redistribution routine.
410	real zdvdif(ngrid,nlayer), zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
411	real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
412	real zdhadj(ngrid,nlayer) ! Convadj routine.
413	REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer) ! condense_n2 routine.
414
415	! For Pressure and Mass :
416	real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
417	real zdmassmr_col(ngrid) ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
418	real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s).
419
420
421
422	! Local variables for LOCAL CALCULATIONS:
423	! ---------------------------------------
424	real zflubid(ngrid)
425	real zplanck(ngrid),zpopsk(ngrid,nlayer)
426	REAL zdum1(ngrid,nlayer)
427	REAL zdum2(ngrid,nlayer)
428	real ztim1,ztim2,ztim3, z1,z2
429	real ztime_fin
430	real zdh(ngrid,nlayer)
431	real gmplanet
432	real taux(ngrid),tauy(ngrid)
433
434
435	! local variables for DIAGNOSTICS : (diagfi & stat)
436	! -------------------------------------------------
437	real ps(ngrid) ! Surface Pressure.
438	real zt(ngrid,nlayer) ! Atmospheric Temperature.
439	real zu(ngrid,nlayer),zv(ngrid,nlayer) ! Zonal and Meridional Winds.
440	real zq(ngrid,nlayer,nq) ! Atmospheric Tracers.
441	real zdtadj(ngrid,nlayer) ! Convadj Diagnostic.
442	real zdtdyn(ngrid,nlayer) ! Dynamical Heating (K/s).
443	real zdudyn(ngrid,nlayer) ! Dynamical Zonal Wind tendency (m.s-2).
444
445	real reff(ngrid,nlayer) ! Effective dust radius (used if doubleq=T).
446	real vmr(ngrid,nlayer) ! volume mixing ratio
447	real time_phys
448
449	real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
450
451	real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
452
453	! Pluto adding variables
454	real vmr_ch4(ngrid) ! vmr ch4
455	real vmr_co(ngrid) ! vmr co
456	real rho(ngrid,nlayer) ! density
457	real zrho_ch4(ngrid,nlayer) ! density methane kg.m-3
458	real zrho_co(ngrid,nlayer) ! density CO kg.m-3
459	real zrho_haze(ngrid,nlayer) ! density haze kg.m-3
460	real zdqrho_photprec(ngrid,nlayer) !photolysis rate kg.m-3.s-1
461	real zq1temp_ch4(ngrid) !
462	real qsat_ch4(ngrid) !
463	real qsat_ch4_l1(ngrid) !
464	! CHARACTER(LEN=20) :: txt ! to temporarly store text for eddy tracers
465	real profmmr(ngrid,nlayer) ! fixed profile of haze if haze_proffix
466	real sensiblehf1(ngrid) ! sensible heat flux
467	real sensiblehf2(ngrid) ! sensible heat flux
468
469	! included by RW for H2O Manabe scheme
470	real rneb_man(ngrid,nlayer) ! H2O cloud fraction (moistadj).
471	real rneb_lsc(ngrid,nlayer) ! H2O cloud fraction (large scale).
472
473
474	! to test energy conservation (RW+JL)
475	real mass(ngrid,nlayer),massarea(ngrid,nlayer)
476	real dEtot, dEtots, AtmToSurf_TurbFlux
477	real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
478	!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
479
480	!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
481
482	real dtmoist_max,dtmoist_min
483	real dItot, dItot_tmp, dVtot, dVtot_tmp
484
485
486	real dWtot, dWtot_tmp, dWtots, dWtots_tmp
487
488	real psat_tmp ! AF24: to remove?
489
490	real qsat_generic(ngrid,nlayer,nq) ! generic condensable tracers (GCS) specific concentration at saturation (kg/kg_of_air).
491	real RH_generic(ngrid,nlayer,nq) ! generic condensable tracers (GCS) Relative humidity.
492	real rneb_generic(ngrid,nlayer,nq) ! GCS cloud fraction (generic condensation).
493	real psat_tmp_generic
494	real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
495	!$OMP THREADPRIVATE(metallicity)
496
497	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
498
499	! For Clear Sky Case. (AF24: deleted)
500
501	real nconsMAX, vdifcncons(ngrid), cadjncons(ngrid) ! Vdfic water conservation test. By RW
502
503	real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW
504
505	real,save,allocatable :: reffcol(:,:)
506	!$OMP THREADPRIVATE(reffcol)
507
508	! Non-oro GW tendencies
509	REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
510	REAL d_t_hin(ngrid,nlayer)
511	! Diagnostics 2D of gw_nonoro
512	REAL zustrhi(ngrid), zvstrhi(ngrid)
513
514
515	real :: tsurf2(ngrid)
516	!! real :: flux_o(ngrid),flux_g(ngrid)
517	real :: flux_g(ngrid)
518	real :: flux_sens_lat(ngrid)
519	real :: qsurfint(ngrid,nq)
520	#ifdef MESOSCALE
521	REAL :: lsf_dt(nlayer)
522	REAL :: lsf_dq(nlayer)
523	#endif
524
525	! local variables for skin depth check
526	real :: inertia_min,inertia_max
527	real :: diurnal_skin ! diurnal skin depth (m)
528	real :: annual_skin ! anuual skin depth (m)
529
530	! when no startfi file is asked for init
531	real alpha,lay1 ! coefficients for building layers
532	integer iloop
533
534	! flags to trigger extra sanity checks
535	logical, save :: check_physics_inputs=.false.
536	logical, save :: check_physics_outputs=.false.
537	!$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)
538
539	! Misc
540	character*2 :: str2
541	character(len=10) :: tmp1
542	character(len=10) :: tmp2
543	!==================================================================================================
544
545	! -----------------
546	! I. INITIALISATION
547	! -----------------
548
549	! --------------------------------
550	! I.1 First Call Initialisation.
551	! --------------------------------
552	if (firstcall) then
553	call getin_p("check_physics_inputs", check_physics_inputs)
554	call getin_p("check_physics_outputs", check_physics_outputs)
555
556	! Allocate saved arrays (except for 1D model, where this has already
557	! been done)
558	#ifndef MESOSCALE
559	if (ngrid>1) call phys_state_var_init(nq)
560	#endif
561
562	! Variables set to 0
563	! ~~~~~~~~~~~~~~~~~~
564	dtrad(:,:) = 0.0
565	fluxrad(:) = 0.0
566	tau_col(:) = 0.0
567	zdtsw(:,:) = 0.0
568	zdtlw(:,:) = 0.0
569
570	! Initialize tracer names, indexes and properties.
571	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
572	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
573	if (tracer) then
574	call initracer(ngrid,nq)
575	! if(photochem) then !AF24: removed
576	endif
577	! Initialize aerosol indexes.
578	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
579	! call iniaerosol
580	! allocate related local arrays
581	! (need be allocated instead of automatic because of "naerkind")
582	allocate(aerosol(ngrid,nlayer,naerkind))
583	allocate(reffcol(ngrid,naerkind))
584
585	#ifdef CPP_XIOS
586	! Initialize XIOS context
587	write(,) "physiq: call wxios_context_init"
588	CALL wxios_context_init
589	#endif
590
591	! Read 'startfi.nc' file.
592	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
593	#ifndef MESOSCALE
594	call phyetat0(startphy_file, &
595	ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
596	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,inertiedat)
597
598	#else
599
600	day_ini = pday
601	#endif
602
603	#ifndef MESOSCALE
604	if (.not.startphy_file) then
605	! starting without startfi.nc and with callsoil
606	! is not yet possible as soildepth default is not defined
607	if (callsoil) then
608	! default mlayer distribution, following a power law:
609	! mlayer(k)=lay1alpha*(k-1/2)
610	lay1=2.e-4
611	alpha=2
612	do iloop=0,nsoilmx-1
613	mlayer(iloop)=lay1(alpha*(iloop-0.5))
614	enddo
615	lay1=sqrt(mlayer(0)*mlayer(1))
616	alpha=mlayer(1)/mlayer(0)
617	do iloop=1,nsoilmx
618	layer(iloop)=lay1(alpha*(iloop-1))
619	enddo
620	endif
621	! additionnal "academic" initialization of physics
622	if (is_master) write(,) "Physiq: initializing tsurf(:) to pt(:,1) !!"
623	tsurf(:)=pt(:,1)
624	if (is_master) write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
625	do isoil=1,nsoilmx
626	tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
627	enddo
628	if (is_master) write(,) "Physiq: initializing day_ini to pday !"
629	day_ini=pday
630	endif
631	#endif
632	if (pday.ne.day_ini) then
633	write(,) "ERROR in physiq.F90:"
634	write(,) "bad synchronization between physics and dynamics"
635	write(,) "dynamics day: ",pday
636	write(,) "physics day: ",day_ini
637	stop
638	endif
639
640	write (,) 'In physiq day_ini =', day_ini
641
642	! Initialize albedo calculation.
643	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
644	albedo(:,:)=0.0
645	albedo_bareground(:)=0.0
646	albedo_snow_SPECTV(:)=0.0
647	albedo_n2_ice_SPECTV(:)=0.0
648
649	ptime0=ptime
650	write (,) 'In physiq ptime0 =', ptime
651
652	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV)
653
654	! Initialize orbital calculation.
655	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
656	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
657
658	savedeclin=0.
659	saveday=pday
660	adjust=0. ! albedo adjustment for convergeps
661
662	! Initialize soil.
663	! ~~~~~~~~~~~~~~~~
664	if (callsoil) then
665	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
666	ptimestep,tsurf,tsoil,capcal,fluxgrd)
667	else ! else of 'callsoil'.
668	print*,'WARNING! Thermal conduction in the soil turned off'
669	capcal(:)=1.e6
670	fluxgrd(:)=intheat
671	print*,'Flux from ground = ',intheat,' W m^-2'
672	endif ! end of 'callsoil'.
673
674	icount=1
675
676	! Initialize surface history variable.
677	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
678	qsurf_hist(:,:)=qsurf(:,:)
679
680	!! call WriteField_phy("post_qsurf_hist_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
681
682	! Initialize variable for dynamical heating and zonal wind tendency diagnostic
683	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
684	ztprevious(:,:)=pt(:,:)
685	zuprevious(:,:)=pu(:,:)
686
687	! Set temperature just above condensation temperature (for Early Mars)
688	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
689	if(nearn2cond) then
690	write(,)' WARNING! Starting at Tcond+1K'
691	do l=1, nlayer
692	do ig=1,ngrid
693	pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4 &
694	-pt(ig,l)) / ptimestep
695	enddo
696	enddo
697	endif
698
699	if(meanOLR)then
700	call system('rm -f rad_bal.out') ! to record global radiative balance.
701	call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
702	call system('rm -f h2o_bal.out') ! to record global hydrological balance.
703	endif
704
705
706	!! Initialize variables for water cycle ! AF24: removed
707
708	! Set metallicity for GCS
709	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
710	metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
711	call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic
712
713	! Set some parameters for the thermal plume model !AF24: removed
714	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
715
716	#ifndef MESOSCALE
717	if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
718	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
719	ptimestep,pday+nday,time_phys,cell_area, &
720	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
721	endif
722
723	!! call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
724	#endif
725	if (corrk) then
726	! We initialise the spectral grid here instead of
727	! at firstcall of callcorrk so we can output XspecIR, XspecVI
728	! when using Dynamico
729	if (is_master) print*, "physiq_mod: Correlated-k data base folder:",trim(datadir)
730	call getin_p("corrkdir",corrkdir)
731	if (is_master) print*,"corrkdir = ", corrkdir
732	write (tmp1, '(i4)') L_NSPECTI
733	write (tmp2, '(i4)') L_NSPECTV
734	banddir=trim(trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2)))
735	banddir=trim(trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir)))
736	call setspi !Basic infrared properties.
737	call setspv ! Basic visible properties.
738	call sugas_corrk ! Set up gaseous absorption properties.
739	if (aerohaze) then
740	call suaer_corrk ! Set up aerosol optical properties.
741	endif
742	endif
743
744	!! call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
745	! XIOS outputs
746	#ifdef CPP_XIOS
747
748	if (is_master) write(,) "physiq: call initialize_xios_output"
749	call initialize_xios_output(pday,ptime,ptimestep,daysec,year_day, &
750	presnivs,pseudoalt,mlayer,WNOI,WNOV)
751	#endif
752
753	!! call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
754
755	write(,) "physiq: end of firstcall"
756	endif ! end of 'firstcall'
757
758	!! call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
759	!! call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
760
761	if (check_physics_inputs) then
762	!check the validity of input fields coming from the dynamics
763	call check_physics_fields("begin physiq:", pt, pu, pv, pplev, pq)
764	endif
765
766	! call writediagfi(ngrid,"pre_physical_rnat"," "," ",2,rnat)
767	! call writediagfi(ngrid,"pre_physical_capcal"," "," ",2,capcal)
768
769	! ------------------------------------------------------
770	! I.2 Initializations done at every physical timestep:
771	! ------------------------------------------------------
772
773	#ifdef CPP_XIOS
774	! update XIOS time/calendar
775	call update_xios_timestep
776	#endif
777
778	! Initialize various variables
779	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
780
781	if ( .not.nearn2cond ) then
782	pdt(1:ngrid,1:nlayer) = 0.0
783	endif
784	zdtsurf(1:ngrid) = 0.0
785	pdq(1:ngrid,1:nlayer,1:nq) = 0.0
786	dqsurf(1:ngrid,1:nq)= 0.0
787	pdu(1:ngrid,1:nlayer) = 0.0
788	pdv(1:ngrid,1:nlayer) = 0.0
789	pdpsrf(1:ngrid) = 0.0
790	zflubid(1:ngrid) = 0.0
791	flux_sens_lat(1:ngrid) = 0.0
792	taux(1:ngrid) = 0.0
793	tauy(1:ngrid) = 0.0
794
795	if (conservn2) then
796	write(,) 'conservn2 iniloop'
797	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
798	endif
799
800	zday=pday+ptime ! Compute time, in sols (and fraction thereof).
801
802	! Compute Stellar Longitude (Ls), and orbital parameters.
803	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
804	if (season) then
805	call stellarlong(zday,zls)
806	else
807	call stellarlong(noseason_day,zls)
808	end if
809
810
811	! Get Lyman alpha flux at specific Ls
812	if (haze) then
813	call lymalpha(zls,zfluxuv)
814	print*, 'Haze lyman-alpha zls,zfluxuv=',zls,zfluxuv
815	end if
816
817
818	IF (triton) then
819	CALL orbitetriton(zls,zday,dist_star,declin)
820	ELSE
821	call orbite(zls,dist_star,declin,right_ascen)
822	ENDIF
823
824
825	if (diurnal) then
826	zlss=-2.pi(zday-.5)
827	else if(diurnal .eqv. .false.) then
828	zlss=9999.
829	endif
830
831
832	glat(:) = g !AF24: removed oblateness
833
834	! Compute geopotential between layers.
835	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
836	zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
837	do l=1,nlayer
838	zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
839	enddo
840
841	zzlev(1:ngrid,1)=0.
842
843	do l=2,nlayer
844	do ig=1,ngrid
845	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
846	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
847	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
848	enddo
849	enddo
850
851	!Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
852
853	zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
854
855	! Compute potential temperature
856	! Note : Potential temperature calculation may not be the same in physiq and dynamic...
857	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
858	do l=1,nlayer
859	do ig=1,ngrid
860	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
861	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
862	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
863	massarea(ig,l)=mass(ig,l)*cell_area(ig)
864	enddo
865	enddo
866
867	! Compute vertical velocity (m/s) from vertical mass flux
868	! w = F / (rhoarea) and rho = P/(rT)
869	! But first linearly interpolate mass flux to mid-layers
870	if (.not.fast) then
871	do l=1,nlayer-1
872	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
873	enddo
874	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
875	do l=1,nlayer
876	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) / &
877	(pplay(1:ngrid,l)*cell_area(1:ngrid))
878	enddo
879	! omega in Pa/s
880	do l=1,nlayer-1
881	omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
882	enddo
883	omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
884	do l=1,nlayer
885	omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid)
886	enddo
887	endif
888
889	if (conservn2) then
890	write(,) 'conservn2 thermo'
891	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
892	endif
893
894
895
896	! --------------------------------------------------------
897	! 1.3 thermosphere
898	! --------------------------------------------------------
899
900	! ajout de la conduction depuis la thermosphere
901	IF (callconduct) THEN
902
903	call conduction (ngrid,nlayer,ptimestep, &
904	pplay,pt,zzlay,zzlev,zdtconduc,tsurf)
905	DO l=1,nlayer
906	DO ig=1,ngrid
907	pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
908	ENDDO
909	ENDDO
910
911	ENDIF
912
913	! ajout de la viscosite moleculaire
914	IF (callmolvis) THEN
915	call molvis(ngrid,nlayer,ptimestep, &
916	pplay,pt,zzlay,zzlev, &
917	zdtconduc,pu,tsurf,zdumolvis)
918	call molvis(ngrid,nlayer,ptimestep, &
919	pplay,pt,zzlay,zzlev, &
920	zdtconduc,pv,tsurf,zdvmolvis)
921
922	DO l=1,nlayer
923	DO ig=1,ngrid
924	! pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
925	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
926	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
927	ENDDO
928	ENDDO
929	ENDIF
930
931	IF (callmoldiff) THEN
932	call moldiff_red(ngrid,nlayer,nq, &
933	pplay,pplev,pt,pdt,pq,pdq,ptimestep, &
934	zzlay,zdtconduc,zdqmoldiff)
935
936	DO l=1,nlayer
937	DO ig=1,ngrid
938	DO iq=1, nq
939	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
940	ENDDO
941	ENDDO
942	ENDDO
943	ENDIF
944
945	if (conservn2) then
946	write(,) 'conservn2 thermo'
947	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
948	endif
949
950
951	!---------------------------------
952	! II. Compute radiative tendencies
953	!---------------------------------
954	! Saving qsurf to compute paleo flux condensation/sublimation
955	DO iq=1, nq
956	DO ig=1,ngrid
957	IF (qsurf(ig,iq).lt.0.) then
958	qsurf(ig,iq)=0.
959	ENDIF
960	qsurf1(ig,iq)=qsurf(ig,iq)
961	ENDDO
962	ENDDO
963
964
965	! Compute local stellar zenith angles
966	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
967	fract = 0
968	if (diurnal) then
969	ztim1=SIN(declin)
970	ztim2=COS(declin)COS(2.pi*(zday-.5))
971	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
972
973	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
974	ztim1,ztim2,ztim3,mu0,fract)
975	else if(diurnal .eqv. .false.) then
976
977	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad)
978	! WARNING: this function appears not to work in 1D
979
980	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
981	mu0 = cos(pi*szangle/180.0)
982	fract= 1/(4*mu0) ! AF24: from pluto.old
983	endif
984
985	endif
986
987
988	! Pluto albedo /IT changes depending on surface ices (only in 3D)
989	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
990	if (ngrid.ne.1) then
991
992	!! Specific to change albedo of N2 so that Psurf
993	!! converges toward 1.4 Pa in "1989" seasons for Triton
994	!! converges toward 1.1 Pa in "2015" seasons for Pluto
995	if (convergeps) then
996	if (triton) then
997	! 1989 declination
998	if (declin180./pi.gt.-46..and.declin180./pi.lt.-45. &
999	.and.zday.gt.saveday+1000. &
1000	.and.declin.lt.savedeclin) then
1001	call globalaverage2d(ngrid,pplev(:,1),globave)
1002	if (globave.gt.1.5) then
1003	adjust=adjust+0.005
1004	else if (globave.lt.1.3) then
1005	adjust=adjust-0.005
1006	endif
1007	saveday=zday
1008	endif
1009	else
1010	! Pluto : 2015 declination current epoch
1011	if (declin180./pi.gt.50.and.declin180./pi.lt.51. &
1012	.and.zday.gt.saveday+10000. &
1013	.and.declin.gt.savedeclin) then
1014	call globalaverage2d(ngrid,pplev(:,1),globave)
1015	if (globave.gt.1.2) then
1016	adjust=adjust+0.005
1017	else if (globave.lt.1.) then
1018	adjust=adjust-0.005
1019	endif
1020	saveday=zday
1021	endif
1022	endif
1023	end if
1024	end if ! if ngrid ne 1
1025
1026	call surfprop(ngrid,nq,fract,qsurf,tsurf, &
1027	capcal,adjust,dist_star,flusurfold,ptimestep,zls,&
1028	albedo,emis,inertiedat)
1029	! do i=2,ngrid
1030	! albedo(i,:) = albedo(1,:)
1031	! enddo
1032
1033	if (firstcall.and.callsoil) then
1034	! AF24 Originally in soil.F, but inertiedat is modified by surfprop
1035	! Additional checks: is the vertical discretization sufficient
1036	! to resolve diurnal and annual waves?
1037	do ig=1,ngrid
1038	! extreme inertia for this column
1039	inertia_min=minval(inertiedat(ig,:))
1040	inertia_max=maxval(inertiedat(ig,:))
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)) ! bug (2): Subscript #3 of the array REFFRAD has value 2 which is greater than the upper bound of 1
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 write_output("Ls","solar longitude","deg",zls*180./pi)
2360	! call write_output("Lss","sub solar longitude","deg",zlss*180./pi)
2361	call write_output("RA","right ascension","deg",right_ascen*180./pi)
2362	call write_output("Declin","solar declination","deg",declin*180./pi)
2363	call write_output("dist_star","dist_star","AU",dist_star)
2364
2365	call write_output("tsurf","Surface temperature","K",tsurf)
2366	call write_output("ps","Surface pressure","Pa",ps)
2367	call write_output("emis","Emissivity","",emis)
2368
2369	if (fast) then
2370	call write_output("globave","surf press","Pa",globave)
2371	call write_output("fluxrad","fluxrad","W m-2",fluxrad)
2372	call write_output("fluxgrd","fluxgrd","W m-2",fluxgrd)
2373	call write_output("capcal","capcal","W.s m-2 K-1",capcal)
2374	! call write_output("dplanck","dplanck","W.s m-2 K-1",dplanck)
2375	call write_output("tsoil","tsoil","K",tsoil)
2376	else
2377	if (check_physics_outputs) then
2378	! Check the validity of updated fields at the end of the physics step
2379	call check_physics_fields("physiq:", zt, zu, zv, pplev, zq)
2380	endif
2381
2382	call write_output("zzlay","Midlayer altitude", "m",zzlay(:,:))
2383	call write_output("zzlev","Interlayer altitude", "m",zzlev(:,1:nlayer))
2384	!call write_output('pphi','Geopotential',' ',pphi)
2385
2386	call write_output("temperature","temperature","K",zt)
2387	call write_output("teta","potential temperature","K",zh)
2388	call write_output("u","Zonal wind","m.s-1",zu)
2389	call write_output("v","Meridional wind","m.s-1",zv)
2390	call write_output("w","Vertical wind","m.s-1",pw)
2391	call write_output("p","Pressure","Pa",pplay)
2392	endif
2393
2394	! Total energy balance diagnostics
2395	if(callrad)then
2396	call write_output("ALB","Surface albedo"," ",albedo_equivalent)
2397	call write_output("ASR","absorbed stellar rad.","W m-2",fluxabs_sw)
2398	call write_output("ISR","incoming stellar rad.","W m-2",fluxtop_dn)
2399	call write_output("OLR","outgoing longwave rad.","W m-2",fluxtop_lw)
2400	call write_output("GND","heat flux from ground","W m-2",fluxgrd)
2401	call write_output("DYN","dynamical heat input","W m-2",fluxdyn)
2402	endif ! end of 'callrad'
2403
2404	if(enertest) then
2405	if (calldifv) then
2406	call write_output("q2","turbulent kinetic energy","J.kg^-1",q2)
2407	call write_output("sensibFlux","sensible heat flux","w.m^-2",sensibFlux)
2408	endif
2409
2410	if (corrk) then
2411	call write_output("dEzradsw","radiative heating","w.m^-2",dEzradsw)
2412	call write_output("dEzradlw","radiative heating","w.m^-2",dEzradlw)
2413	endif
2414
2415	if (generic_condensation) then
2416	call write_output("genericconddE","heat from generic condensation","W m-2",genericconddE)
2417	call write_output("dt_generic_condensation","heating from generic condensation","K s-1",dt_generic_condensation)
2418	endif
2419
2420	endif ! end of 'enertest'
2421
2422	! Diagnostics of optical thickness
2423	! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
2424	if (diagdtau) then
2425	do nw=1,L_NSPECTV
2426	write(str2,'(i2.2)') nw
2427	call write_output('dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',int_dtauv(:,nlayer:1:-1,nw))
2428	enddo
2429	do nw=1,L_NSPECTI
2430	write(str2,'(i2.2)') nw
2431	call write_output('dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',int_dtaui(:,nlayer:1:-1,nw))
2432	enddo
2433	endif
2434
2435	! Temporary inclusions for heating diagnostics.
2436	call write_output("zdtsw","SW heating","T s-1",zdtsw)
2437	call write_output("zdtlw","LW heating","T s-1",zdtlw)
2438	call write_output("dtrad","radiative heating","K s-1",dtrad)
2439	call write_output("zdtdyn","Dyn. heating","T s-1",zdtdyn)
2440
2441	! For Debugging.
2442	!call write_output('rnat','Terrain type',' ',real(rnat))
2443
2444	! Output tracers.
2445	if (tracer) then
2446	! call write_output("zdtc","tendancy T cond N2","K",zdtc)
2447
2448	do iq=1,nq
2449	call write_output(noms(iq),noms(iq),'kg/kg',zq(:,:,iq))
2450	! call write_output(trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2451	! 'kg m^-2',qsurf_hist(1,iq) )
2452	call write_output(trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2453	'kg m^-2',qcol(:,iq) )
2454	call write_output(trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2455	'kg m^-2',qsurf(:,iq) )
2456	enddo ! end of 'nq' loop
2457
2458	!Pluto specific
2459	call write_output('n2_iceflux','n2_iceflux',"kg m^-2 s^-1",flusurf(1,igcm_n2) )
2460	if (haze_radproffix)then
2461	call write_output('haze_reff','haze_reff','m',reffrad(1,1,1))
2462	end if
2463	if (methane) then
2464	! call write_output("rice_ch4","ch4 ice mass mean radius","m",rice_ch4)
2465	! call write_output("zq1temp_ch4"," "," ",zq1temp_ch4)
2466	! call write_output("qsat_ch4"," "," ",qsat_ch4)
2467	! call write_output("qsat_ch4_l1"," "," ",qsat_ch4_l1)
2468
2469	call write_output('ch4_iceflux','ch4_iceflux',&
2470	"kg m^-2 s^-1",flusurf(1,igcm_ch4_ice) )
2471	call write_output("vmr_ch4","vmr_ch4","%",vmr_ch4)
2472	if (.not.fast) then
2473	call write_output("zrho_ch4","zrho_ch4","kg.m-3",zrho_ch4(:,:))
2474	endif
2475
2476	! Tendancies
2477	call write_output("zdqcn2_ch4","zdq condn2 ch4","",&
2478	zdqc(:,:,igcm_ch4_gas))
2479	call write_output("zdqdif_ch4","zdqdif ch4","",&
2480	zdqdif(:,:,igcm_ch4_gas))
2481	call write_output("zdqsdif_ch4_ice","zdqsdif ch4","",&
2482	zdqsdif(:,igcm_ch4_ice))
2483	call write_output("zdqadj_ch4","zdqadj ch4","",&
2484	zdqadj(:,:,igcm_ch4_gas))
2485	if (sedimentation) then
2486	call write_output("zdqsed_ch4","zdqsed ch4","",&
2487	zdqsed(:,:,igcm_ch4_gas))
2488	call write_output("zdqssed_ch4","zdqssed ch4","",&
2489	zdqssed(:,igcm_ch4_gas))
2490	endif
2491	if (metcloud.and.(.not.fast)) then
2492	call write_output("zdtch4cloud","ch4 cloud","T s-1",&
2493	zdtch4cloud)
2494	call write_output("zdqch4cloud","ch4 cloud","T s-1",&
2495	zdqch4cloud(1,1,igcm_ch4_gas))
2496	endif
2497
2498	endif
2499
2500	if (carbox) then
2501	! call write_output("zdtcocloud","tendancy T cocloud","K",zdtcocloud)
2502	call write_output('co_iceflux','co_iceflux',&
2503	"kg m^-2 s^-1",flusurf(1,igcm_co_ice) )
2504	call write_output("vmr_co","vmr_co","%",vmr_co)
2505	if (.not.fast) THEN
2506	call write_output("zrho_co","zrho_co","kg.m-3",zrho_co(:,:))
2507	endif
2508	endif
2509
2510	if (haze) then
2511	! call write_output("zrho_haze","zrho_haze","kg.m-3",zrho_haze(:,:))
2512	call write_output("zdqrho_photprec","zdqrho_photprec",&
2513	"kg.m-3.s-1",zdqrho_photprec(:,:))
2514	call write_output("zdqphot_prec","zdqphot_prec","",&
2515	zdqphot_prec(:,:))
2516	call write_output("zdqhaze_ch4","zdqhaze_ch4","",&
2517	zdqhaze(:,:,igcm_ch4_gas))
2518	call write_output("zdqhaze_prec","zdqhaze_prec","",&
2519	zdqhaze(:,:,igcm_prec_haze))
2520	if (igcm_haze.ne.0) then
2521	call write_output("zdqhaze_haze","zdqhaze_haze","",&
2522	zdqhaze(:,:,igcm_haze))
2523	if (sedimentation) then
2524	call write_output("zdqssed_haze","zdqssed haze",&
2525	"kg/m2/s",zdqssed(:,igcm_haze))
2526	endif
2527	endif
2528	call write_output("zdqphot_ch4","zdqphot_ch4","",&
2529	zdqphot_ch4(:,:))
2530	call write_output("zdqconv_prec","zdqconv_prec","",&
2531	zdqconv_prec(:,:))
2532	! call write_output("zdqhaze_col","zdqhaze col","kg/m2/s",
2533	! & zdqhaze_col(:))
2534	endif
2535
2536	if (aerohaze) then
2537	call write_output("tau_col",&
2538	"Total aerosol optical depth","opacity",tau_col)
2539	endif
2540
2541	endif ! end of 'tracer'
2542
2543
2544	! Output spectrum.
2545	if(specOLR.and.corrk)then
2546	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
2547	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
2548	call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W/m^2/cm^-1",3,GSR_nu)
2549	endif
2550
2551	#else
2552	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
2553	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
2554	comm_ALBEQ(1:ngrid)=albedo_equivalent(1:ngrid)
2555	if (.not.calldifv) comm_LATENT_HF(:)=0.0
2556
2557	if ((tracer).and.(generic_condensation)) then
2558
2559	! If you have set generic_condensation (and not water) and you have set several GCS
2560	! then the outputs given to WRF will be only the ones for the last generic tracer
2561	! (Because it is rewritten every tracer in the loop)
2562	! WRF can take only one moist tracer
2563
2564	do iq=1,nq
2565	call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
2566
2567	if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
2568
2569	reffrad_generic_zeros_for_wrf(:,:) = 1.
2570
2571	comm_CLOUDFRAC(1:ngrid,1:nlayer) = cloudfrac(1:ngrid,1:nlayer)
2572	comm_TOTCLOUDFRAC(1:ngrid) = totcloudfrac(1:ngrid) !??????
2573	comm_SURFRAIN(1:ngrid) = zdqsrain_generic(1:ngrid,iq)
2574	comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_gas)
2575	comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_generic_ice)
2576	! comm_H2OICE_REFF(1:ngrid,1:nlayer) = reffrad_generic_zeros_for_wrf(1:ngrid,1:nlayer) ! for now zeros !
2577	!comm_H2OICE_REFF(1:ngrid,1:nlayer) = 0*zdtrain_generic(1:ngrid,1:nlayer) ! for now zeros !
2578	comm_REEVAP(1:ngrid) = reevap_precip_generic(1:ngrid,iq)
2579	comm_DTRAIN(1:ngrid,1:nlayer) = zdtrain_generic(1:ngrid,1:nlayer)
2580	comm_DTLSC(1:ngrid,1:nlayer) = dt_generic_condensation(1:ngrid,1:nlayer)
2581	comm_RH(1:ngrid,1:nlayer) = RH_generic(1:ngrid,1:nlayer,iq)
2582
2583	endif
2584	end do ! do iq=1,nq loop on tracers
2585
2586	else
2587	comm_CLOUDFRAC(1:ngrid,1:nlayer)=0.
2588	comm_TOTCLOUDFRAC(1:ngrid)=0.
2589	comm_SURFRAIN(1:ngrid)=0.
2590	comm_DQVAP(1:ngrid,1:nlayer)=0.
2591	comm_DQICE(1:ngrid,1:nlayer)=0.
2592	! comm_H2OICE_REFF(1:ngrid,1:nlayer)=0.
2593	comm_REEVAP(1:ngrid)=0.
2594	comm_DTRAIN(1:ngrid,1:nlayer)=0.
2595	comm_DTLSC(1:ngrid,1:nlayer)=0.
2596	comm_RH(1:ngrid,1:nlayer)=0.
2597
2598	endif ! if water, if generic_condensation, else
2599
2600	comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid)
2601	comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid)
2602	comm_FLUXTOP_LW(1:ngrid)=fluxtop_lw(1:ngrid)
2603	comm_FLUXSURF_SW(1:ngrid)=fluxsurf_sw(1:ngrid)
2604	comm_FLUXSURF_LW(1:ngrid)=fluxsurf_lw(1:ngrid)
2605	comm_FLXGRD(1:ngrid)=fluxgrd(1:ngrid)
2606	sensibFlux(1:ngrid) = zflubid(1:ngrid) - capcal(1:ngrid)*zdtsdif(1:ngrid) !!! ????
2607	comm_HR_DYN(1:ngrid,1:nlayer) = zdtdyn(1:ngrid,1:nlayer)
2608	#endif
2609
2610	! XIOS outputs
2611	#ifdef CPP_XIOS
2612	! Send fields to XIOS: (NB these fields must also be defined as
2613	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
2614	CALL send_xios_field("controle",tab_cntrl_mod,1)
2615
2616	CALL send_xios_field("ap",ap,1)
2617	CALL send_xios_field("bp",bp,1)
2618	CALL send_xios_field("aps",aps,1)
2619	CALL send_xios_field("bps",bps,1)
2620
2621	if (lastcall.and.is_omp_master) then
2622	write(,) "physiq: call xios_context_finalize"
2623	call xios_context_finalize
2624	endif
2625	#endif
2626
2627	if (check_physics_outputs) then
2628	! Check the validity of updated fields at the end of the physics step
2629	call check_physics_fields("end of physiq:", zt, zu, zv, pplev, zq)
2630	endif
2631
2632	icount=icount+1
2633
2634	end subroutine physiq
2635
2636	end module physiq_mod

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

Download in other formats:

Original Format