Context Navigation

physiq_mod.F90 @ 3356

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

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

Download in other formats:

Original Format