Context Navigation

physiq_mod.F90 @ 3390

Last change on this file since 3390 was 3390, checked in by tbertrand, 5 months ago
LMDZ.PLUTO resolving some issues in the code for 3D runs TB
Property svn:executable set to ``*
File size: 101.2 KB

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

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

Download in other formats:

Original Format