Context Navigation

physiq_mod.F90 @ 3361

Last change on this file since 3361 was 3361, checked in by tbertrand, 5 months ago
LMDZ.PLUTO Fixing a bug in physiq.F when calling callcorrk_pluto (the old version of callcorrk for pluto) TB
Property svn:executable set to ``*
File size: 101.2 KB

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

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

Download in other formats:

Original Format