Context Navigation

physiq_mod.F90 @ 3642

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

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

Download in other formats:

Original Format