Context Navigation

physiq_mod.F90 @ 3237

Last change on this file since 3237 was 3237, checked in by afalco, 9 months ago
Pluto PCM: Import methane and ch4 clouds from pluto.old AF
Property svn:executable set to ``*
File size: 84.8 KB

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

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

Download in other formats:

Original Format