Context Navigation

physiq_mod.F90 @ 3275

Last change on this file since 3275 was 3275, checked in by afalco, 8 months ago
Pluto PCM: Changed _vap to _gas; Included surfprop.F90; callcorrk includes methane AF
Property svn:executable set to ``*
File size: 100.4 KB

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

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

Download in other formats:

Original Format