Context Navigation

physiq_mod.F90 @ 3196

Last change on this file since 3196 was 3196, checked in by afalco, 10 months ago

Pluto PCM:

Included N2 condensation from PLUTO.old & read haze aerosols.

Property svn:executable set to *

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

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

Download in other formats:

Original Format