Context Navigation

physiq_mod.F90 @ 3691

Last change on this file since 3691 was 3688, checked in by afalco, 10 months ago
Pluto: xios: added missing fields. added explanation of newstart folder in deftank. 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, &
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	zdqc(:,:,:)=0.
521	zdqsc(:,:)=0.
522
523	! Initialize tracer names, indexes and properties.
524	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
525	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
526	if (tracer) then
527	call initracer(ngrid,nq)
528	! if(photochem) then !AF24: removed
529	endif
530	! Initialize aerosol indexes.
531	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
532	! call iniaerosol
533	! allocate related local arrays
534	! (need be allocated instead of automatic because of "naerkind")
535	allocate(dtau_aer(ngrid,nlayer,naerkind))
536
537	#ifdef CPP_XIOS
538	! Initialize XIOS context
539	write(,) "physiq: call wxios_context_init"
540	CALL wxios_context_init
541	#endif
542
543	! Read 'startfi.nc' file.
544	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
545	call phyetat0(startphy_file, &
546	ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
547	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,inertiedat)
548
549	if (.not.startphy_file) then
550	! starting without startfi.nc and with callsoil
551	! is not yet possible as soildepth default is not defined
552	if (callsoil) then
553	! default mlayer distribution, following a power law:
554	! mlayer(k)=lay1alpha*(k-1/2)
555	lay1=2.e-4
556	alpha=2
557	do iloop=0,nsoilmx-1
558	mlayer(iloop)=lay1(alpha*(iloop-0.5))
559	enddo
560	lay1=sqrt(mlayer(0)*mlayer(1))
561	alpha=mlayer(1)/mlayer(0)
562	do iloop=1,nsoilmx
563	layer(iloop)=lay1(alpha*(iloop-1))
564	enddo
565	endif
566	! additionnal "academic" initialization of physics
567	if (is_master) write(,) "Physiq: initializing tsurf(:) to pt(:,1) !!"
568	tsurf(:)=pt(:,1)
569	if (is_master) write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
570	do isoil=1,nsoilmx
571	tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
572	enddo
573	if (is_master) write(,) "Physiq: initializing day_ini to pday !"
574	day_ini=pday
575	endif
576
577	if (pday.ne.day_ini) then
578	write(,) "ERROR in physiq.F90:"
579	write(,) "bad synchronization between physics and dynamics"
580	write(,) "dynamics day: ",pday
581	write(,) "physics day: ",day_ini
582	stop
583	endif
584
585	write (,) 'In physiq day_ini =', day_ini
586
587	! Initialize albedo calculation.
588	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
589	albedo(:,:)=0.0
590	albedo_bareground(:)=0.0
591	albedo_snow_SPECTV(:)=0.0
592	albedo_n2_ice_SPECTV(:)=0.0
593
594	ptime0=ptime
595	write (,) 'In physiq ptime0 =', ptime
596
597	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV)
598
599	! Initialize orbital calculation.
600	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
601	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
602
603	savedeclin=0.
604	saveday=pday
605	adjust=0. ! albedo adjustment for convergeps
606
607	! Initialize soil.
608	! ~~~~~~~~~~~~~~~~
609	if (callsoil) then
610	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
611	ptimestep,tsurf,tsoil,capcal,fluxgrd)
612	else ! else of 'callsoil'.
613	print*,'WARNING! Thermal conduction in the soil turned off'
614	capcal(:)=1.e6
615	fluxgrd(:)=intheat
616	print*,'Flux from ground = ',intheat,' W m^-2'
617	endif ! end of 'callsoil'.
618
619	icount=1
620
621	! Initialize variable for dynamical heating and zonal wind tendency diagnostic
622	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
623	ztprevious(:,:)=pt(:,:)
624	zuprevious(:,:)=pu(:,:)
625
626	if(meanOLR)then
627	call system('rm -f rad_bal.out') ! to record global radiative balance.
628	call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
629	call system('rm -f h2o_bal.out') ! to record global hydrological balance.
630	endif
631
632	if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
633	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
634	ptimestep,pday+nday,time_phys,cell_area, &
635	albedo_bareground,zmea,zstd,zsig,zgam,zthe)
636	endif
637
638	! Initialize correlated-k.
639	! ~~~~~~~~~~~~~~~~~~~~~~~~
640	if (corrk) then
641	! We initialise the spectral grid here instead of
642	! at firstcall of callcorrk so we can output XspecIR, XspecVI
643	! when using Dynamico
644	if (is_master) print*, "physiq_mod: Correlated-k data base folder:",trim(datadir)
645	call getin_p("corrkdir",corrkdir)
646	if (is_master) print*,"corrkdir = ", corrkdir
647	write (tmp1, '(i4)') L_NSPECTI
648	write (tmp2, '(i4)') L_NSPECTV
649	banddir=trim(trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2)))
650	banddir=trim(trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir)))
651	call setspi !Basic infrared properties.
652	call setspv ! Basic visible properties.
653	call sugas_corrk ! Set up gaseous absorption properties.
654	if (optichaze) then
655	call suaer_corrk ! Set up aerosol optical properties.
656	endif
657	endif
658
659	! Initialize microphysics.
660	! ~~~~~~~~~~~~~~~~~~~~~~~~
661	IF (callmufi) THEN
662	! Initialize microphysics arrays.
663	call inimufi(ptimestep)
664	ENDIF ! end callmufi
665
666	!! call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
667	! XIOS outputs
668	#ifdef CPP_XIOS
669
670	if (is_master) write(,) "physiq: call initialize_xios_output"
671	call initialize_xios_output(pday,ptime,ptimestep,daysec,year_day, &
672	presnivs,pseudoalt,mlayer,WNOI,WNOV)
673	#endif
674
675	!! call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
676
677	write(,) "physiq: end of firstcall"
678	endif ! end of 'firstcall'
679
680	!! call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
681	!! call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
682
683	if (check_physics_inputs) then
684	!check the validity of input fields coming from the dynamics
685	call check_physics_fields("begin physiq:", pt, pu, pv, pplev, pq)
686	endif
687
688	! call writediagfi(ngrid,"pre_physical_rnat"," "," ",2,rnat)
689	! call writediagfi(ngrid,"pre_physical_capcal"," "," ",2,capcal)
690
691	! ------------------------------------------------------
692	! I.2 Initializations done at every physical timestep:
693	! ------------------------------------------------------
694
695	#ifdef CPP_XIOS
696	! update XIOS time/calendar
697	call update_xios_timestep
698	#endif
699
700	! Initialize various variables
701	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
702
703	pdt(1:ngrid,1:nlayer) = 0.0
704	zdtsurf(1:ngrid) = 0.0
705	pdq(1:ngrid,1:nlayer,1:nq) = 0.0
706	dqsurf(1:ngrid,1:nq)= 0.0
707	pdu(1:ngrid,1:nlayer) = 0.0
708	pdv(1:ngrid,1:nlayer) = 0.0
709	pdpsrf(1:ngrid) = 0.0
710	zflubid(1:ngrid) = 0.0
711	flux_sens_lat(1:ngrid) = 0.0
712	taux(1:ngrid) = 0.0
713	tauy(1:ngrid) = 0.0
714
715	if (conservn2) then
716	write(,) 'conservn2 iniloop'
717	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
718	endif
719
720	zday=pday+ptime ! Compute time, in sols (and fraction thereof).
721
722	! Compute Stellar Longitude (Ls), and orbital parameters.
723	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
724	if (season) then
725	call stellarlong(zday,zls)
726	else
727	call stellarlong(noseason_day,zls)
728	end if
729
730	! Get Lyman alpha flux at specific Ls
731	if (callmufi.or.haze) then
732	call lymalpha(zls,zfluxuv)
733	print*, 'Haze lyman-alpha zls,zfluxuv=',zls,zfluxuv
734	end if
735
736	IF (triton) then
737	CALL orbitetriton(zls,zday,dist_star,declin)
738	ELSE
739	call orbite(zls,dist_star,declin,right_ascen)
740	ENDIF
741
742	if (diurnal) then
743	zlss=-2.pi(zday-.5)
744	else if(diurnal .eqv. .false.) then
745	zlss=9999.
746	endif
747
748	glat(:) = g !AF24: removed oblateness
749
750	! Compute geopotential between layers.
751	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
752	zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
753	do l=1,nlayer
754	zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
755	enddo
756
757	zzlev(1:ngrid,1)=0.
758
759	do l=2,nlayer
760	do ig=1,ngrid
761	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
762	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
763	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
764	enddo
765	enddo
766
767	!Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
768
769	zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
770
771	! Compute potential temperature
772	! Note : Potential temperature calculation may not be the same in physiq and dynamic...
773	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
774	do l=1,nlayer
775	do ig=1,ngrid
776	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
777	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
778	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
779	massarea(ig,l)=mass(ig,l)*cell_area(ig)
780	enddo
781	enddo
782
783	! Compute vertical velocity (m/s) from vertical mass flux
784	! w = F / (rhoarea) and rho = P/(rT)
785	! But first linearly interpolate mass flux to mid-layers
786	if (.not.fast) then
787	do l=1,nlayer-1
788	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
789	enddo
790	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
791	do l=1,nlayer
792	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) / &
793	(pplay(1:ngrid,l)*cell_area(1:ngrid))
794	enddo
795	! omega in Pa/s
796	do l=1,nlayer-1
797	omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
798	enddo
799	omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
800	do l=1,nlayer
801	omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid)
802	enddo
803	endif
804
805	if (conservn2) then
806	write(,) 'conservn2 thermo'
807	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
808	endif
809
810	! Compute variations of g with latitude (to do).
811	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
812	gzlat(:,:) = g
813
814	! Initialize microphysical diagnostics.
815	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
816	IF (callmufi) THEN
817	! Initialize intensive to extensive factor (kg_air/m3: X/kg_air --> X/m3).
818	int2ext(:,:) = (pplev(:,1:nlayer)-pplev(:,2:nlayer+1)) / gzlat(:,1:nlayer) / (zzlev(:,2:nlayer+1)-zzlev(:,1:nlayer))
819
820	! Initialize microphysics diagnostics arrays.
821	call inimufi_diag(ngrid,nlayer,nq,pq,int2ext)
822	ENDIF ! end callmufi
823
824	! --------------------------------------------------------
825	! II.1 Thermosphere
826	! --------------------------------------------------------
827
828	! ajout de la conduction depuis la thermosphere
829	IF (callconduct) THEN
830
831	call conduction (ngrid,nlayer,ptimestep, &
832	pplay,pt,zzlay,zzlev,zdtconduc,tsurf)
833	DO l=1,nlayer
834	DO ig=1,ngrid
835	pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
836	ENDDO
837	ENDDO
838
839	ENDIF
840
841	! ajout de la viscosite moleculaire
842	IF (callmolvis) THEN
843	call molvis(ngrid,nlayer,ptimestep, &
844	pplay,pt,zzlay,zzlev, &
845	zdtconduc,pu,tsurf,zdumolvis)
846	call molvis(ngrid,nlayer,ptimestep, &
847	pplay,pt,zzlay,zzlev, &
848	zdtconduc,pv,tsurf,zdvmolvis)
849
850	DO l=1,nlayer
851	DO ig=1,ngrid
852	! pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
853	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
854	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
855	ENDDO
856	ENDDO
857	ENDIF
858
859	IF (callmoldiff) THEN
860	call moldiff_red(ngrid,nlayer,nq, &
861	pplay,pplev,pt,pdt,pq,pdq,ptimestep, &
862	zzlay,zdtconduc,zdqmoldiff)
863
864	DO l=1,nlayer
865	DO ig=1,ngrid
866	DO iq=1, nq
867	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
868	ENDDO
869	ENDDO
870	ENDDO
871	ENDIF
872
873	if (conservn2) then
874	write(,) 'conservn2 thermosphere'
875	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
876	endif
877
878
879	!---------------------------------
880	! II.2 Compute radiative tendencies
881	!---------------------------------
882	! Saving qsurf to compute paleo flux condensation/sublimation
883	DO iq=1, nq
884	DO ig=1,ngrid
885	IF (qsurf(ig,iq).lt.0.) then
886	qsurf(ig,iq)=0.
887	ENDIF
888	qsurf1(ig,iq)=qsurf(ig,iq)
889	ENDDO
890	ENDDO
891
892
893	! Compute local stellar zenith angles
894	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
895	fract = 0
896	if (diurnal) then
897	ztim1=SIN(declin)
898	ztim2=COS(declin)COS(2.pi*(zday-.5))
899	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
900
901	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
902	ztim1,ztim2,ztim3,mu0,fract)
903	else if(diurnal .eqv. .false.) then
904
905	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad)
906	! WARNING: this function appears not to work in 1D
907
908	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
909	mu0 = cos(pi*szangle/180.0)
910	fract= 1/(4*mu0) ! AF24: from pluto.old
911	endif
912
913	endif
914
915
916	! Pluto albedo /IT changes depending on surface ices (only in 3D)
917	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
918	if (ngrid.ne.1) then
919
920	!! Specific to change albedo of N2 so that Psurf
921	!! converges toward 1.4 Pa in "1989" seasons for Triton
922	!! converges toward 1.1 Pa in "2015" seasons for Pluto
923	if (convergeps) then
924	if (triton) then
925	! 1989 declination
926	if (declin180./pi.gt.-46..and.declin180./pi.lt.-45. &
927	.and.zday.gt.saveday+1000. &
928	.and.declin.lt.savedeclin) then
929	call planetwide_sumval(pplev(:,1)*cell_area(:)/totarea_planet,globave)
930	if (globave.gt.1.5) then
931	adjust=adjust+0.005
932	else if (globave.lt.1.3) then
933	adjust=adjust-0.005
934	endif
935	saveday=zday
936	endif
937	else
938	! Pluto : 2015 declination current epoch
939	if (declin180./pi.gt.50.and.declin180./pi.lt.51. &
940	.and.zday.gt.saveday+10000. &
941	.and.declin.gt.savedeclin) then
942	call planetwide_sumval(pplev(:,1)*cell_area(:)/totarea_planet,globave)
943	if (globave.gt.1.2) then
944	adjust=adjust+0.005
945	else if (globave.lt.1.) then
946	adjust=adjust-0.005
947	endif
948	saveday=zday
949	endif
950	endif
951	end if
952	end if ! if ngrid ne 1
953
954	call surfprop(ngrid,nq,fract,qsurf,tsurf, &
955	capcal,adjust,dist_star,flusurfold,ptimestep,zls,&
956	albedo,emis,therm_inertia)
957	! do i=2,ngrid
958	! albedo(i,:) = albedo(1,:)
959	! enddo
960	! AF24: TODO check albedo has been initialized here
961
962	if (firstcall.and.callsoil) then
963	! AF24 Originally in soil.F, but therm_inertia is modified by surfprop
964	! Additional checks: is the vertical discretization sufficient
965	! to resolve diurnal and annual waves?
966	do ig=1,ngrid
967	! extreme inertia for this column
968	inertia_min=minval(therm_inertia(ig,:))
969	inertia_max=maxval(therm_inertia(ig,:))
970	! diurnal and annual skin depth
971	diurnal_skin=(inertia_min/volcapa)*sqrt(daysec/pi)
972	annual_skin=(inertia_max/volcapa)sqrt(year_daydaysec/pi)
973	if (0.5*diurnal_skin<layer(1)) then
974	! one should have the fist layer be at least half of diurnal skin depth
975	write(,) "soil Error: grid point ig=",ig
976	write(,) " longitude=",longitude(ig)*(180./pi)
977	write(,) " latitude=",latitude(ig)*(180./pi)
978	write(,) " first soil layer depth ",layer(1)
979	write(,) " not small enough for a diurnal skin depth of ", &
980	diurnal_skin
981	write(,) " change soil layer distribution (comsoil_h.F90)"
982	call abort_physic("soil","change soil layer distribution (comsoil_h.F90)",1)
983	endif
984	if (2.*annual_skin>layer(nsoilmx)) then
985	! one should have the full soil be at least twice the diurnal skin depth
986	write(,) "soil Error: grid point ig=",ig
987	write(,) " longitude=",longitude(ig)*(180./pi)
988	write(,) " latitude=",latitude(ig)*(180./pi)
989	write(,) " total soil layer depth ",layer(nsoilmx)
990	write(,) " not large enough for an annual skin depth of ", &
991	annual_skin
992	write(,) " change soil layer distribution (comsoil_h.F90)"
993	call abort_physic("soil","change soil layer distribution (comsoil_h.F90)",1)
994	endif
995	enddo ! of do ig=1,ngrid
996
997	end if ! callsoil
998
999	if (callrad) then
1000	if( mod(icount-1,iradia).eq.0.or.lastcall) then
1001
1002	! Eclipse incoming sunlight !AF24: removed
1003
1004	!! call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1005	!! call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1006
1007
1008	if (corrk) then
1009
1010	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1011	! II.a Call correlated-k radiative transfer scheme
1012	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1013	if(kastprof)then
1014	print*,'kastprof should not = true here'
1015	call abort
1016	endif
1017
1018	! standard callcorrk
1019	if (oldplutocorrk) then
1020	call callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf, &
1021	albedo(:,1),emis,mu0,pplev,pplay,pt, &
1022	zzlay,zzlev,tsurf,fract,dist_star,dtau_aer, &
1023	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, &
1024	fluxtop_sw,fluxtop_dn,reffrad, &
1025	int_dtaui,int_dtauv,tau_col, &
1026	ptime,pday,firstcall,lastcall)
1027	albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
1028	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+ &
1029	fluxsurf_sw(1:ngrid)*(1.-albedo(1:ngrid,1))
1030	fluxabs_sw(1:ngrid)=fluxtop_dn(1:ngrid)-fluxtop_sw(1:ngrid)
1031	else
1032	muvar(1:ngrid,1:nlayer+1)=mugaz
1033	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
1034	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
1035	zzlay,zzlev,tsurf,fract,dist_star,dtau_aer,muvar, &
1036	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, &
1037	fluxsurfabs_sw,fluxtop_lw, &
1038	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,GSR_nu, &
1039	int_dtaui,int_dtauv, &
1040	tau_col,firstcall,lastcall)
1041	! Radiative flux from the sky absorbed by the surface (W.m-2).
1042	GSR=0.0
1043	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+ &
1044	fluxsurfabs_sw(1:ngrid)
1045	endif ! oldplutocorrk
1046	!GG (feb2021): Option to "artificially" decrease the raditive time scale in
1047	!the deep atmosphere press > 0.1 bar. Suggested by MT.
1048	!! COEFF_RAD to be "tuned" to facilitate convergence of tendency
1049
1050	!coeff_rad=0. ! 0 values, it doesn't accelerate the convergence
1051	!coeff_rad=0.5
1052	!coeff_rad=1.
1053	!do l=1, nlayer
1054	! do ig=1,ngrid
1055	! if(pplay(ig,l).ge.1.d4) then
1056	! zdtsw(ig,l)=zdtsw(ig,l)(pplay(ig,l)/1.d4)*coeff_rad
1057	! zdtlw(ig,l)=zdtlw(ig,l)(pplay(ig,l)/1.d4)*coeff_rad
1058	! endif
1059	! enddo
1060	!enddo
1061
1062	! AF24: removed CLFvarying Option
1063
1064
1065	!if(noradsurf)then ! no lower surface; SW flux just disappears
1066	! GSR = SUM(fluxsurf_sw(1:ngrid)*cell_area(1:ngrid))/totarea
1067	! fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
1068	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
1069	!endif
1070
1071	! Net atmospheric radiative heating rate (K.s-1)
1072	dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
1073
1074	! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that !
1075	if (firstcall .and. albedo_spectral_mode) then
1076	call spectral_albedo_calc(albedo_snow_SPECTV)
1077	endif
1078
1079	else
1080	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1081	! II.b Atmosphere has no radiative effect
1082	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1083	fluxtop_dn(1:ngrid) = fract(1:ngrid)mu0(1:ngrid)Fat1AU/dist_star**2
1084	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
1085	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
1086	endif
1087	fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
1088	print*,'------------WARNING---WARNING------------' ! by MT2015.
1089	print*,'You are in corrk=false mode, '
1090	print*,'and the surface albedo is taken equal to the first visible spectral value'
1091
1092	albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
1093	fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
1094	fluxabs_sw(1:ngrid)=fluxsurfabs_sw(1:ngrid)
1095	fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid)
1096	fluxtop_lw(1:ngrid) = emis(1:ngrid)sigmatsurf(1:ngrid)**4
1097
1098	dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating
1099
1100	endif ! end of corrk
1101
1102	endif ! of if(mod(icount-1,iradia).eq.0)
1103
1104
1105	! Transformation of the radiative tendencies
1106	! ------------------------------------------
1107	zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
1108	zplanck(1:ngrid)=emis(1:ngrid)sigmazplanck(1:ngrid)*zplanck(1:ngrid)
1109	fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
1110	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
1111
1112	! Test of energy conservation
1113	!----------------------------
1114	if(enertest)then
1115	call planetwide_sumval(cppmassarea(:,:)zdtsw(:,:)/totarea_planet,dEtotSW)
1116	call planetwide_sumval(cppmassarea(:,:)zdtlw(:,:)/totarea_planet,dEtotLW)
1117	!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
1118	call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
1119	call planetwide_sumval((fluxsurf_lw(:)emis(:)-zplanck(:))cell_area(:)/totarea_planet,dEtotsLW)
1120	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
1121	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
1122	if (is_master) then
1123	print*,'---------------------------------------------------------------'
1124	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
1125	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
1126	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
1127	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
1128	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
1129	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
1130	endif
1131	endif ! end of 'enertest'
1132
1133	endif ! of if (callrad)
1134
1135	!! call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1136	!! call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1137
1138	if (conservn2) then
1139	write(,) 'conservn2 radiat'
1140	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
1141	endif
1142
1143	! --------------------------------------------
1144	! III. Vertical diffusion (turbulent mixing) :
1145	! --------------------------------------------
1146
1147	if (calldifv) then
1148
1149	zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
1150
1151	if (oldplutovdifc) then
1152	zdum1(:,:) = 0.
1153	zdum2(:,:) = 0.
1154	zdh(:,:)=pdt(:,:)/zpopsk(:,:)
1155
1156	! Calling vdif (Martian version WITH N2 condensation)
1157	CALL vdifc_pluto(ngrid,nlayer,nq,zpopsk, &
1158	ptimestep,capcal,lwrite, &
1159	pplay,pplev,zzlay,zzlev,z0, &
1160	pu,pv,zh,pq,pt,tsurf,emis,qsurf, &
1161	zdum1,zdum2,zdh,pdq,pdt,zflubid, &
1162	zdudif,zdvdif,zdhdif,zdtsdif,q2, &
1163	zdqdif,zdqsdif,qsat_ch4,qsat_ch4_l1) !,zq1temp_ch4,qsat_ch4)
1164
1165	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
1166
1167	bcond=1./tcond1p4Pa
1168	acond=r/lw_n2
1169
1170	! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
1171	else if (UseTurbDiff) then
1172
1173	call turbdiff(ngrid,nlayer,nq, &
1174	ptimestep,capcal, &
1175	pplay,pplev,zzlay,zzlev,z0, &
1176	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
1177	pdt,pdq,zflubid, &
1178	zdudif,zdvdif,zdtdif,zdtsdif, &
1179	sensibFlux,q2,zdqdif,zdqevap,zdqsdif, &
1180	taux,tauy)
1181
1182	else ! if .not. (oldplutovdifc) .and. (UseTurbDiff)
1183
1184	zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1185
1186	call vdifc(ngrid,nlayer,nq,zpopsk, &
1187	ptimestep,capcal,lwrite, &
1188	pplay,pplev,zzlay,zzlev,z0, &
1189	pu,pv,zh,pq,tsurf,emis,qsurf, &
1190	zdh,pdq,zflubid, &
1191	zdudif,zdvdif,zdhdif,zdtsdif, &
1192	sensibFlux,q2,zdqdif,zdqsdif)
1193
1194	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1195	zdqevap(1:ngrid,1:nlayer)=0.
1196
1197	end if !end of 'UseTurbDiff'
1198
1199	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
1200
1201	!!! this is always done, except for turbulence-resolving simulations
1202	if (.not. turb_resolved) then
1203	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
1204	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
1205	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
1206	endif
1207
1208	if (tracer) then
1209	pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
1210	dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
1211	end if ! of if (tracer)
1212
1213	! test energy conservation
1214	!-------------------------
1215	if(enertest)then
1216
1217	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
1218	do ig = 1, ngrid
1219	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
1220	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
1221	enddo
1222
1223	call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
1224	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
1225	call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
1226	call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
1227
1228	if (is_master) then
1229
1230	if (UseTurbDiff) then
1231	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1232	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
1233	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1234	else
1235	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1236	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
1237	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1238	end if
1239	endif ! end of 'is_master'
1240
1241	! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
1242	endif ! end of 'enertest'
1243
1244	else ! calldifv
1245
1246	ztim1=4.sigmaptimestep
1247	DO ig=1,ngrid
1248	ztim2=ztim1emis(ig)tsurf(ig)**3
1249	z1=capcal(ig)*tsurf(ig)+ &
1250	ztim2tsurf(ig)+ (fluxrad(ig)+fluxgrd(ig))ptimestep
1251	z2= capcal(ig)+ztim2
1252	zdtsurf(ig)=(z1/z2 - tsurf(ig))/ptimestep
1253
1254	! for output:
1255	!dplanck(ig)=4.stephanptimestepemis(ig)tsurf(ig)**3
1256	ENDDO
1257
1258	! if(.not.newtonian)then
1259	!zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
1260
1261	! ------------------------------------------------------------------
1262	! Methane surface sublimation and condensation in fast model (nogcm)
1263	! ------------------------------------------------------------------
1264	if ((methane).and.(fast).and.condmetsurf) THEN
1265
1266	call ch4surf(ngrid,nlayer,nq,ptimestep,capcal, &
1267	tsurf,zdtsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf, &
1268	zdqch4fast,zdqsch4fast)
1269
1270	dqsurf(1:ngrid,igcm_ch4_ice)= dqsurf(1:ngrid,igcm_ch4_ice) + &
1271	zdqsch4fast(1:ngrid)
1272	pdq(1:ngrid,1,igcm_ch4_gas)= pdq(1:ngrid,1,igcm_ch4_gas) + &
1273	zdqch4fast(1:ngrid)
1274	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+lw_ch4*zdqsch4fast(1:ngrid)/capcal(1:ngrid)
1275	end if
1276	! ------------------------------------------------------------------
1277	! CO surface sublimation and condensation in fast model (nogcm)
1278	! ------------------------------------------------------------------
1279	if ((carbox).and.(fast).and.condcosurf) THEN
1280
1281	call cosurf(ngrid,nlayer,nq,ptimestep, &
1282	tsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf, &
1283	zdqcofast,zdqscofast)
1284
1285	dqsurf(1:ngrid,igcm_co_ice)= dqsurf(1:ngrid,igcm_co_ice) + &
1286	zdqscofast(1:ngrid)
1287	pdq(1:ngrid,1,igcm_co_gas)= pdq(1:ngrid,1,igcm_co_gas) + &
1288	zdqcofast(1:ngrid)
1289	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+lw_co*zdqscofast(1:ngrid)/capcal(1:ngrid)
1290	end if
1291
1292
1293	endif ! end of 'calldifv'
1294
1295	if (conservn2) then
1296	write(,) 'conservn2 calldifv'
1297	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)+ &
1298	dqsurf(:,1)*ptimestep)
1299	endif
1300	if (methane.and.conservch4) then
1301	write(,) 'conservch4 calldifv'
1302	if (fast) then
1303	call testconservfast(ngrid,nlayer,nq,pq(:,1,igcm_ch4_gas),pdq(:,1,igcm_ch4_gas), &
1304	qsurf(:,igcm_ch4_ice),dqsurf(:,igcm_ch4_ice), &
1305	ptimestep,pplev,zdqch4fast,zdqsch4fast,'CH4',' vdifc ')
1306	else
1307	call testconserv(ngrid,nlayer,nq,pq,pdq,qsurf,dqsurf, &
1308	igcm_ch4_gas,igcm_ch4_ice, &
1309	ptimestep,pplev,zdqdif,zdqsdif,'CH4',' vdifc ')
1310	endif
1311	endif
1312
1313	!-------------------
1314	! IV. Convection :
1315	!-------------------
1316
1317	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1318	! IV.a Dry convective adjustment :
1319	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1320
1321	if(calladj) then
1322
1323	zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1324	zduadj(1:ngrid,1:nlayer)=0.0
1325	zdvadj(1:ngrid,1:nlayer)=0.0
1326	zdhadj(1:ngrid,1:nlayer)=0.0
1327
1328
1329	call convadj(ngrid,nlayer,nq,ptimestep, &
1330	pplay,pplev,zpopsk, &
1331	pu,pv,zh,pq, &
1332	pdu,pdv,zdh,pdq, &
1333	zduadj,zdvadj,zdhadj, &
1334	zdqadj)
1335
1336	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
1337	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer)
1338	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
1339	zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1340
1341	if(tracer) then
1342	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
1343	end if
1344
1345	! Test energy conservation
1346	if(enertest)then
1347	call planetwide_sumval(cppmassarea(:,:)zdtadj(:,:)/totarea_planet,dEtot)
1348	if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
1349	endif
1350
1351	! ! Test water conservation !AF24: removed
1352
1353	endif ! end of 'calladj'
1354
1355	!-----------------------------------------------
1356	! V. Nitrogen condensation-sublimation :
1357	!-----------------------------------------------
1358
1359	if (n2cond) then
1360
1361	if (.not.tracer) then
1362	print*,'We need a N2 ice tracer to condense N2'
1363	call abort
1364	endif
1365
1366	call condense_n2(ngrid,nlayer,nq,ptimestep, &
1367	capcal,pplay,pplev,tsurf,pt, &
1368	pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, &
1369	qsurf(1,igcm_n2),albedo,emis, &
1370	zdtc,zdtsurfc,pdpsrf,zduc,zdvc, &
1371	zdqc,zdqsc(1,igcm_n2))
1372
1373	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
1374	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer)+zdvc(1:ngrid,1:nlayer)
1375	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer)+zduc(1:ngrid,1:nlayer)
1376	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
1377
1378	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq)
1379	dqsurf(1:ngrid,igcm_n2) = dqsurf(1:ngrid,igcm_n2) + zdqsc(1:ngrid,igcm_n2)
1380
1381	!! call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1382	!! call writediagfi(ngrid,"condense_n2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1383
1384	! test energy conservation
1385	if(enertest)then
1386	call planetwide_sumval(cppmassarea(:,:)zdtc(:,:)/totarea_planet,dEtot)
1387	call planetwide_sumval(capcal(:)zdtsurfc(:)cell_area(:)/totarea_planet,dEtots)
1388	if (is_master) then
1389	print*,'In n2cloud atmospheric energy change =',dEtot,' W m-2'
1390	print*,'In n2cloud surface energy change =',dEtots,' W m-2'
1391	endif
1392	endif
1393
1394	endif ! end of 'n2cond'
1395
1396	if (conservn2) then
1397	write(,) 'conservn2 n2cond'
1398	call testconservmass(ngrid,nlayer,pplev(:,1)+ &
1399	pdpsrf(:)ptimestep,qsurf(:,1)+dqsurf(:,1)ptimestep)
1400	endif
1401	if (methane.and.conservch4) then
1402	write(,) 'conservch4 n2cond'
1403	if (fast) then
1404	call testconservfast(ngrid,nlayer,nq,pq(:,1,igcm_ch4_gas),pdq(:,1,igcm_ch4_gas), &
1405	qsurf(:,igcm_ch4_ice),dqsurf(:,igcm_ch4_ice), &
1406	ptimestep,pplev,zdqch4fast,zdqsch4fast,'CH4',' n2cond')
1407	else
1408	call testconserv(ngrid,nlayer,nq,pq,pdq,qsurf,dqsurf, &
1409	igcm_ch4_gas,igcm_ch4_ice, &
1410	ptimestep,pplev,zdqc,zdqsc,'CH4',' n2cond')
1411	endif
1412	endif
1413
1414	!---------------------------------------------
1415	! VI. Specific parameterizations for tracers
1416	!---------------------------------------------
1417
1418	if (tracer) then
1419
1420	! ---------------------------------------
1421	! Methane ice condensation in the atmosphere
1422	! ----------------------------------------
1423	rice_ch4(:,:)=0 ! initialization needed for callsedim
1424	zdqch4cloud(:,:,:)=0.
1425	if ((methane).and.(metcloud).and.(.not.fast)) THEN
1426	call ch4cloud(ngrid,nlayer,ptimestep, &
1427	pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, &
1428	pq,pdq,zdqch4cloud,zdqsch4cloud,zdtch4cloud, &
1429	nq,rice_ch4)
1430
1431	DO l=1,nlayer
1432	DO ig=1,ngrid
1433	pdq(ig,l,igcm_ch4_gas)=pdq(ig,l,igcm_ch4_gas)+ &
1434	zdqch4cloud(ig,l,igcm_ch4_gas)
1435	pdq(ig,l,igcm_ch4_ice)=pdq(ig,l,igcm_ch4_ice)+ &
1436	zdqch4cloud(ig,l,igcm_ch4_ice)
1437	ENDDO
1438	ENDDO
1439
1440	! Increment methane ice surface tracer tendency
1441	DO ig=1,ngrid
1442	dqsurf(ig,igcm_ch4_ice)=dqsurf(ig,igcm_ch4_ice)+ &
1443	zdqsch4cloud(ig,igcm_ch4_ice)
1444	ENDDO
1445
1446	! update temperature tendancy
1447	DO ig=1,ngrid
1448	DO l=1,nlayer
1449	pdt(ig,l)=pdt(ig,l)+zdtch4cloud(ig,l)
1450	ENDDO
1451	ENDDO
1452	end if
1453
1454	! ---------------------------------------
1455	! CO ice condensation in the atmosphere
1456	! ----------------------------------------
1457	zdqcocloud(:,:,:)=0.
1458	IF ((carbox).and.(monoxcloud).and.(.not.fast)) THEN
1459	call cocloud(ngrid,nlayer,ptimestep, &
1460	pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, &
1461	pq,pdq,zdqcocloud,zdqscocloud,zdtcocloud, &
1462	nq,rice_co,qsurf(1,igcm_n2),dqsurf(1,igcm_n2))
1463
1464	DO l=1,nlayer
1465	DO ig=1,ngrid
1466	pdq(ig,l,igcm_co_gas)=pdq(ig,l,igcm_co_gas)+ &
1467	zdqcocloud(ig,l,igcm_co_gas)
1468	pdq(ig,l,igcm_co_ice)=pdq(ig,l,igcm_co_ice)+ &
1469	zdqcocloud(ig,l,igcm_co_ice)
1470	ENDDO
1471	ENDDO
1472
1473	! Increment CO ice surface tracer tendency
1474	DO ig=1,ngrid
1475	dqsurf(ig,igcm_co_ice)=dqsurf(ig,igcm_co_ice)+ &
1476	zdqscocloud(ig,igcm_co_ice)
1477	ENDDO
1478
1479	! update temperature tendancy
1480	DO ig=1,ngrid
1481	DO l=1,nlayer
1482	pdt(ig,l)=pdt(ig,l)+zdtcocloud(ig,l)
1483	ENDDO
1484	ENDDO
1485	ELSE
1486	rice_co(:,:)=0 ! initialization needed for callsedim
1487	END IF ! of IF (carbox)
1488
1489	! ----------------------------------------
1490	! VI.1. Microphysics / Aerosol particles
1491	! ----------------------------------------
1492	! Call of microphysics
1493	IF (callmufi) THEN
1494
1495	! Production for microphysics
1496	IF (call_haze_prod_pCH4) THEN
1497	zdqphot_prec(:,:) = 0.
1498	zdqphot_ch4(:,:) = 0.
1499	pdqmufi_prod(:,:,:) = 0.
1500	call hazecloud(ngrid,nlayer,nq,ptimestep,zday, &
1501	pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,pdqmufi_prod, &
1502	zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
1503	pdq(:,:,:) = pdq(:,:,:) + pdqmufi_prod(:,:,:) ! Should be updated
1504	ENDIF ! end call_haze_prod_pCH4
1505
1506	pdqmufi(:,:,:) = 0.
1507
1508	call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,pdqmufi_prod,pdqmufi)
1509
1510	pdq(:,:,:) = pdq(:,:,:) + pdqmufi(:,:,:)
1511
1512	ELSE
1513	IF (haze) THEN
1514	zdqphot_prec(:,:) = 0.
1515	zdqphot_ch4(:,:) = 0.
1516	zdqhaze(:,:,:) = 0.
1517	! Forcing to a fixed haze profile if haze_proffix
1518	if (haze_proffix.and.i_haze.gt.0.) then
1519	call haze_prof(ngrid,nlayer,zzlay,pplay,pt, &
1520	reffrad,profmmr)
1521	zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze))/ptimestep
1522	else
1523	call hazecloud(ngrid,nlayer,nq,ptimestep,zday, &
1524	pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,zdqhaze, &
1525	zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
1526	endif
1527	pdq(:,:,:) = pdq(:,:,:) + zdqhaze(:,:,:) ! Should be updated
1528	ENDIF ! end haze
1529
1530	IF (fast.and.fasthaze) THEN
1531	call prodhaze(ngrid,nlayer,nq,ptimestep,pplev,pq,pdq,dist_star, &
1532	mu0,declin,zdqprodhaze,zdqsprodhaze,gradflux,fluxbot, &
1533	fluxlym_sol_bot,fluxlym_ipm_bot,flym_sol,flym_ipm)
1534	DO ig=1,ngrid
1535	pdq(ig,1,igcm_ch4_gas)=pdq(ig,1,igcm_ch4_gas)+ &
1536	zdqprodhaze(ig,igcm_ch4_gas)
1537	pdq(ig,1,igcm_prec_haze)=pdq(ig,1,igcm_prec_haze)+ &
1538	zdqprodhaze(ig,igcm_prec_haze)
1539	pdq(ig,1,igcm_haze)=abs(pdq(ig,1,igcm_haze)+ &
1540	zdqprodhaze(ig,igcm_haze))
1541	qsurf(ig,igcm_haze)= qsurf(ig,igcm_haze)+ &
1542	zdqsprodhaze(ig)*ptimestep
1543	ENDDO
1544	ENDIF ! end fast.and.fasthaze
1545
1546	! Sedimentation.
1547	if (sedimentation) then
1548	zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
1549	zdqssed(1:ngrid,1:nq) = 0.0
1550	if (oldplutosedim)then
1551	call callsedim_pluto(ngrid,nlayer,ptimestep, &
1552	pplev,zzlev,pt,pdt,rice_ch4,rice_co, &
1553	pq,pdq,zdqsed,zdqssed,nq,pphi)
1554	else
1555	call callsedim(ngrid,nlayer,ptimestep, &
1556	pplev,zzlev,pt,pdt,pq,pdq, &
1557	zdqsed,zdqssed,nq,pphi)
1558	endif
1559	! Whether it falls as rain or snow depends only on the surface temperature
1560	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
1561	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
1562	end if ! end of 'sedimentation'
1563
1564	ENDIF ! end callmufi
1565
1566	! ---------------
1567	! VI.2. Updates
1568	! ---------------
1569
1570	! Updating Atmospheric Mass and Tracers budgets.
1571	if(mass_redistrib) then
1572
1573	zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * 0
1574	! ( zdqevap(1:ngrid,1:nlayer) &
1575	! ! + zdqrain(1:ngrid,1:nlayer,igcm_h2o_gas) &
1576	! ! + dqmoist(1:ngrid,1:nlayer,igcm_h2o_gas) &
1577	! + dqvaplscale(1:ngrid,1:nlayer) )
1578
1579	do ig = 1, ngrid
1580	zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
1581	enddo
1582
1583	! call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
1584	! call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
1585	call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
1586
1587	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
1588	capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
1589	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
1590	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
1591
1592	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
1593	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
1594	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
1595	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
1596	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
1597	pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
1598	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
1599
1600	endif
1601
1602	! call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1603
1604	! -----------------------------
1605	! VI.3. Surface Tracer Update
1606	! -----------------------------
1607
1608	qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
1609
1610	endif! end of if 'tracer'
1611
1612	if (conservn2) then
1613	write(,) 'conservn2 tracer'
1614	call testconservmass(ngrid,nlayer,pplev(:,1)+ &
1615	pdpsrf(:)*ptimestep,qsurf(:,1))
1616	endif
1617
1618	DO ig=1,ngrid
1619	flusurf(ig,igcm_n2)=(qsurf(ig,igcm_n2)- &
1620	qsurf1(ig,igcm_n2))/ptimestep
1621	flusurfold(ig,igcm_n2)=flusurf(ig,igcm_n2)
1622	if (methane) then
1623	flusurf(ig,igcm_ch4_ice)=(qsurf(ig,igcm_ch4_ice)- &
1624	qsurf1(ig,igcm_ch4_ice))/ptimestep
1625	flusurfold(ig,igcm_ch4_ice)=flusurf(ig,igcm_ch4_ice)
1626	endif
1627	if (carbox) then
1628	flusurf(ig,igcm_co_ice)=(qsurf(ig,igcm_co_ice)- &
1629	qsurf1(ig,igcm_co_ice))/ptimestep
1630	!flusurfold(ig,igcm_co_ice)=flusurf(ig,igcm_co_ice)
1631	endif
1632	ENDDO
1633
1634	!! Special source of haze particle !
1635	! todo: should be placed in haze and use tendency of n2 instead of flusurf
1636	IF (source_haze) THEN
1637	write(,) "Source haze not supported yet."
1638	stop
1639	! call hazesource(ngrid,nlayer,nq,ptimestep, &
1640	! pplev,flusurf,mu0,zdq_source)
1641
1642	DO iq=1, nq
1643	DO l=1,nlayer
1644	DO ig=1,ngrid
1645	pdq(ig,l,iq)=pdq(ig,l,iq)+zdq_source(ig,l,iq)
1646	ENDDO
1647	ENDDO
1648	ENDDO
1649	ENDIF
1650
1651	!------------------------------------------------
1652	! VII. Surface and sub-surface soil temperature
1653	!------------------------------------------------
1654
1655	! For diagnostic
1656	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1657	if (.not.fast) then
1658	DO ig=1,ngrid
1659	rho(ig,1) = pplay(ig,1)/(r*pt(ig,1))
1660	sensiblehf1(ig)=rho(ig,1)cpp(0.4/log(zzlay(ig,1)/z0))*2 &
1661	(pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1))*0.5 &
1662	(tsurf(ig)-pt(ig,1))
1663	if (calldifv) then
1664	sensiblehf2(ig)=zflubid(ig)-capcal(ig)*zdtsdif(ig)
1665	end if
1666	ENDDO
1667	endif
1668
1669
1670	! VII.1 Increment surface temperature
1671	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1672	tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
1673
1674	! Prevent surface (.e.g. non volatile ch4) to exceed max temperature
1675	! Lellouch et al., 2000,2011
1676	IF (tsurfmax) THEN
1677	DO ig=1,ngrid
1678	if (albedo_equivalent(ig).gt.albmin_ch4.and. &
1679	qsurf(ig,igcm_n2).eq.0.) then
1680	tsurf(ig)=min(tsurf(ig),54.)
1681	endif
1682	ENDDO
1683	ENDIF
1684
1685	! VII.2 Compute soil temperatures and subsurface heat flux.
1686	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1687	if (callsoil) then
1688	call soil(ngrid,nsoilmx,.false.,lastcall,therm_inertia, &
1689	ptimestep,tsurf,tsoil,capcal,fluxgrd)
1690	endif
1691
1692	! ! For output :
1693	! tidat_out(:,:)=0.
1694	! DO l=1,nsoilmx
1695	! tidat_out(:,l)=therm_inertia(:,l)
1696	! ENDDO
1697
1698	! Test energy conservation
1699	if(enertest)then
1700	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf(:)/totarea_planet,dEtots)
1701	if (is_master) print*,'Surface energy change =',dEtots,' W m-2'
1702	endif
1703
1704
1705
1706	! VII.3 multiply tendencies of cond/subli for paleo loop only in the
1707	! last Pluto year of the simulation
1708	! Year day must be adapted in the startfi for each object
1709	! Paleo uses year_day to calculate the annual mean tendancies
1710	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1711	IF (paleo) then
1712	if (zday.gt.day_ini+ptime0+nday-year_day) then
1713	DO iq=1,nq
1714	DO ig=1,ngrid
1715	qsurfyear(ig,iq)=qsurfyear(ig,iq)+ &
1716	(qsurf(ig,iq)-qsurf1(ig,iq)) !kg m-2 !ptimestep
1717	ENDDO
1718	ENDDO
1719	endif
1720	endif
1721
1722	! VII.4 Glacial flow at each timestep glastep or at lastcall
1723	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1724	IF (fast.and.glaflow) THEN
1725	if ((mod(zday-day_ini-ptime0,glastep)).lt.1. &
1726	.or.lastcall) then
1727	IF (lastcall) then
1728	dstep=mod(zday-day_ini-ptime0,glastep)*daysec
1729	else
1730	dstep=glastep*daysec
1731	endif
1732	zdqflow(:,:)=qsurf(:,:)
1733	IF (paleo) then
1734	call spreadglacier_paleo(ngrid,nq,qsurf, &
1735	phisfinew,dstep,tsurf)
1736	else
1737	call spreadglacier_simple(ngrid,nq,qsurf,dstep)
1738	endif
1739	zdqflow(:,:)=(zdqflow(:,:)-qsurf(:,:))/dstep
1740
1741	if (conservn2) then
1742	write(,) 'conservn2 glaflow'
1743	call testconservmass(ngrid,nlayer,pplev(:,1)+ &
1744	pdpsrf(:)*ptimestep,qsurf(:,1))
1745	endif
1746
1747	endif
1748	ENDIF
1749
1750	!---------------------------------------------------
1751	! VIII. Perform diagnostics and write output files
1752	!---------------------------------------------------
1753
1754	! Note : For output only: the actual model integration is performed in the dynamics.
1755
1756
1757	! Temperature, zonal and meridional winds.
1758	zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
1759	zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
1760	zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
1761
1762	!! Recast thermal plume vertical velocity array for outputs
1763	!! AF24: removed
1764
1765	! Diagnostic.
1766	zdtdyn(1:ngrid,1:nlayer) = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep
1767	ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
1768
1769	zdudyn(1:ngrid,1:nlayer) = (pu(1:ngrid,1:nlayer)-zuprevious(1:ngrid,1:nlayer)) / ptimestep
1770	zuprevious(1:ngrid,1:nlayer) = zu(1:ngrid,1:nlayer)
1771
1772	if(firstcall)then
1773	zdtdyn(1:ngrid,1:nlayer)=0.0
1774	zdudyn(1:ngrid,1:nlayer)=0.0
1775	endif
1776
1777	! Dynamical heating diagnostic
1778	fluxdyn(:)=0.0
1779	if (.not.fast) then
1780	do ig=1,ngrid
1781	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp
1782	enddo
1783	endif
1784
1785	! Tracers.
1786	zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
1787
1788	! Surface pressure.
1789	ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
1790	call planetwide_sumval(ps(:)*cell_area(:)/totarea_planet,globave)
1791
1792	! pressure density !pluto specific
1793	IF (.not.fast) then !
1794	do ig=1,ngrid
1795	do l=1,nlayer
1796	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1797	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1798	rho(ig,l) = zplay(ig,l)/(r*zt(ig,l))
1799	enddo
1800	zplev(ig,nlayer+1)=pplev(ig,nlayer+1)/pplev(ig,1)*ps(ig)
1801	enddo
1802	ENDIF
1803
1804
1805	! Surface and soil temperature information
1806	call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
1807	call planetwide_minval(tsurf(:),Ts2)
1808	call planetwide_maxval(tsurf(:),Ts3)
1809	if(callsoil)then
1810	TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
1811	if (is_master) then
1812	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
1813	print*,Ts1,Ts2,Ts3,TsS
1814	end if
1815	else
1816	if (is_master) then
1817	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
1818	print*,Ts1,Ts2,Ts3
1819	endif
1820	end if
1821
1822
1823	! Check the energy balance of the simulation during the run
1824	if(corrk)then
1825
1826	call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
1827	call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
1828	call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
1829	call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
1830	call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
1831	do ig=1,ngrid
1832	if(fluxtop_dn(ig).lt.0.0)then
1833	print*,'fluxtop_dn has gone crazy'
1834	print*,'fluxtop_dn=',fluxtop_dn(ig)
1835	print*,'tau_col=',tau_col(ig)
1836	print*,'dtau_aer=',dtau_aer(ig,:,:)
1837	print*,'temp= ',pt(ig,:)
1838	print*,'pplay= ',pplay(ig,:)
1839	call abort
1840	endif
1841	end do
1842
1843	if(ngrid.eq.1)then
1844	DYN=0.0
1845	endif
1846
1847	if (is_master) then
1848	print*,' ISR ASR OLR GND DYN [W m^-2]'
1849	print*, ISR,ASR,OLR,GND,DYN
1850	endif
1851
1852	if(enertest .and. is_master)then
1853	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
1854	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
1855	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
1856	endif
1857
1858	if(meanOLR .and. is_master)then
1859	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
1860	! to record global radiative balance
1861	open(92,file="rad_bal.out",form='formatted',position='append')
1862	write(92,*) zday,ISR,ASR,OLR
1863	close(92)
1864	open(93,file="tem_bal.out",form='formatted',position='append')
1865	if(callsoil)then
1866	write(93,*) zday,Ts1,Ts2,Ts3,TsS
1867	else
1868	write(93,*) zday,Ts1,Ts2,Ts3
1869	endif
1870	close(93)
1871	endif
1872	endif
1873
1874	endif ! end of 'corrk'
1875
1876
1877	! Diagnostic to test radiative-convective timescales in code.
1878	if(testradtimes)then
1879	open(38,file="tau_phys.out",form='formatted',position='append')
1880	ig=1
1881	do l=1,nlayer
1882	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1883	enddo
1884	close(38)
1885	print*,'As testradtimes enabled,'
1886	print*,'exiting physics on first call'
1887	call abort
1888	endif
1889
1890
1891	! Compute column amounts (kg m-2) if tracers are enabled.
1892	if(tracer)then
1893	qcol(1:ngrid,1:nq)=0.0
1894	do iq=1,nq
1895	do ig=1,ngrid
1896	qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
1897	enddo
1898	enddo
1899
1900	endif ! end of 'tracer'
1901
1902	if (methane) then
1903	IF (fast) then ! zq is the mixing ratio supposingly mixed in all atmosphere
1904	DO ig=1,ngrid
1905	vmr_ch4(ig)=zq(ig,1,igcm_ch4_gas)* &
1906	mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.
1907	ENDDO
1908	ELSE
1909	DO ig=1,ngrid
1910	! compute vmr methane
1911	vmr_ch4(ig)=qcol(ig,igcm_ch4_gas)* &
1912	g/ps(ig)mmol(igcm_n2)/mmol(igcm_ch4_gas)100.
1913	! compute density methane
1914	DO l=1,nlayer
1915	zrho_ch4(ig,l)=zq(ig,l,igcm_ch4_gas)*rho(ig,l)
1916	ENDDO
1917	ENDDO
1918	ENDIF
1919	endif
1920
1921	if (carbox) then
1922	IF (fast) then
1923	DO ig=1,ngrid
1924	vmr_co(ig)=zq(ig,1,igcm_co_gas)* &
1925	mmol(igcm_n2)/mmol(igcm_co_gas)*100.
1926	ENDDO
1927	ELSE
1928	DO ig=1,ngrid
1929	! compute vmr CO
1930	vmr_co(ig)=qcol(ig,igcm_co_gas)* &
1931	g/ps(ig)mmol(igcm_n2)/mmol(igcm_co_gas)100.
1932	! compute density CO
1933	DO l=1,nlayer
1934	zrho_co(ig,l)=zq(ig,l,igcm_co_gas)*rho(ig,l)
1935	ENDDO
1936	ENDDO
1937	ENDIF
1938	endif
1939
1940	zrho_haze(:,:)=0.
1941	zdqrho_photprec(:,:)=0.
1942	IF (haze.and.optichaze) then
1943	DO ig=1,ngrid
1944	DO l=1,nlayer
1945	zrho_haze(ig,l)=zq(ig,l,igcm_haze)*rho(ig,l)
1946	zdqrho_photprec(ig,l)=zdqphot_prec(ig,l)*rho(ig,l)
1947	ENDDO
1948	ENDDO
1949	ENDIF
1950
1951	IF (fasthaze) then
1952	DO ig=1,ngrid
1953	qcol(ig,igcm_haze)=zq(ig,1,igcm_haze)*pplev(ig,1)/g
1954	qcol(ig,igcm_prec_haze)=zq(ig,1,igcm_prec_haze)*pplev(ig,1)/g
1955	ENDDO
1956	ENDIF
1957
1958	! Info about Ls, declin...
1959	IF (fast) THEN
1960	if (is_master) write(,),'Ls=',zls180./pi,' dec=',declin180./pi
1961	if (is_master) write(,),'zday=',zday,' ps=',globave
1962	IF (lastcall) then
1963	if (is_master) write(,),'lastcall'
1964	ENDIF
1965	ELSE
1966	if (is_master) write(,),'Ls=',zls180./pi,'decli=',declin180./pi,'zday=',zday
1967	ENDIF
1968
1969	lecttsoil=0 ! default value for lecttsoil
1970	call getin_p("lecttsoil",lecttsoil)
1971	IF (lastcall.and.(ngrid.EQ.1).and.(lecttsoil.eq.1)) THEN
1972	! save tsoil temperature profile for 1D profile
1973	OPEN(13,file='proftsoil.out',form='formatted')
1974	DO i=1,nsoilmx
1975	write(13,*) tsoil(1,i)
1976	ENDDO
1977	CLOSE(13)
1978	ENDIF
1979
1980	if (is_master) print,'--> Ls =',zls180./pi
1981
1982	if(lastcall) then
1983	IF (grid_type==unstructured) THEN !IF DYNAMICO
1984	! DYNAMICO: no need to add a dynamics time step to ztime_fin
1985	ztime_fin = ptime
1986	ELSE ! LMDZ
1987	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1988	ENDIF ! of IF (grid_type==unstructured)
1989	!! Update surface ice distribution to iterate to steady state if requested
1990	!! AF24: removed
1991
1992	! endif
1993	if (paleo) then
1994	! time range for tendencies of ice flux qsurfyear
1995	zdt_tot=year_day ! Last year of simulation
1996
1997	masslost(:)=0.
1998	massacc(:)=0.
1999
2000	DO ig=1,ngrid
2001	! update new reservoir of ice on the surface
2002	DO iq=1,nq
2003	! kg/m2 to be sublimed or condensed during paleoyears
2004	qsurfyear(ig,iq)=qsurfyear(ig,iq)* &
2005	paleoyears365.25/(zdt_totdaysec/86400.)
2006
2007	! special case if we sublime the entire reservoir
2008	!! AF: TODO : fix following lines (real_area), using line below:
2009	! call planetwide_sumval((-qsurfyear(:,iq)-qsurf(:,iq))*cell_area(:),masslost)
2010
2011	! IF (-qsurfyear(ig,iq).gt.qsurf(ig,iq)) THEN
2012	! masslost(iq)=masslost(iq)+real_area(ig)* &
2013	! (-qsurfyear(ig,iq)-qsurf(ig,iq))
2014	! qsurfyear(ig,iq)=-qsurf(ig,iq)
2015	! ENDIF
2016
2017	! IF (qsurfyear(ig,iq).gt.0.) THEN
2018	! massacc(iq)=massacc(iq)+real_area(ig)*qsurfyear(ig,iq)
2019	! ENDIF
2020
2021
2022	ENDDO
2023	ENDDO
2024
2025	DO ig=1,ngrid
2026	DO iq=1,nq
2027	qsurfpal(ig,iq)=qsurf(ig,iq)+qsurfyear(ig,iq)
2028	IF (qsurfyear(ig,iq).gt.0.) THEN
2029	qsurfpal(ig,iq)=qsurfpal(ig,iq)- &
2030	qsurfyear(ig,iq)*masslost(iq)/massacc(iq)
2031	ENDIF
2032	ENDDO
2033	ENDDO
2034	! Finally ensure conservation of qsurf
2035	DO iq=1,nq
2036	call planetwide_sumval(qsurf(:,iq)*cell_area(:)/totarea_planet,globaveice(iq))
2037	call planetwide_sumval(qsurfpal(:,iq)*cell_area(:)/totarea_planet,globavenewice(iq))
2038	IF (globavenewice(iq).gt.0.) THEN
2039	qsurfpal(:,iq)=qsurfpal(:,iq)* &
2040	globaveice(iq)/globavenewice(iq)
2041	ENDIF
2042	ENDDO
2043
2044	! update new geopotential depending on the ice reservoir
2045	phisfipal(:)=phisfinew(:)+qsurfpal(:,igcm_n2)*g/1000.
2046	!phisfipal(ig)=phisfi(ig)
2047
2048	if (kbo.or.triton) then ! case of Triton : we do not change the orbital parameters
2049	pdaypal=pday ! no increment of pdaypal to keep same evolution of the subsolar point
2050	eccpal=1.-periastr/((periastr+apoastr)/2.) !no change of ecc
2051	peri_daypal=peri_day ! no change
2052	oblipal=obliquit ! no change
2053	tpalnew=tpal
2054	adjustnew=adjust
2055
2056	else ! Pluto
2057	! update new pday and tpal (Myr) to be set in startfi controle
2058	pdaypal=int(day_ini+paleoyears*365.25/6.3872)
2059	tpalnew=tpal+paleoyears*1.e-6 ! Myrs
2060
2061	! update new N2 ice adjustment (not tested yet on Pluto)
2062	adjustnew=adjust
2063
2064	! update milankovitch parameters : obliquity,Lsp,ecc
2065	call calcmilank(tpalnew,oblipal,peri_daypal,eccpal)
2066	!peri_daypal=peri_day
2067	!eccpal=0.009
2068	endif
2069
2070	if (is_master) write(,) "Paleo peri=",peri_daypal," tpal=",tpalnew
2071	if (is_master) write(,) "Paleo eccpal=",eccpal," tpal=",tpalnew
2072
2073	!----------------------------------------------------------------------
2074	! Writing NetCDF file "RESTARTFI" at the end of the run
2075	!----------------------------------------------------------------------
2076	! Note: 'restartfi' is stored just before dynamics are stored
2077	! in 'restart'. Between now and the writing of 'restart',
2078	! there will have been the itau=itau+1 instruction and
2079	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
2080	! thus we store for time=time+dtvr
2081
2082	! create restartfi
2083	if (ngrid.ne.1) then
2084	print*, "physdem1pal not yet implemented"
2085	stop
2086	!TODO: import this routine from pluto.old
2087	! call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq, &
2088	! ptimestep,pdaypal, &
2089	! ztime_fin,tsurf,tsoil,emis,q2,qsurfpal, &
2090	! cell_area,albedodat,therm_inertia,zmea,zstd,zsig, &
2091	! zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal, &
2092	! peri_daypal)
2093	endif
2094	else ! 'paleo'
2095
2096	if (ngrid.ne.1) then
2097	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
2098
2099	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
2100	ptimestep,ztime_fin,tsurf, &
2101	tsoil,therm_inertia,emis,albedo,q2,qsurf)
2102	endif
2103
2104	endif ! end of 'paleo'
2105	endif ! end of 'lastcall'
2106
2107	!------------------------------------------------------------------------------
2108	! OUTPUT in netcdf file "DIAGFI.NC",
2109	! containing any variable for diagnostic
2110	!
2111	! Note 1 : output with period "ecritphy", set in "run.def"
2112	! Note 2 : writediagfi can also be called from any other subroutine
2113	! for any variable, but its preferable to keep all the
2114	! calls in one place ...
2115	!------------------------------------------------------------------------------
2116
2117	!-------- General 1D variables
2118
2119	call write_output("Ls","solar longitude","deg",zls*180./pi)
2120	! call write_output("Lss","sub solar longitude","deg",zlss*180./pi)
2121	call write_output("RA","right ascension","deg",right_ascen*180./pi)
2122	call write_output("Declin","solar declination","deg",declin*180./pi)
2123	call write_output("dist_star","dist_star","AU",dist_star)
2124	call write_output("globave","surf press","Pa",globave)
2125
2126	!-------- General 2D variables
2127
2128	call write_output("tsurf","Surface temperature","K",tsurf)
2129	call write_output("ps","Surface pressure","Pa",ps)
2130	call write_output("emis","Emissivity","",emis)
2131	if (grid_type == regular_lonlat) then
2132	call write_output("area","Mesh area","m2", &
2133	cell_area_for_lonlat_outputs)
2134	else ! unstructured grid (e.g. dynamico)
2135	call write_output("area","Mesh area","m2",cell_area)
2136	endif
2137
2138	if (fast) then
2139	call write_output("fluxrad","fluxrad","W m-2",fluxrad)
2140	call write_output("fluxgrd","fluxgrd","W m-2",fluxgrd)
2141	! call write_output("dplanck","dplanck","W.s m-2 K-1",dplanck)
2142	! "soil" variables
2143	call write_output("capcal","capcal","W.s m-2 K-1",capcal)
2144	call write_output("tsoil","tsoil","K",tsoil)
2145	call write_output("therm_inertia","therm_inertia","S.I.",therm_inertia)
2146	endif
2147
2148	! Total energy balance diagnostics
2149	if(callrad)then
2150	call write_output("ALB","Surface albedo"," ",albedo_equivalent)
2151	call write_output("ASR","absorbed stellar rad.","W m-2",fluxabs_sw)
2152	call write_output("ISR","incoming stellar rad.","W m-2",fluxtop_dn)
2153	call write_output("OLR","outgoing longwave rad.","W m-2",fluxtop_lw)
2154	call write_output("GND","heat flux from ground","W m-2",fluxgrd)
2155	if (.not.fast) then
2156	call write_output("DYN","dynamical heat input","W m-2",fluxdyn)
2157	endif
2158	endif ! end of 'callrad'
2159
2160	!-------- General 3D variables
2161
2162	if (.not.fast) then
2163	if (check_physics_outputs) then
2164	! Check the validity of updated fields at the end of the physics step
2165	call check_physics_fields("physiq:", zt, zu, zv, pplev, zq)
2166	endif
2167
2168	call write_output("zzlay","Midlayer altitude", "m",zzlay(:,:))
2169	call write_output("zzlev","Interlayer altitude", "m",zzlev(:,1:nlayer))
2170	!call write_output('pphi','Geopotential',' ',pphi)
2171
2172	call write_output("temperature","temperature","K",zt)
2173	call write_output("teta","potential temperature","K",zh)
2174	call write_output("u","Zonal wind","m.s-1",zu)
2175	call write_output("v","Meridional wind","m.s-1",zv)
2176	call write_output("w","Vertical wind","m.s-1",pw)
2177	call write_output("p","Pressure","Pa",pplay)
2178	call write_output("omega","omega","Pa/s",omega)
2179	endif
2180
2181	if(enertest) then
2182	if (calldifv) then
2183	call write_output("q2","turbulent kinetic energy","J.kg^-1",q2)
2184	call write_output("sensibFlux","sensible heat flux","w.m^-2",sensibFlux)
2185	endif
2186	if (corrk) then
2187	call write_output("dEzradsw","radiative heating","w.m^-2",dEzRadsw)
2188	call write_output("dEzradlw","radiative heating","w.m^-2",dEzRadlw)
2189	endif
2190	endif ! end of 'enertest'
2191
2192	! Diagnostics of optical thickness
2193	! Warning this is exp(-tau), I let you postproc with -log to have tau itself
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
2203	! Temporary inclusions for heating diagnostics.
2204	if (.not.fast) then
2205	call write_output("zdtsw","SW heating","K s-1",zdtsw)
2206	call write_output("zdtlw","LW heating","K s-1",zdtlw)
2207	call write_output("dtrad","radiative heating","K s-1",dtrad)
2208	call write_output("zdtdyn","Dyn. heating","K s-1",zdtdyn)
2209	call write_output("zdudyn","Dyn. U","m s-2",zdudyn)
2210	call write_output("zdtconduc","tendancy conduc","K s-1",zdtconduc)
2211	call write_output("zdumolvis","tendancy molvis","m s-1",zdumolvis)
2212	call write_output("zdvmolvis","tendancy molvis","m s-1",zdvmolvis)
2213	call write_output("zdtdif","tendancy T diff","K s-1",zdtdif)
2214	call write_output("zdtsdif","tendancy Ts diff","K s-1",zdtsdif)
2215	call write_output("zdtadj","tendancy T adj","K s-1",zdtadj)
2216	endif
2217
2218	! Output tracers.
2219	if (tracer) then
2220
2221	do iq=1,nq
2222	if (.not.fast) then
2223	call write_output(noms(iq),noms(iq),'kg/kg',zq(:,:,iq))
2224	endif
2225	call write_output(trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2226	'kg m^-2',qcol(:,iq) )
2227	call write_output(trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2228	'kg m^-2',qsurf(:,iq) )
2229	enddo ! end of 'nq' loop
2230
2231	! N2 cycle
2232	call write_output('n2_iceflux','n2_iceflux',"kg m^-2 s^-1",flusurf(:,igcm_n2) )
2233	if (.not.fast) then
2234	call write_output("zdtc","tendancy T cond N2","K s-1",zdtc)
2235	call write_output("zdtsurfc","tendancy Ts cond N2","K s-1",zdtsurfc)
2236	call write_output("zduc","tendancy U cond N2","m s-1",zduc)
2237	call write_output("zdvc","tendancy V cond N2","m s-1",zdvc)
2238	call write_output("zdqc_n2","tendancy tracer cond N2","kg kg-1 s-1",zdqc(:,:,1))
2239	call write_output("zdqsc_n2","tendancy tracer surf cond N2","kg kg-1 s-1",zdqsc(:,1))
2240	call write_output("zdqdif_n2","tendancy tracer diff","kg kg-1 s-1",zdqdif(:,:,1))
2241	call write_output("zdqsdif_n2","tendancy tracer surf diff","kg kg-1 s-1",zdqsdif(:,1))
2242	call write_output("zdqadj_n2","tendancy tracer adj","K s-1",zdqadj(:,:,1))
2243	endif
2244
2245	! CH4 cycle
2246	if (methane) then
2247
2248	call write_output('ch4_iceflux','ch4_iceflux',&
2249	"kg m^-2 s^-1",flusurf(:,igcm_ch4_ice) )
2250	call write_output("vmr_ch4","vmr_ch4","%",vmr_ch4)
2251
2252	if (.not.fast) then
2253	call write_output("zrho_ch4","zrho_ch4","kg.m-3",zrho_ch4(:,:))
2254	!call write_output("rice_ch4","ch4 ice mass mean radius","m",rice_ch4)
2255	!call write_output("zq1temp_ch4"," "," ",zq1temp_ch4)
2256	!call write_output("qsat_ch4"," "," ",qsat_ch4)
2257	!call write_output("qsat_ch4_l1"," "," ",qsat_ch4_l1)
2258
2259	! 3D Tendancies
2260	call write_output("zdqcn2_ch4","zdq condn2 ch4","",&
2261	zdqc(:,:,igcm_ch4_gas))
2262	call write_output("zdqdif_ch4","zdqdif ch4","",&
2263	zdqdif(:,:,igcm_ch4_gas))
2264	call write_output("zdqsdif_ch4_ice","zdqsdif ch4","",&
2265	zdqsdif(:,igcm_ch4_ice))
2266	call write_output("zdqadj_ch4","zdqadj ch4","",&
2267	zdqadj(:,:,igcm_ch4_gas))
2268	endif
2269
2270	if (sedimentation) then
2271	call write_output("zdqsed_ch4","zdqsed ch4","",&
2272	zdqsed(:,:,igcm_ch4_gas))
2273	call write_output("zdqssed_ch4","zdqssed ch4","",&
2274	zdqssed(:,igcm_ch4_gas))
2275	endif
2276
2277	if (metcloud.and.(.not.fast)) then
2278	call write_output("zdtch4cloud","ch4 cloud","K s-1",&
2279	zdtch4cloud)
2280	call write_output("zdqch4cloud_gas","ch4 cloud","kg kg-1 s-1",&
2281	zdqch4cloud(:,:,igcm_ch4_gas))
2282	call write_output("zdqch4cloud_ice","ch4 cloud","kg kg-1 s-1",&
2283	zdqch4cloud(:,:,igcm_ch4_ice))
2284	endif
2285
2286	endif
2287
2288	! CO cycle
2289	if (carbox) then
2290	! call write_output("zdtcocloud","tendancy T cocloud","K",zdtcocloud)
2291	call write_output('co_iceflux','co_iceflux',&
2292	"kg m^-2 s^-1",flusurf(:,igcm_co_ice) )
2293	call write_output("vmr_co","vmr_co","%",vmr_co)
2294	if (.not.fast) THEN
2295	call write_output("zrho_co","zrho_co","kg.m-3",zrho_co(:,:))
2296	endif
2297	endif
2298
2299	! Haze
2300	if (haze) then
2301
2302	if (haze_radproffix)then
2303	call write_output('haze_reff','haze_reff','m',reffrad(:,:,1))
2304	end if
2305	!call write_output("zrho_haze","zrho_haze","kg.m-3",zrho_haze(:,:))
2306	!call write_output("zdqhaze_col","zdqhaze col","kg/m2/s",&
2307	! zdqhaze_col(:))
2308
2309	! 3D Tendencies
2310	call write_output("zdqrho_photprec","zdqrho_photprec",&
2311	"kg.m-3.s-1",zdqrho_photprec(:,:))
2312	call write_output("zdqphot_prec","zdqphot_prec","",&
2313	zdqphot_prec(:,:))
2314	call write_output("zdqhaze_ch4","zdqhaze_ch4","",&
2315	zdqhaze(:,:,igcm_ch4_gas))
2316	call write_output("zdqhaze_prec","zdqhaze_prec","",&
2317	zdqhaze(:,:,igcm_prec_haze))
2318	call write_output("zdqphot_ch4","zdqphot_ch4","",&
2319	zdqphot_ch4(:,:))
2320	call write_output("zdqconv_prec","zdqconv_prec","",&
2321	zdqconv_prec(:,:))
2322
2323	if (igcm_haze.ne.0) then
2324	call write_output("zdqhaze_haze","zdqhaze_haze","",&
2325	zdqhaze(:,:,igcm_haze))
2326	if (sedimentation) then
2327	call write_output("zdqssed_haze","zdqssed haze",&
2328	"kg/m2/s",zdqssed(:,igcm_haze))
2329	endif
2330	endif
2331
2332	if (optichaze) then
2333	call write_output("tau_col",&
2334	"Total aerosol optical depth","opacity",tau_col)
2335	endif
2336
2337	endif ! end haze
2338
2339	if (callmufi) then
2340	! Tracers:
2341	call write_output("m0as","Density number of spherical aerosols","m-3",zq(:,:,micro_indx(1))*int2ext(:,:))
2342	call write_output("m3as","Volume of spherical aerosols","m3.m-3",zq(:,:,micro_indx(2))*int2ext(:,:))
2343	call write_output("m0af","Density number of fractal aerosols","m-3",zq(:,:,micro_indx(3))*int2ext(:,:))
2344	call write_output("m3af","Volume of fractal aerosols","m3.m-3",zq(:,:,micro_indx(4))*int2ext(:,:))
2345
2346	! Diagnostics:
2347	call write_output("rcs","Characteristic radius of spherical aerosols","m",mp2m_rc_sph(:,:))
2348	call write_output("rcf","Characteristic radius of fractal aerosols","m",mp2m_rc_fra(:,:))
2349
2350	if (optichaze) then
2351	call write_output("tau_col",&
2352	"Total aerosol optical depth","opacity",tau_col)
2353	endif ! end optichaze
2354	endif ! end callmufi
2355
2356	endif ! end of 'tracer'
2357
2358	! Output spectrum.
2359	if(specOLR.and.corrk)then
2360	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
2361	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
2362	call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W/m^2/cm^-1",3,GSR_nu)
2363	endif
2364
2365	! XIOS outputs
2366	#ifdef CPP_XIOS
2367	! Send fields to XIOS: (NB these fields must also be defined as
2368	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
2369	CALL send_xios_field("controle",tab_cntrl_mod,1)
2370
2371	CALL send_xios_field("ap",ap,1)
2372	CALL send_xios_field("bp",bp,1)
2373	CALL send_xios_field("aps",aps,1)
2374	CALL send_xios_field("bps",bps,1)
2375
2376	if (lastcall.and.is_omp_master) then
2377	write(,) "physiq: call xios_context_finalize"
2378	call xios_context_finalize
2379	endif
2380	#endif
2381
2382	if (check_physics_outputs) then
2383	! Check the validity of updated fields at the end of the physics step
2384	call check_physics_fields("end of physiq:", zt, zu, zv, pplev, zq)
2385	endif
2386
2387	icount=icount+1
2388
2389	end subroutine physiq
2390
2391	end module physiq_mod

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

Download in other formats:

Original Format