Context Navigation

physiq_mod.F90 @ 3279

Last change on this file since 3279 was 3278, checked in by jleconte, 22 months ago
Correction on moistadj_generic tendencies (Noe C)
Property svn:executable set to ``*
File size: 121.6 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 watercommon_h, only : RLVTT, Psat_water,epsi,su_watercycle, RV, T_h2o_ice_liq
21	use generic_cloud_common_h, only : epsi_generic, Psat_generic
22	use thermcell_mod, only: init_thermcell_mod
23	use gases_h, only: gnom, gfrac, ngasmx
24	use radcommon_h, only: sigma, glat, grav, BWNV, WNOI, DWNI, DWNV, WNOV
25	use suaer_corrk_mod, only: suaer_corrk
26	use radii_mod, only: h2o_reffrad, co2_reffrad
27	use aerosol_mod, only: iniaerosol, iaero_co2, iaero_h2o
28	use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, &
29	dryness
30	use comdiurn_h, only: coslat, sinlat, coslon, sinlon
31	use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
32	use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
33	use geometry_mod, only: latitude, longitude, cell_area
34	USE comgeomfi_h, only: totarea, totarea_planet
35	USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
36	alpha_lift, alpha_devil, qextrhor, &
37	igcm_h2o_ice, igcm_h2o_vap, igcm_dustbin, &
38	igcm_co2_ice, nesp, is_chim, is_condensable,constants_epsi_generic
39	use time_phylmdz_mod, only: ecritphy, iphysiq, nday
40	use phyetat0_mod, only: phyetat0
41	use wstats_mod, only: callstats, wstats, mkstats
42	use phyredem, only: physdem0, physdem1
43	!! use slab_ice_h, only: capcalocean, capcalseaice,capcalsno, &
44	!! noceanmx
45	!! use ocean_slab_mod, only :ocean_slab_init, ocean_slab_ice, &
46	!! capcalocean, capcalseaice,capcalsno, nslay, &
47	!! ocean_slab_final
48	!! use surf_heat_transp_mod,only :init_masquv
49	use ocean_slab_mod, only :ocean_slab_init, ocean_slab_noice, ocean_slab_ice, &
50	ocean_slab_frac, ocean_slab_get_vars, &
51	capcalocean, capcalseaice,capcalsno, nslay, knon, &
52	ocean_slab_final
53	use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
54	use mod_phys_lmdz_para, only : is_master
55	use planete_mod, only: apoastr, periastr, year_day, peri_day, &
56	obliquit, nres, z0
57	use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
58	use time_phylmdz_mod, only: daysec
59	use callkeys_mod, only: albedo_spectral_mode, calladj, calldifv, &
60	calllott_nonoro, callrad, callsoil, nosurf, &
61	calltherm, CLFvarying, co2cond, corrk, diagdtau, &
62	diurnal, enertest, fat1au, flatten, j2, &
63	hydrology, intheat, iradia, kastprof, &
64	lwrite, mass_redistrib, massplanet, meanOLR, &
65	nearco2cond, newtonian, noseason_day, oblate, &
66	ok_slab_ocean, photochem, rings_shadow, rmean, &
67	season, sedimentation,generic_condensation, &
68	specOLR, callthermos, &
69	startphy_file, testradtimes, tlocked, &
70	tracer, UseTurbDiff, water, watercond, &
71	waterrain, generic_rain, global1d, szangle, &
72	moistadjustment, moistadjustment_generic, varspec, varspec_data, nvarlayer
73	use generic_tracer_index_mod, only: generic_tracer_index
74	use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
75	use check_fields_mod, only: check_physics_fields
76	use conc_mod, only: rnew, cpnew, ini_conc_mod
77	use phys_state_var_mod
78	use callcorrk_mod, only: callcorrk
79	use conduction_mod
80	use pindex_mod, only: pindex
81	use vdifc_mod, only: vdifc
82	use turbdiff_mod, only: turbdiff
83	use turb_mod, only : q2,sensibFlux,turb_resolved
84	use mass_redistribution_mod, only: mass_redistribution
85	use condensation_generic_mod, only: condensation_generic
86	use datafile_mod, only: datadir
87	#ifndef MESOSCALE
88	use vertical_layers_mod, only: presnivs, pseudoalt
89	use mod_phys_lmdz_omp_data, ONLY: is_omp_master
90	#else
91	use comm_wrf, only : comm_HR_SW, comm_HR_LW, &
92	comm_CLOUDFRAC,comm_TOTCLOUDFRAC, comm_RH, &
93	comm_SURFRAIN,comm_REEVAP,comm_HR_DYN, &
94	comm_DQICE,comm_DQVAP,comm_ALBEQ, &
95	comm_FLUXTOP_DN,comm_FLUXABS_SW, &
96	comm_FLUXTOP_LW,comm_FLUXSURF_SW, &
97	comm_FLUXSURF_LW,comm_FLXGRD, &
98	comm_DTRAIN,comm_DTLSC,comm_H2OICE_REFF, &
99	comm_LATENT_HF
100
101	#endif
102
103	#ifdef CPP_XIOS
104	use xios_output_mod, only: initialize_xios_output, &
105	update_xios_timestep, &
106	send_xios_field
107	use wxios, only: wxios_context_init, xios_context_finalize
108	#endif
109
110	implicit none
111
112
113	!==================================================================
114	!
115	! Purpose
116	! -------
117	! Central subroutine for all the physics parameterisations in the
118	! universal model. Originally adapted from the Mars LMDZ model.
119	!
120	! The model can be run without or with tracer transport
121	! depending on the value of "tracer" in file "callphys.def".
122	!
123	!
124	! It includes:
125	!
126	! I. Initialization :
127	! I.1 Firstcall initializations.
128	! I.2 Initialization for every call to physiq.
129	!
130	! II. Compute radiative transfer tendencies (longwave and shortwave) :
131	! II.a Option 1 : Call correlated-k radiative transfer scheme.
132	! II.b Option 2 : Call Newtonian cooling scheme.
133	! II.c Option 3 : Atmosphere has no radiative effect.
134	!
135	! III. Vertical diffusion (turbulent mixing) :
136	!
137	! IV. Convection :
138	! IV.a Thermal plume model
139	! IV.b Dry convective adjusment
140	!
141	! V. Condensation and sublimation of gases (currently just CO2).
142	!
143	! VI. Tracers
144	! VI.1. Water and water ice.
145	! VI.2 Photochemistry
146	! VI.3. Aerosols and particles.
147	! VI.4. Updates (pressure variations, surface budget).
148	! VI.5. Slab Ocean.
149	! VI.6. Surface Tracer Update.
150	!
151	! VII. Surface and sub-surface soil temperature.
152	!
153	! VIII. Perform diagnostics and write output files.
154	!
155	!
156	! arguments
157	! ---------
158	!
159	! INPUT
160	! -----
161	!
162	! ngrid Size of the horizontal grid.
163	! nlayer Number of vertical layers.
164	! nq Number of advected fields.
165	!
166	! firstcall True at the first call.
167	! lastcall True at the last call.
168	!
169	! pday Number of days counted from the North. Spring equinoxe.
170	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
171	! ptimestep timestep (s).
172	!
173	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa).
174	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
175	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2.s-2).
176	!
177	! pu(ngrid,nlayer) u, zonal component of the wind (ms-1).
178	! pv(ngrid,nlayer) v, meridional component of the wind (ms-1).
179	!
180	! pt(ngrid,nlayer) Temperature (K).
181	!
182	! pq(ngrid,nlayer,nq) Advected fields.
183	!
184	! pudyn(ngrid,nlayer) \
185	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
186	! ptdyn(ngrid,nlayer) / corresponding variables.
187	! pqdyn(ngrid,nlayer,nq) /
188	! flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
189	!
190	! OUTPUT
191	! ------
192	!
193	! pdu(ngrid,nlayer) \
194	! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding
195	! pdt(ngrid,nlayer) / variables due to physical processes.
196	! pdq(ngrid,nlayer) /
197	! pdpsrf(ngrid) /
198	!
199	!
200	! Authors
201	! -------
202	! Frederic Hourdin 15/10/93
203	! Francois Forget 1994
204	! Christophe Hourdin 02/1997
205	! Subroutine completely rewritten by F. Forget (01/2000)
206	! Water ice clouds: Franck Montmessin (update 06/2003)
207	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
208	! New correlated-k radiative scheme: R. Wordsworth (2009)
209	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
210	! Improved water cycle: R. Wordsworth / B. Charnay (2010)
211	! To F90: R. Wordsworth (2010)
212	! New turbulent diffusion scheme: J. Leconte (2012)
213	! Loops converted to F90 matrix format: J. Leconte (2012)
214	! No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
215	! Purge of the code : M. Turbet (2015)
216	! Photochemical core developped by F. Lefevre: B. Charnay (2017)
217	!==================================================================
218
219
220	! 0. Declarations :
221	! ------------------
222
223	include "netcdf.inc"
224
225	! Arguments :
226	! -----------
227
228	! INPUTS:
229	! -------
230
231	integer,intent(in) :: ngrid ! Number of atmospheric columns.
232	integer,intent(in) :: nlayer ! Number of atmospheric layers.
233	integer,intent(in) :: nq ! Number of tracers.
234
235	logical,intent(in) :: firstcall ! Signals first call to physics.
236	logical,intent(in) :: lastcall ! Signals last call to physics.
237
238	real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0.
239	real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
240	real,intent(in) :: ptimestep ! Physics timestep (s).
241	real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
242	real,intent(in) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
243	real,intent(in) :: pphi(ngrid,nlayer) ! Geopotential at mid-layer (m2s-2).
244	real,intent(in) :: pu(ngrid,nlayer) ! Zonal wind component (m/s).
245	real,intent(in) :: pv(ngrid,nlayer) ! Meridional wind component (m/s).
246	real,intent(in) :: pt(ngrid,nlayer) ! Temperature (K).
247	real,intent(in) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
248	real,intent(in) :: flxw(ngrid,nlayer) ! Vertical mass flux (ks/s) at lower boundary of layer
249
250	! OUTPUTS:
251	! --------
252
253	! Physical tendencies :
254
255	real,intent(out) :: pdu(ngrid,nlayer) ! Zonal wind tendencies (m/s/s).
256	real,intent(out) :: pdv(ngrid,nlayer) ! Meridional wind tendencies (m/s/s).
257	real,intent(out) :: pdt(ngrid,nlayer) ! Temperature tendencies (K/s).
258	real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
259	real,intent(out) :: pdpsrf(ngrid) ! Surface pressure tendency (Pa/s).
260
261	! Local saved variables:
262	! ----------------------
263	integer,save :: day_ini ! Initial date of the run (sol since Ls=0).
264	integer,save :: icount ! Counter of calls to physiq during the run.
265	!$OMP THREADPRIVATE(day_ini,icount)
266
267	! Local variables :
268	! -----------------
269
270	! Aerosol (dust or ice) extinction optical depth at reference wavelength
271	! for the "naerkind" optically active aerosols:
272
273	real,save,allocatable :: aerosol(:,:,:) ! Aerosols
274	!$OMP THREADPRIVATE(aerosol)
275	real zh(ngrid,nlayer) ! Potential temperature (K).
276	real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
277	real omega(ngrid,nlayer) ! omega velocity (Pa/s, >0 when downward)
278
279	integer l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_vap, igcm_generic_ice
280	logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
281
282	real zls ! Solar longitude (radians).
283	real zlss ! Sub solar point longitude (radians).
284	real zday ! Date (time since Ls=0, calculated in sols).
285	real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers.
286	real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries.
287
288	! VARIABLES for the thermal plume model
289
290	real f(ngrid) ! Mass flux norm
291	real fm(ngrid,nlayer+1) ! Mass flux
292	real fm_bis(ngrid,nlayer) ! Recasted fm
293	real entr(ngrid,nlayer) ! Entrainment
294	real detr(ngrid,nlayer) ! Detrainment
295	real dqevap(ngrid,nlayer,nq) ! water tracer mass mixing ratio variations due to evaporation
296	real dtevap(ngrid,nlayer) ! temperature variation due to evaporation
297	real zqtherm(ngrid,nlayer,nq) ! vapor mass mixing ratio after evaporation
298	real zttherm(ngrid,nlayer) ! temperature after evaporation
299	real fraca(ngrid,nlayer+1) ! Fraction of the surface that plumes occupies
300	real zw2(ngrid,nlayer+1) ! Vertical speed
301	real zw2_bis(ngrid,nlayer) ! Recasted zw2
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_co2 routine.
309	real zdtsdif(ngrid) ! Turbdiff/vdifc routines.
310	real zdtsurf_hyd(ngrid) ! Hydrol routine.
311
312	! For Thermosphere : (K/s)
313	real zdtconduc(ngrid,nlayer)
314
315	! For Atmospheric Temperatures : (K/s)
316	real dtlscale(ngrid,nlayer) ! Largescale routine.
317	real dt_generic_condensation(ngrid,nlayer) ! condensation_generic routine.
318	real zdtc(ngrid,nlayer) ! Condense_co2 routine.
319	real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
320	real zdttherm(ngrid,nlayer) ! Calltherm routine.
321	real zdtmr(ngrid,nlayer) ! Mass_redistribution routine.
322	real zdtrain(ngrid,nlayer) ! Rain routine.
323	real zdtrain_generic(ngrid,nlayer) ! Rain_generic routine.
324	real dtmoist(ngrid,nlayer) ! Moistadj routine.
325	real dt_ekman(ngrid,nslay), dt_hdiff(ngrid,nslay) ! Slab_ocean routine.
326	real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
327	real zdtchim(ngrid,nlayer) ! Calchim routine.
328
329	! For Surface Tracers : (kg/m2/s)
330	real dqsurf(ngrid,nq) ! Cumulated tendencies.
331	real zdqsurfc(ngrid) ! Condense_co2 routine.
332	real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines.
333	real zdqssed(ngrid,nq) ! Callsedim routine.
334	real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine.
335	real zdqsrain(ngrid), zdqssnow(ngrid) ! Rain routine.
336	real zdqsrain_generic(ngrid,nq), zdqssnow_generic(ngrid,nq) ! Rain_generic routine.
337	real dqs_hyd(ngrid,nq) ! Hydrol routine.
338	real reevap_precip(ngrid) ! re-evaporation flux of precipitation (integrated over the atmospheric column)
339	real reevap_precip_generic(ngrid,nq)
340
341	! For Tracers : (kg/kg_of_air/s)
342	real zdqc(ngrid,nlayer,nq) ! Condense_co2 routine.
343	real zdqadj(ngrid,nlayer,nq) ! Convadj routine.
344	real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines.
345	real zdqevap(ngrid,nlayer) ! Turbdiff routine.
346	real zdqtherm(ngrid,nlayer,nq) ! Calltherm routine.
347	real zdqsed(ngrid,nlayer,nq) ! Callsedim routine.
348	real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine.
349	real zdqrain(ngrid,nlayer,nq) ! Rain routine.
350	real dqmoist(ngrid,nlayer,nq) ! Moistadj routine.
351	real dqvaplscale(ngrid,nlayer) ! Largescale routine.
352	real dqcldlscale(ngrid,nlayer) ! Largescale routine.
353	real dqvaplscale_generic(ngrid,nlayer,nq) ! condensation_generic routine.
354	real dqcldlscale_generic(ngrid,nlayer,nq) ! condensation_generic routine.
355	real dq_rain_generic_vap(ngrid,nlayer,nq) ! rain_generic routine
356	real dq_rain_generic_cld(ngrid,nlayer,nq) ! rain_generic routine
357	REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine
358	REAL,allocatable,save :: zdqschim(:,:) ! Calchim_asis routine
359	!$OMP THREADPRIVATE(zdqchim,zdqschim)
360
361	REAL array_zero1(ngrid)
362	REAL array_zero2(ngrid,nlayer)
363
364	! For Winds : (m/s/s)
365	real zdvadj(ngrid,nlayer), zduadj(ngrid,nlayer) ! Convadj routine.
366	real zdutherm(ngrid,nlayer), zdvtherm(ngrid,nlayer) ! Calltherm routine.
367	real zdumr(ngrid,nlayer), zdvmr(ngrid,nlayer) ! Mass_redistribution routine.
368	real zdvdif(ngrid,nlayer), zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
369	real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
370	real zdhadj(ngrid,nlayer) ! Convadj routine.
371
372	! For Pressure and Mass :
373	real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
374	real zdmassmr_col(ngrid) ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
375	real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s).
376
377
378
379	! Local variables for LOCAL CALCULATIONS:
380	! ---------------------------------------
381	real zflubid(ngrid)
382	real zplanck(ngrid),zpopsk(ngrid,nlayer)
383	real ztim1,ztim2,ztim3, z1,z2
384	real ztime_fin
385	real zdh(ngrid,nlayer)
386	real gmplanet
387	real taux(ngrid),tauy(ngrid)
388
389
390	! local variables for DIAGNOSTICS : (diagfi & stat)
391	! -------------------------------------------------
392	real ps(ngrid) ! Surface Pressure.
393	real zt(ngrid,nlayer) ! Atmospheric Temperature.
394	real zu(ngrid,nlayer),zv(ngrid,nlayer) ! Zonal and Meridional Winds.
395	real zq(ngrid,nlayer,nq) ! Atmospheric Tracers.
396	real zdtadj(ngrid,nlayer) ! Convadj Diagnostic.
397	real zdtdyn(ngrid,nlayer) ! Dynamical Heating (K/s).
398	real zdudyn(ngrid,nlayer) ! Dynamical Zonal Wind tendency (m.s-2).
399
400	real reff(ngrid,nlayer) ! Effective dust radius (used if doubleq=T).
401	real vmr(ngrid,nlayer) ! volume mixing ratio
402	real time_phys
403
404	real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
405
406	real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
407
408	! included by RW for H2O Manabe scheme
409	real rneb_man(ngrid,nlayer) ! H2O cloud fraction (moistadj).
410	real rneb_lsc(ngrid,nlayer) ! H2O cloud fraction (large scale).
411
412
413	! to test energy conservation (RW+JL)
414	real mass(ngrid,nlayer),massarea(ngrid,nlayer)
415	real dEtot, dEtots, AtmToSurf_TurbFlux
416	real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
417	!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
418
419	!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
420
421	real dtmoist_max,dtmoist_min
422	real dItot, dItot_tmp, dVtot, dVtot_tmp
423
424
425	real h2otot ! Total Amount of water. For diagnostic.
426	real icesrf,liqsrf,icecol,vapcol ! Total Amounts of water (ice,liq,vap). For diagnostic.
427	real dWtot, dWtot_tmp, dWtots, dWtots_tmp
428	logical,save :: watertest
429	!$OMP THREADPRIVATE(watertest)
430
431	real qsat(ngrid,nlayer) ! Water Vapor Volume Mixing Ratio at saturation (kg/kg_of_air).
432	real RH(ngrid,nlayer) ! Relative humidity.
433	real psat_tmp
434
435	real qsat_generic(ngrid,nlayer,nq) ! generic condensable tracers (GCS) specific concentration at saturation (kg/kg_of_air).
436	real RH_generic(ngrid,nlayer,nq) ! generic condensable tracers (GCS) Relative humidity.
437	real rneb_generic(ngrid,nlayer,nq) ! GCS cloud fraction (generic condensation).
438	real psat_tmp_generic
439	real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
440	!$OMP THREADPRIVATE(metallicity)
441
442	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
443
444	logical clearsky ! For double radiative transfer call. By BC
445
446	! For Clear Sky Case.
447	real fluxsurfabs_sw1(ngrid) ! For SW/LW flux.
448	real tau_col1(ngrid) ! For aerosol optical depth diagnostic.
449	real OLR_nu1(ngrid,L_NSPECTI) ! Clear sky TOA LW radiation in each IR band
450	real OSR_nu1(ngrid,L_NSPECTV) ! Clear sky TOA SW radiation in each VI band
451	real GSR_nu1(ngrid,L_NSPECTV) ! Clear sky Surface SW radiation in each VI band
452	real int_dtaui1(ngrid,nlayer,L_NSPECTI),int_dtauv1(ngrid,nlayer,L_NSPECTV) ! For optical thickness diagnostics.
453	real tf, ntf
454
455	real net_fluxsurf_lw(ngrid) ! net longwave flux at the surface (for diagnostics only)
456
457	real nconsMAX, vdifcncons(ngrid), cadjncons(ngrid) ! Vdfic water conservation test. By RW
458
459	real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW
460
461	! For fixed variable molar mass
462	real, dimension(:),allocatable,save :: p_var
463	real, dimension(:),allocatable,save :: mu_var
464	real, dimension(:,:),allocatable,save :: frac_var
465	!$OMP THREADPRIVATE(p_var,mu_var,frac_var)
466	real :: mu_vari(nlayer)
467	real :: muvari(ngrid,nlayer)
468	integer ios,k
469	character*100 dt_file
470
471	real,save,allocatable :: reffcol(:,:)
472	!$OMP THREADPRIVATE(reffcol)
473
474	! Non-oro GW tendencies
475	REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
476	REAL d_t_hin(ngrid,nlayer)
477	! Diagnostics 2D of gw_nonoro
478	REAL zustrhi(ngrid), zvstrhi(ngrid)
479
480
481	real :: tsurf2(ngrid)
482	!! real :: flux_o(ngrid),flux_g(ngrid)
483	real :: flux_g(ngrid)
484	real :: flux_sens_lat(ngrid)
485	real :: qsurfint(ngrid,nq)
486	#ifdef MESOSCALE
487	REAL :: lsf_dt(nlayer)
488	REAL :: lsf_dq(nlayer)
489	#endif
490
491	! flags to trigger extra sanity checks
492	logical, save :: check_physics_inputs=.false.
493	logical, save :: check_physics_outputs=.false.
494	!$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)
495
496	! Misc
497	character*2 :: str2
498	character(len=10) :: tmp1
499	character(len=10) :: tmp2
500	!==================================================================================================
501
502	! -----------------
503	! I. INITIALISATION
504	! -----------------
505
506	! --------------------------------
507	! I.1 First Call Initialisation.
508	! --------------------------------
509	if (firstcall) then
510	call getin_p("check_physics_inputs", check_physics_inputs)
511	call getin_p("check_physics_outputs", check_physics_outputs)
512
513	! Allocate saved arrays (except for 1D model, where this has already
514	! been done)
515	#ifndef MESOSCALE
516	if (ngrid>1) call phys_state_var_init(nq)
517	#endif
518
519	! Variables set to 0
520	! ~~~~~~~~~~~~~~~~~~
521	dtrad(:,:) = 0.0
522	fluxrad(:) = 0.0
523	tau_col(:) = 0.0
524	zdtsw(:,:) = 0.0
525	zdtlw(:,:) = 0.0
526	zdtconduc(:,:)= 0.0
527
528	! Initialize fixed variable mu
529	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
530
531	if(varspec) then
532	IF (.NOT.ALLOCATED(p_var)) ALLOCATE(p_var(nvarlayer))
533	IF (.NOT.ALLOCATED(mu_var)) ALLOCATE(mu_var(nvarlayer))
534	IF (.NOT.ALLOCATED(frac_var)) ALLOCATE(frac_var(nvarlayer,ngasmx))
535	p_var = 0.0
536	mu_var = 0.0
537	frac_var = 0.0
538	muvari = 0.0
539	dt_file=TRIM(varspec_data)
540	open(33,file=dt_file,form='formatted',status='old',iostat=ios)
541	if (ios.ne.0) then ! file not found
542	write(,) 'Error from varspec'
543	write(,) 'Data file ',trim(varspec_data),' not found.'
544	write(,) 'Check that the data is in your run repository.'
545	call abort_physic !a verifier
546	else
547	do k=1,nvarlayer
548	read(33,*) p_var(k), mu_var(k),frac_var(k,1:ngasmx)
549	!The order of columns in frac_var must correspond to order of molecules gases.def
550	!The format of your file must be:
551	! pressure(k) molar_mass(k), molar_fraction_of_gas_1(k), molar_fraction_of_gas_2(k),..., molar_fraction_of_gas_ngasmx(k)
552	enddo
553	endif
554	close(33)
555	else
556	IF (.NOT.ALLOCATED(p_var)) ALLOCATE(p_var(nlayer))
557	IF (.NOT.ALLOCATED(mu_var)) ALLOCATE(mu_var(nlayer))
558	IF (.NOT.ALLOCATED(frac_var)) ALLOCATE(frac_var(nlayer,ngasmx))
559	endif
560
561	! Initialize tracer names, indexes and properties.
562	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
563	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
564	if (tracer) then
565	call initracer(ngrid,nq)
566	if(photochem) then
567	call ini_conc_mod(ngrid,nlayer)
568	IF (.NOT.ALLOCATED(zdqchim)) ALLOCATE(zdqchim(ngrid,nlayer,nesp))
569	IF (.NOT.ALLOCATED(zdqschim)) ALLOCATE(zdqschim(ngrid,nesp))
570	endif
571	endif
572	! Initialize aerosol indexes.
573	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
574	call iniaerosol
575	! allocate related local arrays
576	! (need be allocated instead of automatic because of "naerkind")
577	allocate(aerosol(ngrid,nlayer,naerkind))
578	allocate(reffcol(ngrid,naerkind))
579
580	#ifdef CPP_XIOS
581	! Initialize XIOS context
582	write(,) "physiq: call wxios_context_init"
583	CALL wxios_context_init
584	#endif
585
586	! Read 'startfi.nc' file.
587	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
588	#ifndef MESOSCALE
589	call phyetat0(startphy_file, &
590	ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
591	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, &
592	cloudfrac,totcloudfrac,hice, &
593	rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
594
595	#else
596
597	day_ini = pday
598	#endif
599
600	#ifndef MESOSCALE
601	if (.not.startphy_file) then
602	! additionnal "academic" initialization of physics
603	if (is_master) write(,) "Physiq: initializing tsurf(:) to pt(:,1) !!"
604	tsurf(:)=pt(:,1)
605	if (is_master) write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
606	do isoil=1,nsoilmx
607	tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
608	enddo
609	if (is_master) write(,) "Physiq: initializing day_ini to pday !"
610	day_ini=pday
611	endif
612	#endif
613	if (pday.ne.day_ini) then
614	write(,) "ERROR in physiq.F90:"
615	write(,) "bad synchronization between physics and dynamics"
616	write(,) "dynamics day: ",pday
617	write(,) "physics day: ",day_ini
618	stop
619	endif
620
621	write (,) 'In physiq day_ini =', day_ini
622
623	! Initialize albedo calculation.
624	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
625	albedo(:,:)=0.0
626	albedo_bareground(:)=0.0
627	albedo_snow_SPECTV(:)=0.0
628	albedo_co2_ice_SPECTV(:)=0.0
629	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
630
631	! Initialize orbital calculation.
632	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
633	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
634
635
636	if(tlocked)then
637	print*,'Planet is tidally locked at resonance n=',nres
638	print*,'Make sure you have the right rotation rate!!!'
639	endif
640
641	! Initialize oceanic variables.
642	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
643
644	if (ok_slab_ocean)then
645
646	!! call ocean_slab_init(ngrid,ptimestep, tslab, &
647	!! sea_ice, pctsrf_sic)
648
649	!! call ini_surf_heat_transp_mod()
650
651	knindex(:) = 0
652	knon = 0
653
654	do ig=1,ngrid
655	zmasq(ig)=1
656	! knindex(ig) = ig
657	if (nint(rnat(ig)).eq.0) then ! rnat = 0 (ocean), 1 (continent)
658	zmasq(ig)=0
659	knon=knon+1
660	knindex(knon)=ig
661	endif
662	enddo
663
664	call ocean_slab_init(ptimestep, pctsrf_sic, tslab, tsea_ice, sea_ice, zmasq)
665
666	!! CALL init_masquv(ngrid,zmasq)
667
668	endif ! end of 'ok_slab_ocean'.
669
670
671	! Initialize soil.
672	! ~~~~~~~~~~~~~~~~
673	if (callsoil) then
674
675	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
676	ptimestep,tsurf,tsoil,capcal,fluxgrd)
677
678	if (ok_slab_ocean) then
679	do ig=1,ngrid
680	if (nint(rnat(ig)).eq.2) then
681	capcal(ig)=capcalocean
682	if (pctsrf_sic(ig).gt.0.5) then
683	capcal(ig)=capcalseaice
684	if (qsurf(ig,igcm_h2o_ice).gt.0.) then
685	capcal(ig)=capcalsno
686	endif
687	endif
688	endif
689	enddo
690	endif ! end of 'ok_slab_ocean'.
691
692	else ! else of 'callsoil'.
693
694	print*,'WARNING! Thermal conduction in the soil turned off'
695	capcal(:)=1.e6
696	fluxgrd(:)=intheat
697	print*,'Flux from ground = ',intheat,' W m^-2'
698
699	endif ! end of 'callsoil'.
700
701	icount=1
702
703	! Here is defined the type of the surface : Continent or Ocean.
704	! BC2014 : This is now already done in newstart.
705	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
706	if (.not.ok_slab_ocean) then
707
708	rnat(:)=1.
709	if (.not.nosurf) then ! inertiedat only defined if there is a surface
710	do ig=1,ngrid
711	if(inertiedat(ig,1).gt.1.E4)then
712	rnat(ig)=0
713	endif
714	enddo
715
716	print*,'WARNING! Surface type currently decided by surface inertia'
717	print*,'This should be improved e.g. in newstart.F'
718	endif
719	endif ! end of 'ok_slab_ocean'.
720
721
722	! Initialize surface history variable.
723	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
724	qsurf_hist(:,:)=qsurf(:,:)
725
726	! Initialize variable for dynamical heating and zonal wind tendency diagnostic
727	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
728	ztprevious(:,:)=pt(:,:)
729	zuprevious(:,:)=pu(:,:)
730
731	! Set temperature just above condensation temperature (for Early Mars)
732	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
733	if(nearco2cond) then
734	write(,)' WARNING! Starting at Tcond+1K'
735	do l=1, nlayer
736	do ig=1,ngrid
737	pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4 &
738	-pt(ig,l)) / ptimestep
739	enddo
740	enddo
741	endif
742
743	if(meanOLR)then
744	call system('rm -f rad_bal.out') ! to record global radiative balance.
745	call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
746	call system('rm -f h2o_bal.out') ! to record global hydrological balance.
747	endif
748
749
750	watertest=.false.
751	if(water)then ! Initialize variables for water cycle
752
753	if(enertest)then
754	watertest = .true.
755	endif
756
757	endif
758
759	! Set metallicity for GCS
760	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
761	metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
762	call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic
763
764	! Set some parameters for the thermal plume model
765	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
766	if (calltherm) then
767	CALL init_thermcell_mod(g, rcp, r, pi, T_h2o_ice_liq, RV)
768	endif
769
770	call su_watercycle ! even if we don't have a water cycle, we might
771	! need epsi for the wvp definitions in callcorrk.F
772	! or RETV, RLvCp for the thermal plume model
773
774	#ifndef MESOSCALE
775	if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
776	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
777	ptimestep,pday+nday,time_phys,cell_area, &
778	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
779	endif
780
781	#endif
782	if (corrk) then
783	! We initialise the spectral grid here instead of
784	! at firstcall of callcorrk so we can output XspecIR, XspecVI
785	! when using Dynamico
786	print*, "physiq_mod: Correlated-k data base folder:",trim(datadir)
787	call getin_p("corrkdir",corrkdir)
788	print*,"corrkdir = ", corrkdir
789	write (tmp1, '(i4)') L_NSPECTI
790	write (tmp2, '(i4)') L_NSPECTV
791	banddir=trim(trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2)))
792	banddir=trim(trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir)))
793	call setspi !Basic infrared properties.
794	call setspv ! Basic visible properties.
795	call sugas_corrk ! Set up gaseous absorption properties.
796	call suaer_corrk ! Set up aerosol optical properties.
797	endif
798
799	! XIOS outputs
800	#ifdef CPP_XIOS
801
802	write(,) "physiq: call initialize_xios_output"
803	call initialize_xios_output(pday,ptime,ptimestep,daysec, &
804	year_day,presnivs,pseudoalt,WNOI,WNOV)
805	#endif
806
807	write(,) "physiq: end of firstcall"
808	endif ! end of 'firstcall'
809
810
811	if (check_physics_inputs) then
812	!check the validity of input fields coming from the dynamics
813	call check_physics_fields("begin physiq:", pt, pu, pv, pplev, pq)
814	endif
815
816
817	! ------------------------------------------------------
818	! I.2 Initializations done at every physical timestep:
819	! ------------------------------------------------------
820
821	#ifdef CPP_XIOS
822	! update XIOS time/calendar
823	call update_xios_timestep
824	#endif
825
826	! Initialize various variables
827	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
828
829	if ( .not.nearco2cond ) then
830	pdt(1:ngrid,1:nlayer) = 0.0
831	endif
832	zdtsurf(1:ngrid) = 0.0
833	pdq(1:ngrid,1:nlayer,1:nq) = 0.0
834	dqsurf(1:ngrid,1:nq)= 0.0
835	pdu(1:ngrid,1:nlayer) = 0.0
836	pdv(1:ngrid,1:nlayer) = 0.0
837	pdpsrf(1:ngrid) = 0.0
838	zflubid(1:ngrid) = 0.0
839	flux_sens_lat(1:ngrid) = 0.0
840	taux(1:ngrid) = 0.0
841	tauy(1:ngrid) = 0.0
842
843	zday=pday+ptime ! Compute time, in sols (and fraction thereof).
844
845	! Compute Stellar Longitude (Ls), and orbital parameters.
846	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
847	if (season) then
848	call stellarlong(zday,zls)
849	else
850	call stellarlong(noseason_day,zls)
851	end if
852
853	call orbite(zls,dist_star,declin,right_ascen)
854
855	if (tlocked) then
856	zlss=Mod(-(2.pi(zday/year_day)nres - right_ascen),2.pi)
857	elseif (diurnal) then
858	zlss=-2.pi(zday-.5)
859	else if(diurnal .eqv. .false.) then
860	zlss=9999.
861	endif
862
863
864	! Compute variations of g with latitude (oblate case).
865	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
866	if (oblate .eqv. .false.) then
867	glat(:) = g
868	else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
869	print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
870	call abort
871	else
872	gmplanet = MassPlanetgrav1e24
873	do ig=1,ngrid
874	glat(ig)= gmplanet/(Rmean*2) (1.D0 + 0.75 J2 - 2.0flatten/3. + (2.flatten - 15./4. J2) * cos(2. * (pi/2. - latitude(ig))))
875	end do
876	endif
877
878	! Compute geopotential between layers.
879	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
880	zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
881	do l=1,nlayer
882	zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
883	enddo
884
885	zzlev(1:ngrid,1)=0.
886
887	do l=2,nlayer
888	do ig=1,ngrid
889	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
890	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
891	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
892	enddo
893	enddo
894
895	!Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
896
897	zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
898
899	! Compute potential temperature
900	! Note : Potential temperature calculation may not be the same in physiq and dynamic...
901	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
902	do l=1,nlayer
903	do ig=1,ngrid
904	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
905	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
906	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
907	massarea(ig,l)=mass(ig,l)*cell_area(ig)
908	enddo
909	enddo
910
911	! Compute vertical velocity (m/s) from vertical mass flux
912	! w = F / (rhoarea) and rho = P/(rT)
913	! But first linearly interpolate mass flux to mid-layers
914	do l=1,nlayer-1
915	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
916	enddo
917	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
918	do l=1,nlayer
919	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) / &
920	(pplay(1:ngrid,l)*cell_area(1:ngrid))
921	enddo
922	! omega in Pa/s
923	do l=1,nlayer-1
924	omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
925	enddo
926	omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
927	do l=1,nlayer
928	omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid)
929	enddo
930
931	! ----------------------------------------------------------------
932	! Compute mean mass, cp, and R
933	! --------------------------------
934	#ifndef MESOSCALE
935	if(photochem) then
936	call concentrations(ngrid,nlayer,nq,pplay,pt,pdt,pq,pdq,ptimestep)
937	endif
938	#endif
939
940	!---------------------------------
941	! II. Compute radiative tendencies
942	!---------------------------------
943
944	! Compute local stellar zenith angles
945	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
946	if (tlocked) then
947	! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb
948	ztim1=SIN(declin)
949	ztim2=COS(declin)*COS(zlss)
950	ztim3=COS(declin)*SIN(zlss)
951
952	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
953	ztim1,ztim2,ztim3,mu0,fract, flatten)
954
955	elseif (diurnal) then
956	ztim1=SIN(declin)
957	ztim2=COS(declin)COS(2.pi*(zday-.5))
958	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
959
960	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
961	ztim1,ztim2,ztim3,mu0,fract, flatten)
962	else if(diurnal .eqv. .false.) then
963
964	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad,flatten)
965	! WARNING: this function appears not to work in 1D
966
967	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
968	mu0 = cos(pi*szangle/180.0)
969	endif
970
971	endif
972
973	if (callrad) then
974	if( mod(icount-1,iradia).eq.0.or.lastcall) then
975
976	! Eclipse incoming sunlight (e.g. Saturn ring shadowing).
977	if(rings_shadow) then
978	call call_rings(ngrid, ptime, pday, diurnal)
979	endif
980
981
982	if (corrk) then
983
984	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
985	! II.a Call correlated-k radiative transfer scheme
986	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
987	if(kastprof)then
988	print*,'kastprof should not = true here'
989	call abort
990	endif
991	if(water) then
992	! take into account water effect on mean molecular weight
993	muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap))
994	muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap))
995	elseif(generic_condensation) then
996	! take into account generic condensable specie (GCS) effect on mean molecular weight
997	do iq=1,nq
998
999	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
1000
1001	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
1002
1003	epsi_generic=constants_epsi_generic(iq)
1004
1005	muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap))
1006	muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap))
1007
1008	endif
1009	end do ! do iq=1,nq loop on tracers
1010	elseif(varspec) then
1011	!take into account fixed variable mean molecular weight
1012	do ig=1,ngrid
1013	call pindex(p_var,mu_var,pplay(ig,:),nvarlayer,nlayer,mu_vari)
1014	muvari(ig,:) = mu_vari
1015	enddo
1016
1017	muvar(1:ngrid,1:nlayer) = muvari(1:ngrid,1:nlayer)
1018	muvar(1:ngrid,nlayer+1) = muvari(1:ngrid,nlayer)
1019
1020	else
1021	muvar(1:ngrid,1:nlayer+1)=mugaz
1022	endif
1023
1024
1025	if(ok_slab_ocean) then
1026	tsurf2(:)=tsurf(:)
1027	do ig=1,ngrid
1028	if (nint(rnat(ig))==0) then
1029	tsurf(ig)=((1.-pctsrf_sic(ig))tslab(ig,1)4+pctsrf_sic(ig)tsea_ice(ig)4)0.25
1030	endif
1031	enddo
1032	endif !(ok_slab_ocean)
1033
1034	! standard callcorrk
1035	clearsky=.false.
1036	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
1037	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
1038	tsurf,fract,dist_star,aerosol,muvar, &
1039	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, &
1040	fluxsurfabs_sw,fluxtop_lw, &
1041	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,GSR_nu, &
1042	int_dtaui,int_dtauv, &
1043	tau_col,cloudfrac,totcloudfrac, &
1044	clearsky,p_var,frac_var,firstcall,lastcall)
1045
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	! Option to call scheme once more for clear regions
1063	if(CLFvarying)then
1064
1065	! ---> PROBLEMS WITH ALLOCATED ARRAYS : temporary solution in callcorrk: do not deallocate if CLFvarying ...
1066	clearsky=.true.
1067	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
1068	albedo,albedo_equivalent1,emis,mu0,pplev,pplay,pt, &
1069	tsurf,fract,dist_star,aerosol,muvar, &
1070	zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1, &
1071	fluxsurfabs_sw1,fluxtop_lw1, &
1072	fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,GSR_nu1, &
1073	int_dtaui1,int_dtauv1, &
1074	tau_col1,cloudfrac,totcloudfrac, &
1075	clearsky,p_var,frac_var,firstcall,lastcall)
1076	clearsky = .false. ! just in case.
1077
1078	! Sum the fluxes and heating rates from cloudy/clear cases
1079	do ig=1,ngrid
1080	tf=totcloudfrac(ig)
1081	ntf=1.-tf
1082	fluxsurf_lw(ig) = ntffluxsurf_lw1(ig) + tffluxsurf_lw(ig)
1083	fluxsurf_sw(ig) = ntffluxsurf_sw1(ig) + tffluxsurf_sw(ig)
1084	albedo_equivalent(ig) = ntfalbedo_equivalent1(ig) + tfalbedo_equivalent(ig)
1085	fluxsurfabs_sw(ig) = ntffluxsurfabs_sw1(ig) + tffluxsurfabs_sw(ig)
1086	fluxtop_lw(ig) = ntffluxtop_lw1(ig) + tffluxtop_lw(ig)
1087	fluxabs_sw(ig) = ntffluxabs_sw1(ig) + tffluxabs_sw(ig)
1088	tau_col(ig) = ntftau_col1(ig) + tftau_col(ig)
1089
1090	zdtlw(ig,1:nlayer) = ntfzdtlw1(ig,1:nlayer) + tfzdtlw(ig,1:nlayer)
1091	zdtsw(ig,1:nlayer) = ntfzdtsw1(ig,1:nlayer) + tfzdtsw(ig,1:nlayer)
1092
1093	OSR_nu(ig,1:L_NSPECTV) = ntfOSR_nu1(ig,1:L_NSPECTV) + tfOSR_nu(ig,1:L_NSPECTV)
1094	GSR_nu(ig,1:L_NSPECTV) = ntfGSR_nu1(ig,1:L_NSPECTV) + tfGSR_nu(ig,1:L_NSPECTV)
1095	OLR_nu(ig,1:L_NSPECTI) = ntfOLR_nu1(ig,1:L_NSPECTI) + tfOLR_nu(ig,1:L_NSPECTI)
1096	if (diagdtau) then
1097	int_dtauv(ig,:,1:L_NSPECTV) = ntfint_dtauv1(ig,:,1:L_NSPECTV) + tfint_dtauv(ig,:,1:L_NSPECTV)
1098	int_dtaui(ig,:,1:L_NSPECTI) = ntfint_dtaui1(ig,:,1:L_NSPECTI) + tfint_dtaui(ig,:,1:L_NSPECTI)
1099	endif
1100	enddo
1101
1102	endif ! end of CLFvarying.
1103
1104	if(ok_slab_ocean) then
1105	tsurf(:)=tsurf2(:)
1106	endif
1107
1108
1109	! Radiative flux from the sky absorbed by the surface (W.m-2).
1110	GSR=0.0
1111	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurfabs_sw(1:ngrid)
1112	net_fluxsurf_lw(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
1113
1114	!if(noradsurf)then ! no lower surface; SW flux just disappears
1115	! GSR = SUM(fluxsurf_sw(1:ngrid)*cell_area(1:ngrid))/totarea
1116	! fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
1117	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
1118	!endif
1119
1120	! Net atmospheric radiative heating rate (K.s-1)
1121	dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
1122
1123	! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that !
1124	if (firstcall .and. albedo_spectral_mode) then
1125	call spectral_albedo_calc(albedo_snow_SPECTV)
1126	endif
1127
1128	elseif(newtonian)then
1129
1130	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1131	! II.b Call Newtonian cooling scheme
1132	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1133	call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
1134
1135	zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep
1136	! e.g. surface becomes proxy for 1st atmospheric layer ?
1137
1138	else
1139
1140	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1141	! II.c Atmosphere has no radiative effect
1142	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1143	fluxtop_dn(1:ngrid) = fract(1:ngrid)mu0(1:ngrid)Fat1AU/dist_star**2
1144	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
1145	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
1146	endif
1147	fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
1148	print*,'------------WARNING---WARNING------------' ! by MT2015.
1149	print*,'You are in corrk=false mode, '
1150	print*,'and the surface albedo is taken equal to the first visible spectral value'
1151
1152	fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
1153	fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid)
1154	fluxtop_lw(1:ngrid) = emis(1:ngrid)sigmatsurf(1:ngrid)**4
1155
1156	dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating
1157
1158	endif ! end of corrk
1159
1160	endif ! of if(mod(icount-1,iradia).eq.0)
1161
1162
1163	! Transformation of the radiative tendencies
1164	! ------------------------------------------
1165	zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
1166	zplanck(1:ngrid)=emis(1:ngrid)sigmazplanck(1:ngrid)*zplanck(1:ngrid)
1167	net_fluxsurf_lw(1:ngrid)=net_fluxsurf_lw(1:ngrid)-zplanck(1:ngrid)
1168	fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
1169	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
1170
1171	! Test of energy conservation
1172	!----------------------------
1173	if(enertest)then
1174	call planetwide_sumval(cppmassarea(:,:)zdtsw(:,:)/totarea_planet,dEtotSW)
1175	call planetwide_sumval(cppmassarea(:,:)zdtlw(:,:)/totarea_planet,dEtotLW)
1176	!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
1177	call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
1178	call planetwide_sumval((fluxsurf_lw(:)emis(:)-zplanck(:))cell_area(:)/totarea_planet,dEtotsLW)
1179	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
1180	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
1181	if (is_master) then
1182	print*,'---------------------------------------------------------------'
1183	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
1184	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
1185	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
1186	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
1187	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
1188	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
1189	endif
1190	endif ! end of 'enertest'
1191
1192	endif ! of if (callrad)
1193
1194
1195	! --------------------------------------------
1196	! III. Vertical diffusion (turbulent mixing) :
1197	! --------------------------------------------
1198
1199	if (calldifv) then
1200
1201	zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
1202
1203	! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
1204	if (UseTurbDiff) then
1205
1206	call turbdiff(ngrid,nlayer,nq,rnat, &
1207	ptimestep,capcal, &
1208	pplay,pplev,zzlay,zzlev,z0, &
1209	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
1210	pdt,pdq,zflubid, &
1211	zdudif,zdvdif,zdtdif,zdtsdif, &
1212	sensibFlux,q2,zdqdif,zdqevap,zdqsdif, &
1213	taux,tauy)
1214
1215	else
1216
1217	zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1218
1219	call vdifc(ngrid,nlayer,nq,rnat,zpopsk, &
1220	ptimestep,capcal,lwrite, &
1221	pplay,pplev,zzlay,zzlev,z0, &
1222	pu,pv,zh,pq,tsurf,emis,qsurf, &
1223	zdh,pdq,zflubid, &
1224	zdudif,zdvdif,zdhdif,zdtsdif, &
1225	sensibFlux,q2,zdqdif,zdqsdif)
1226
1227	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1228	zdqevap(1:ngrid,1:nlayer)=0.
1229
1230	end if !end of 'UseTurbDiff'
1231
1232	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
1233
1234	!!! this is always done, except for turbulence-resolving simulations
1235	if (.not. turb_resolved) then
1236	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
1237	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
1238	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
1239	endif
1240
1241	if(ok_slab_ocean)then
1242	flux_sens_lat(1:ngrid)=(zdtsdif(1:ngrid)*capcal(1:ngrid)-fluxrad(1:ngrid))
1243	endif
1244
1245	if (tracer) then
1246	pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
1247	dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
1248	end if ! of if (tracer)
1249
1250
1251	! test energy conservation
1252	!-------------------------
1253	if(enertest)then
1254
1255	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
1256	do ig = 1, ngrid
1257	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
1258	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
1259	enddo
1260
1261	call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
1262	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
1263	call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
1264	call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
1265
1266	if (is_master) then
1267
1268	if (UseTurbDiff) then
1269	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1270	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
1271	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1272	else
1273	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1274	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
1275	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1276	end if
1277	endif ! end of 'is_master'
1278
1279	! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
1280	endif ! end of 'enertest'
1281
1282
1283	! Test water conservation.
1284	if(watertest.and.water)then
1285
1286	call planetwide_sumval(massarea(:,:)zdqdif(:,:,igcm_h2o_vap)ptimestep/totarea_planet,dWtot_tmp)
1287	call planetwide_sumval(zdqsdif(:,igcm_h2o_vap)cell_area(:)ptimestep/totarea_planet,dWtots_tmp)
1288	do ig = 1, ngrid
1289	vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
1290	enddo
1291	call planetwide_sumval(massarea(:,:)zdqdif(:,:,igcm_h2o_ice)ptimestep/totarea_planet,dWtot)
1292	call planetwide_sumval(zdqsdif(:,igcm_h2o_ice)cell_area(:)ptimestep/totarea_planet,dWtots)
1293	dWtot = dWtot + dWtot_tmp
1294	dWtots = dWtots + dWtots_tmp
1295	do ig = 1, ngrid
1296	vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
1297	enddo
1298	call planetwide_maxval(vdifcncons(:),nconsMAX)
1299
1300	if (is_master) then
1301	print*,'---------------------------------------------------------------'
1302	print*,'In difv atmospheric water change =',dWtot,' kg m-2'
1303	print*,'In difv surface water change =',dWtots,' kg m-2'
1304	print*,'In difv non-cons factor =',dWtot+dWtots,' kg m-2'
1305	print*,'In difv MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
1306	endif
1307
1308	endif ! end of 'watertest'
1309	!-------------------------
1310
1311	else ! calldifv
1312
1313	if(.not.newtonian)then
1314
1315	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
1316
1317	endif
1318
1319	endif ! end of 'calldifv'
1320
1321
1322	!-------------------
1323	! IV. Convection :
1324	!-------------------
1325
1326	! ~~~~~~~~~~~~~~~~~~~~~~~~~~
1327	! IV.a Thermal plume model :
1328	! ~~~~~~~~~~~~~~~~~~~~~~~~~~
1329
1330	IF (calltherm) THEN
1331
1332	! AB: We need to evaporate ice before calling thermcell_main.
1333	IF (water) THEN
1334	CALL evap(ngrid,nlayer,nq,ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,zqtherm,zttherm)
1335	ELSE
1336	zttherm(:,:) = pt(:,:) + pdt(:,:) * ptimestep
1337	zqtherm(:,:,:) = pq(:,:,:) + pdq(:,:,:) * ptimestep
1338	ENDIF
1339
1340	CALL thermcell_main(ngrid, nlayer, nq, ptimestep, firstcall, &
1341	pplay, pplev, pphi, zpopsk, &
1342	pu, pv, zttherm, zqtherm, &
1343	zdutherm, zdvtherm, zdttherm, zdqtherm, &
1344	fm, entr, detr, zw2, fraca)
1345
1346	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdutherm(1:ngrid,1:nlayer)
1347	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvtherm(1:ngrid,1:nlayer)
1348	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdttherm(1:ngrid,1:nlayer)
1349
1350	IF (tracer) THEN
1351	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqtherm(1:ngrid,1:nlayer,1:nq)
1352	ENDIF
1353
1354	ENDIF ! end of 'calltherm'
1355
1356	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1357	! IV.b Dry convective adjustment :
1358	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1359
1360	if(calladj) then
1361
1362	zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1363	zduadj(1:ngrid,1:nlayer)=0.0
1364	zdvadj(1:ngrid,1:nlayer)=0.0
1365	zdhadj(1:ngrid,1:nlayer)=0.0
1366
1367
1368	call convadj(ngrid,nlayer,nq,ptimestep, &
1369	pplay,pplev,zpopsk, &
1370	pu,pv,zh,pq, &
1371	pdu,pdv,zdh,pdq, &
1372	zduadj,zdvadj,zdhadj, &
1373	zdqadj)
1374
1375	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
1376	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer)
1377	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
1378	zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1379
1380	if(tracer) then
1381	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
1382	end if
1383
1384	! Test energy conservation
1385	if(enertest)then
1386	call planetwide_sumval(cppmassarea(:,:)zdtadj(:,:)/totarea_planet,dEtot)
1387	if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
1388	endif
1389
1390	! Test water conservation
1391	if(watertest)then
1392	call planetwide_sumval(massarea(:,:)zdqadj(:,:,igcm_h2o_vap)ptimestep/totarea_planet,dWtot_tmp)
1393	do ig = 1, ngrid
1394	cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
1395	enddo
1396	call planetwide_sumval(massarea(:,:)zdqadj(:,:,igcm_h2o_ice)ptimestep/totarea_planet,dWtot)
1397	dWtot = dWtot + dWtot_tmp
1398	do ig = 1, ngrid
1399	cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
1400	enddo
1401	call planetwide_maxval(cadjncons(:),nconsMAX)
1402
1403	if (is_master) then
1404	print*,'In convadj atmospheric water change =',dWtot,' kg m-2'
1405	print*,'In convadj MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
1406	endif
1407
1408	endif ! end of 'watertest'
1409
1410	endif ! end of 'calladj'
1411	!----------------------------------------------
1412	! Non orographic Gravity Waves:
1413	!---------------------------------------------
1414	IF (calllott_nonoro) THEN
1415
1416	CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep, &
1417	cpnew, rnew, &
1418	pplay, &
1419	zmax_th, &! max altitude reached by thermals (m)
1420	pt, pu, pv, &
1421	pdt, pdu, pdv, &
1422	zustrhi,zvstrhi, &
1423	d_t_hin, d_u_hin, d_v_hin)
1424
1425	! Update tendencies
1426	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) &
1427	+ d_t_hin(1:ngrid,1:nlayer)
1428	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) &
1429	+ d_u_hin(1:ngrid,1:nlayer)
1430	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) &
1431	+ d_v_hin(1:ngrid,1:nlayer)
1432	print*,'du_nonoro: max ', maxval(d_u_hin), 'min ', minval(d_u_hin)
1433	print*,'dv_nonoro: max ', maxval(d_v_hin), 'min ', minval(d_v_hin)
1434
1435	ENDIF ! of IF (calllott_nonoro)
1436
1437
1438
1439	!-----------------------------------------------
1440	! V. Carbon dioxide condensation-sublimation :
1441	!-----------------------------------------------
1442
1443	if (co2cond) then
1444
1445	if (.not.tracer) then
1446	print*,'We need a CO2 ice tracer to condense CO2'
1447	call abort
1448	endif
1449	call condense_co2(ngrid,nlayer,nq,ptimestep, &
1450	capcal,pplay,pplev,tsurf,pt, &
1451	pdt,zdtsurf,pq,pdq, &
1452	qsurf,zdqsurfc,albedo,emis, &
1453	albedo_bareground,albedo_co2_ice_SPECTV, &
1454	zdtc,zdtsurfc,pdpsrf,zdqc)
1455
1456	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
1457	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
1458
1459	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq)
1460	dqsurf(1:ngrid,igcm_co2_ice) = dqsurf(1:ngrid,igcm_co2_ice) + zdqsurfc(1:ngrid)
1461
1462
1463	! test energy conservation
1464	if(enertest)then
1465	call planetwide_sumval(cppmassarea(:,:)zdtc(:,:)/totarea_planet,dEtot)
1466	call planetwide_sumval(capcal(:)zdtsurfc(:)cell_area(:)/totarea_planet,dEtots)
1467	if (is_master) then
1468	print*,'In co2cloud atmospheric energy change =',dEtot,' W m-2'
1469	print*,'In co2cloud surface energy change =',dEtots,' W m-2'
1470	endif
1471	endif
1472
1473	endif ! end of 'co2cond'
1474
1475
1476	!---------------------------------------------
1477	! VI. Specific parameterizations for tracers
1478	!---------------------------------------------
1479
1480	if (tracer) then
1481
1482	! ---------------------
1483	! VI.1. Water and ice
1484	! ---------------------
1485	if (water) then
1486
1487	! Water ice condensation in the atmosphere
1488	if(watercond.and.(RLVTT.gt.1.e-8))then
1489
1490	if ((.not.calltherm).and.moistadjustment) then
1491	dqmoist(1:ngrid,1:nlayer,1:nq)=0.
1492	dtmoist(1:ngrid,1:nlayer)=0.
1493
1494	! Moist Convective Adjustment.
1495	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1496
1497	call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
1498
1499	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) &
1500	+ dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)
1501	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) &
1502	+ dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice)
1503	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer)
1504
1505	! Test energy conservation.
1506	if(enertest)then
1507	call planetwide_sumval(cppmassarea(:,:)dtmoist(:,:)/totarea_planet,dEtot)
1508	call planetwide_maxval(dtmoist(:,:),dtmoist_max)
1509	call planetwide_minval(dtmoist(:,:),dtmoist_min)
1510	madjdEz(:,:)=cppmass(:,:)dtmoist(:,:)
1511
1512	do ig=1,ngrid
1513	madjdE(ig) = cppSUM(mass(:,:)dtmoist(:,:))
1514	enddo
1515
1516	if (is_master) then
1517	print*,'In moistadj atmospheric energy change =',dEtot,' W m-2'
1518	print,'In moistadj MAX atmospheric energy change =',dtmoist_maxptimestep,'K/step'
1519	print,'In moistadj MIN atmospheric energy change =',dtmoist_minptimestep,'K/step'
1520	endif
1521
1522	call planetwide_sumval(massarea(:,:)dqmoist(:,:,igcm_h2o_vap)ptimestep/totarea_planet+ &
1523	massarea(:,:)dqmoist(:,:,igcm_h2o_ice)ptimestep/totarea_planet,dWtot)
1524	if (is_master) print*,'In moistadj atmospheric water change =',dWtot,' kg m-2'
1525
1526	endif ! end of 'enertest'
1527	else
1528	! rneb_man, dqmoist are output of moistadj, but used later on
1529	! so set them to zero if moistadj is not called
1530	rneb_man(:,:)=0
1531	dqmoist(:,:,:)=0
1532	endif ! end of '(.not.calltherm).and.moistadjustment'
1533
1534	! Large scale condensation/evaporation.
1535	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1536	call largescale(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
1537
1538	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtlscale(1:ngrid,1:nlayer)
1539	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)+dqvaplscale(1:ngrid,1:nlayer)
1540	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice)+dqcldlscale(1:ngrid,1:nlayer)
1541
1542	! Test energy conservation.
1543	if(enertest)then
1544	lscaledEz(:,:) = cppmass(:,:)dtlscale(:,:)
1545	do ig=1,ngrid
1546	lscaledE(ig) = cppSUM(mass(:,:)dtlscale(:,:))
1547	enddo
1548	call planetwide_sumval(cppmassarea(:,:)dtlscale(:,:)/totarea_planet,dEtot)
1549
1550	if (is_master) print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
1551
1552	! Test water conservation.
1553	call planetwide_sumval(massarea(:,:)dqvaplscale(:,:)ptimestep/totarea_planet+ &
1554	SUM(massarea(:,:)dqcldlscale(:,:))ptimestep/totarea_planet,dWtot)
1555
1556	if (is_master) print*,'In largescale atmospheric water change =',dWtot,' kg m-2'
1557	endif ! end of 'enertest'
1558
1559	! Compute cloud fraction.
1560	do l = 1, nlayer
1561	do ig=1,ngrid
1562	cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
1563	enddo
1564	enddo
1565
1566	endif ! end of 'watercond'
1567
1568	! Water ice / liquid precipitation.
1569	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1570	zdqrain(1:ngrid,1:nlayer,1:nq) = 0.0 !JL18 need to do that everytimestep if mass redis is on.
1571
1572	if(waterrain)then
1573
1574	zdqsrain(1:ngrid) = 0.0
1575	zdqssnow(1:ngrid) = 0.0
1576
1577	call rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,pphi,pt,pdt,pq,pdq, &
1578	zdtrain,zdqrain,zdqsrain,zdqssnow,reevap_precip,cloudfrac)
1579
1580	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) &
1581	+ zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)
1582	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) &
1583	+ zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice)
1584	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain(1:ngrid,1:nlayer)
1585
1586	dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid)
1587	dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
1588
1589	! Test energy conservation.
1590	if(enertest)then
1591
1592	call planetwide_sumval(cppmassarea(:,:)zdtrain(:,:)/totarea_planet,dEtot)
1593	if (is_master) print*,'In rain atmospheric T energy change =',dEtot,' W m-2'
1594	call planetwide_sumval(massarea(:,:)zdqrain(:,:,igcm_h2o_ice)/totarea_planetRLVTT/cpp,dItot_tmp)
1595	call planetwide_sumval(cell_area(:)zdqssnow(:)/totarea_planetRLVTT/cpp,dItot)
1596	dItot = dItot + dItot_tmp
1597	call planetwide_sumval(massarea(:,:)zdqrain(:,:,igcm_h2o_vap)ptimestep/totarea_planet,dVtot_tmp)
1598	call planetwide_sumval(cell_area(:)zdqsrain(:)/totarea_planetRLVTT/cpp,dVtot)
1599	dVtot = dVtot + dVtot_tmp
1600	dEtot = dItot + dVtot
1601
1602	if (is_master) then
1603	print*,'In rain dItot =',dItot,' W m-2'
1604	print*,'In rain dVtot =',dVtot,' W m-2'
1605	print*,'In rain atmospheric L energy change =',dEtot,' W m-2'
1606	endif
1607
1608	! Test water conservation
1609	call planetwide_sumval(massarea(:,:)zdqrain(:,:,igcm_h2o_vap)ptimestep/totarea_planet+ &
1610	massarea(:,:)zdqrain(:,:,igcm_h2o_ice)ptimestep/totarea_planet,dWtot)
1611	call planetwide_sumval((zdqsrain(:)+zdqssnow(:))cell_area(:)ptimestep/totarea_planet,dWtots)
1612
1613	if (is_master) then
1614	print*,'In rain atmospheric water change =',dWtot,' kg m-2'
1615	print*,'In rain surface water change =',dWtots,' kg m-2'
1616	print*,'In rain non-cons factor =',dWtot+dWtots,' kg m-2'
1617	endif
1618
1619	endif ! end of 'enertest'
1620
1621	end if ! enf of 'waterrain'
1622
1623	end if ! end of 'water'
1624
1625	! -------------------------
1626	! VI.2. Photochemistry
1627	! -------------------------
1628
1629	#ifndef MESOSCALE
1630	IF (photochem) then
1631
1632	DO ig=1,ngrid
1633	array_zero1(ig)=0.0
1634	DO l=1,nlayer
1635	array_zero2(ig,l)=0.
1636	ENDDO
1637	ENDDO
1638
1639	call calchim_asis(ngrid,nlayer,nq, &
1640	ptimestep,pplay,pplev,pt,pdt,dist_star,mu0, &
1641	fract,zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim, &
1642	array_zero1,array_zero1, &
1643	pu,pdu,pv,pdv,array_zero2,array_zero2,icount,zdtchim)
1644
1645	! increment values of tracers:
1646	iesp = 0
1647	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1648	! tracers is zero anyways
1649	! September 2020: flag is_chim to increment only on chemical species
1650	IF (is_chim(iq)==1) THEN
1651	iesp = iesp + 1
1652	DO l=1,nlayer
1653	DO ig=1,ngrid
1654	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iesp)
1655	ENDDO
1656	ENDDO
1657	ENDIF
1658	ENDDO ! of DO iq=1,nq
1659
1660
1661	! increment surface values of tracers:
1662	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1663	! tracers is zero anyways
1664	DO ig=1,ngrid
1665	! dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
1666	ENDDO
1667	ENDDO ! of DO iq=1,nq
1668
1669	! increment values of temperature:
1670	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtchim(1:ngrid,1:nlayer)
1671
1672	END IF ! of IF (photochem)
1673	#endif
1674
1675
1676	! ------------------------------------------
1677	! VI.3.a Generic Condensable Species (GCS)
1678	! ------------------------------------------
1679
1680	if (generic_condensation) then
1681	if ((.not.water).and.(.not.calltherm).and.moistadjustment_generic) then
1682	write(,) "moist adjustment for GCS"
1683	dqmoist(1:ngrid,1:nlayer,1:nq)=0.
1684	dtmoist(1:ngrid,1:nlayer)=0.
1685
1686	! Moist Convective Adjustment.
1687	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1688	call moistadj_generic(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
1689
1690	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) &
1691	+ dqmoist(1:ngrid,1:nlayer,1:nq)
1692	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer)
1693
1694	! Test energy conservation.
1695	if(enertest)then
1696	call planetwide_sumval(cppmassarea(:,:)dtmoist(:,:)/totarea_planet,dEtot)
1697	call planetwide_maxval(dtmoist(:,:),dtmoist_max)
1698	call planetwide_minval(dtmoist(:,:),dtmoist_min)
1699	madjdEz(:,:)=cppmass(:,:)dtmoist(:,:)
1700
1701	do ig=1,ngrid
1702	madjdE(ig) = cppSUM(mass(:,:)dtmoist(:,:))
1703	enddo
1704
1705	if (is_master) then
1706	print*,'In moistadj_generic atmospheric energy change =',dEtot,' W m-2'
1707	print,'In moistadj_generic MAX atmospheric energy change =',dtmoist_maxptimestep,'K/step'
1708	print,'In moistadj_generic MIN atmospheric energy change =',dtmoist_minptimestep,'K/step'
1709	endif
1710
1711	! igcm_generic_vap & igcm_generic_ice are not declared below
1712	call planetwide_sumval(massarea(:,:)dqmoist(:,:,1)ptimestep/totarea_planet+ &
1713	massarea(:,:)dqmoist(:,:,2)ptimestep/totarea_planet,dWtot)
1714	if (is_master) print*,'In moistadj_generic atmospheric GCS change =',dWtot,' kg m-2'
1715
1716	endif ! end of 'enertest'
1717	else
1718	! rneb_man, dqmoist are output of moistadj, but used later on
1719	! so set them to zero if moistadj is not called
1720	rneb_man(:,:)=0
1721	dqmoist(:,:,:)=0
1722	endif ! end of '(.not.calltherm).and.moistadjustment_generic'
1723
1724	call condensation_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay, &
1725	pt,pq,pdt,pdq,dt_generic_condensation, &
1726	dqvaplscale_generic,dqcldlscale_generic,rneb_generic)
1727	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dt_generic_condensation(1:ngrid,1:nlayer)
1728	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+dqvaplscale_generic(1:ngrid,1:nlayer,1:nq)
1729	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+dqcldlscale_generic(1:ngrid,1:nlayer,1:nq)
1730
1731	if(enertest)then
1732	do ig=1,ngrid
1733	genericconddE(ig) = cppSUM(mass(:,:)dt_generic_condensation(:,:))
1734	enddo
1735
1736	call planetwide_sumval(cppmassarea(:,:)dt_generic_condensation(:,:)/totarea_planet,dEtot)
1737
1738	if (is_master) print*,'In generic condensation atmospheric energy change =',dEtot,' W m-2'
1739	end if
1740
1741	if (.not. water) then
1742	! Compute GCS (Generic Condensable Specie) cloud fraction. For now we can not have both water cloud fraction and GCS cloud fraction
1743	! Water is the priority
1744	! If you have set water and generic_condensation, then cloudfrac will be water cloudfrac
1745	!
1746	! If you have set generic_condensation (and not water) and you have set several GCS
1747	! then cloudfrac will be only the cloudfrac of the last generic tracer
1748	! (Because it is rewritten every tracer in the loop)
1749	!
1750	! Maybe one should create a cloudfrac_generic(ngrid,nlayer,nq) with 3 dimensions, the last one for tracers
1751
1752	! Let's loop on tracers
1753	cloudfrac(:,:)=0.0
1754	do iq=1,nq
1755	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
1756	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
1757	do l = 1, nlayer
1758	do ig=1,ngrid
1759	cloudfrac(ig,l)=rneb_generic(ig,l,iq)
1760	enddo
1761	enddo
1762	endif
1763	end do ! do iq=1,nq loop on tracers
1764	endif ! .not. water
1765
1766	endif !generic_condensation
1767
1768	!Generic Rain
1769
1770	if (generic_rain) then
1771
1772	zdqsrain_generic(1:ngrid,1:nq) = 0.0
1773	zdqssnow_generic(1:ngrid,1:nq) = 0.0
1774
1775	call rain_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay,pphi,pt,pdt,pq,pdq, &
1776	zdtrain_generic,dq_rain_generic_vap,dq_rain_generic_cld, &
1777	zdqsrain_generic,zdqssnow_generic,reevap_precip_generic,cloudfrac)
1778
1779	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) &
1780	+ dq_rain_generic_vap(1:ngrid,1:nlayer,1:nq)
1781	! only the parts with indexes of generic vapor tracers are filled in dq_rain_generic_vap, other parts are 0.
1782	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) &
1783	+ dq_rain_generic_cld(1:ngrid,1:nlayer,1:nq)
1784	! only the parts with indexes of generic ice(cloud) tracers are filled in dq_rain_generic_cld, other parts are 0.
1785
1786	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain_generic(1:ngrid,1:nlayer)
1787
1788	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq)+zdqsrain_generic(1:ngrid,1:nq)
1789
1790	! Test energy conservation.
1791	if(enertest)then
1792
1793	call planetwide_sumval(cppmassarea(:,:)zdtrain_generic(:,:)/totarea_planet,dEtot)
1794	if (is_master) print*,'In rain_generic atmospheric T energy change =',dEtot,' W m-2'
1795
1796	! Test conservationfor each generic condensable specie (GCS) tracer
1797
1798	! Let's loop on tracers
1799
1800	do iq=1,nq
1801
1802	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
1803
1804	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
1805
1806	call planetwide_sumval(massarea(:,:)dq_rain_generic_cld(:,:,igcm_generic_ice)/totarea_planetRLVTT/cpp,dItot_tmp)
1807	call planetwide_sumval(cell_area(:)zdqssnow_generic(:,igcm_generic_ice)/totarea_planetRLVTT/cpp,dItot)
1808	dItot = dItot + dItot_tmp
1809	call planetwide_sumval(massarea(:,:)dq_rain_generic_vap(:,:,igcm_generic_vap)/totarea_planetRLVTT/cpp,dVtot_tmp)
1810	call planetwide_sumval(cell_area(:)zdqsrain_generic(:,igcm_generic_ice)/totarea_planetRLVTT/cpp,dVtot)
1811	dVtot = dVtot + dVtot_tmp
1812	dEtot = dItot + dVtot
1813
1814	if (is_master) then
1815	print*,'In rain_generic dItot =',dItot,' W m-2'
1816	print*,'In rain_generic dVtot =',dVtot,' W m-2'
1817	print*,'In rain_generic atmospheric L energy change =',dEtot,' W m-2'
1818	endif
1819
1820	call planetwide_sumval(massarea(:,:)dq_rain_generic_vap(:,:,igcm_generic_vap)ptimestep/totarea_planet+ &
1821	massarea(:,:)dq_rain_generic_cld(:,:,igcm_generic_ice)ptimestep/totarea_planet,dWtot)
1822	call planetwide_sumval((zdqsrain_generic(:,igcm_generic_ice)+zdqssnow_generic(:,igcm_generic_ice))cell_area(:)ptimestep/totarea_planet,dWtots)
1823
1824	if (is_master) then
1825	print*,'In rain_generic atmospheric generic tracer change =',dWtot,' kg m-2'
1826	print*,'In rain_generic surface generic tracer change =',dWtots,' kg m-2'
1827	print*,'In rain_generic non-cons factor =',dWtot+dWtots,' kg m-2'
1828	endif
1829
1830	endif
1831
1832	end do ! do iq=1,nq loop on tracers
1833
1834	endif ! end of 'enertest'
1835
1836	endif !generic_rain
1837
1838	! -------------------------
1839	! VI.3.b. Sedimentation.
1840	! -------------------------
1841	if (sedimentation) then
1842
1843	zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
1844	zdqssed(1:ngrid,1:nq) = 0.0
1845
1846	if(watertest)then
1847
1848	iq=igcm_h2o_ice
1849	call planetwide_sumval(massarea(:,:)pq(:,:,iq)ptimestep/totarea_planet,dWtot)
1850	call planetwide_sumval(massarea(:,:)pdq(:,:,iq)ptimestep/totarea_planet,dWtots)
1851	if (is_master) then
1852	print*,'Before sedim pq =',dWtot,' kg m-2'
1853	print*,'Before sedim pdq =',dWtots,' kg m-2'
1854	endif
1855	endif
1856
1857	call callsedim(ngrid,nlayer,ptimestep, &
1858	pplev,zzlev,pt,pdt,pq,pdq, &
1859	zdqsed,zdqssed,nq)
1860
1861	if(watertest)then
1862	iq=igcm_h2o_ice
1863	call planetwide_sumval(massarea(:,:)pq(:,:,iq)ptimestep/totarea_planet,dWtot)
1864	call planetwide_sumval(massarea(:,:)pdq(:,:,iq)ptimestep/totarea_planet,dWtots)
1865	if (is_master) then
1866	print*,'After sedim pq =',dWtot,' kg m-2'
1867	print*,'After sedim pdq =',dWtots,' kg m-2'
1868	endif
1869	endif
1870
1871	! Whether it falls as rain or snow depends only on the surface temperature
1872	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
1873	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
1874
1875	! Test water conservation
1876	if(watertest)then
1877	call planetwide_sumval(massarea(:,:)(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice))ptimestep/totarea_planet,dWtot)
1878	call planetwide_sumval((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))cell_area(:)ptimestep/totarea_planet,dWtots)
1879	if (is_master) then
1880	print*,'In sedim atmospheric ice change =',dWtot,' kg m-2'
1881	print*,'In sedim surface ice change =',dWtots,' kg m-2'
1882	print*,'In sedim non-cons factor =',dWtot+dWtots,' kg m-2'
1883	endif
1884	endif
1885
1886	end if ! end of 'sedimentation'
1887
1888
1889	! ---------------
1890	! VI.4. Updates
1891	! ---------------
1892
1893	! Updating Atmospheric Mass and Tracers budgets.
1894	if(mass_redistrib) then
1895
1896	zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * &
1897	( zdqevap(1:ngrid,1:nlayer) &
1898	+ zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap) &
1899	+ dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap) &
1900	+ dqvaplscale(1:ngrid,1:nlayer) )
1901
1902	do ig = 1, ngrid
1903	zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
1904	enddo
1905
1906	call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
1907	call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
1908	call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
1909
1910	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
1911	rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
1912	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
1913	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
1914
1915	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
1916	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
1917	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
1918	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
1919	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
1920	pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
1921	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
1922
1923	endif
1924
1925
1926	! ------------------
1927	! VI.5. Slab Ocean
1928	! ------------------
1929
1930	if (ok_slab_ocean)then
1931
1932	do ig=1,ngrid
1933	qsurfint(:,igcm_h2o_ice)=qsurf(:,igcm_h2o_ice)
1934	enddo
1935
1936	!! These two commented functions are from BC2014 (only kept for record)
1937	!!
1938	!! call ocean_slab_ice(ptimestep, &
1939	!! ngrid, knindex, tsea_ice, fluxrad, &
1940	!! zdqssnow, qsurf(:,igcm_h2o_ice), &
1941	!! - zdqsdif(:,igcm_h2o_vap), &
1942	!! flux_sens_lat,tsea_ice, pctsrf_sic, &
1943	!! taux,tauy,icount)
1944	!! call ocean_slab_get_vars(ngrid,tslab, &
1945	!! sea_ice, flux_o, &
1946	!! flux_g, dt_hdiff, &
1947	!! dt_ekman)
1948
1949	call ocean_slab_noice(icount, ptimestep, knon, knindex, &
1950	zdqssnow, tsea_ice, &
1951	fluxrad, qsurf(:,igcm_h2o_ice), flux_sens_lat, &
1952	tsea_ice, taux, tauy, zmasq)
1953
1954	call ocean_slab_ice(icount, ptimestep, knon, knindex, &
1955	zdqssnow, tsea_ice, &
1956	fluxrad, qsurf(:,igcm_h2o_ice), flux_sens_lat, &
1957	tsea_ice, -zdqsdif(:,igcm_h2o_vap), taux, tauy, zmasq)
1958
1959	call ocean_slab_frac(pctsrf_sic, zmasq)
1960
1961	call ocean_slab_get_vars(ngrid, tslab, sea_ice, flux_g, &
1962	dt_hdiff, dt_ekman)
1963
1964	!! sea_ice defines the sea ice thickness in kg/m2
1965	!! pctsrf_sic defines the percentage of the oceanic grid that is covered by sea ice (between 0 and almost 1)
1966	!! qsurf(:,igcm_h2o_ice) defines the amount of snow that accumulates on top of the sea ice (in kg/m2) ; note that the snow fraction depends on the amount of snow.
1967	!! therefore, beware that qsurf(:,igcm_h2o_ice) has its own meaning for oceanic grid points (when ok_slab_ocean activated)
1968
1969
1970	do ig=1,ngrid
1971	if (nint(rnat(ig)).eq.1)then
1972	tslab(ig,1) = 0.
1973	tslab(ig,2) = 0.
1974	tsea_ice(ig) = 0.
1975	! tice(ig) = 0.
1976	sea_ice(ig) = 0.
1977	pctsrf_sic(ig) = 0.
1978	qsurf(ig,igcm_h2o_ice) = qsurfint(ig,igcm_h2o_ice)
1979	endif
1980	enddo
1981
1982	endif ! end of 'ok_slab_ocean'
1983
1984
1985	! -----------------------------
1986	! VI.6. Surface Tracer Update
1987	! -----------------------------
1988
1989	if(hydrology)then
1990
1991	call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, &
1992	capcal,albedo,albedo_bareground, &
1993	albedo_snow_SPECTV,albedo_co2_ice_SPECTV, &
1994	mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic, &
1995	sea_ice)
1996
1997	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurf_hyd(1:ngrid)
1998	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + dqs_hyd(1:ngrid,1:nq)
1999
2000
2001	qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
2002
2003	! Test energy conservation
2004	if(enertest)then
2005	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf_hyd(:)/totarea_planet,dEtots)
2006	if (is_master) print*,'In hydrol surface energy change =',dEtots,' W m-2'
2007	endif
2008
2009	! test water conservation
2010	if(watertest)then
2011	call planetwide_sumval(dqs_hyd(:,igcm_h2o_ice)cell_area(:)ptimestep/totarea_planet,dWtots)
2012	if (is_master) print*,'In hydrol surface ice change =',dWtots,' kg m-2'
2013	call planetwide_sumval(dqs_hyd(:,igcm_h2o_vap)cell_area(:)ptimestep/totarea_planet,dWtots)
2014	if (is_master) then
2015	print*,'In hydrol surface water change =',dWtots,' kg m-2'
2016	print*,'---------------------------------------------------------------'
2017	endif
2018	endif
2019
2020	else ! of if (hydrology)
2021
2022	qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
2023
2024	end if ! of if (hydrology)
2025
2026	! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
2027	! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
2028	qsurf_hist(:,:) = qsurf(:,:)
2029
2030	endif! end of if 'tracer'
2031
2032
2033	!------------------------------------------------
2034	! VII. Surface and sub-surface soil temperature
2035	!------------------------------------------------
2036
2037
2038	! Increment surface temperature
2039	if(ok_slab_ocean)then
2040	do ig=1,ngrid
2041	if (nint(rnat(ig)).eq.0)then
2042	zdtsurf(ig)= (tslab(ig,1) &
2043	+ pctsrf_sic(ig)*(tsea_ice(ig)-tslab(ig,1))-tsurf(ig))/ptimestep
2044	endif
2045	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
2046	enddo
2047
2048	else
2049	tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
2050	endif ! end of 'ok_slab_ocean'
2051
2052
2053	! Compute soil temperatures and subsurface heat flux.
2054	if (callsoil) then
2055	call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, &
2056	ptimestep,tsurf,tsoil,capcal,fluxgrd)
2057	endif
2058
2059
2060	if (ok_slab_ocean) then
2061
2062	do ig=1,ngrid
2063
2064	fluxgrdocean(ig)=fluxgrd(ig)
2065	if (nint(rnat(ig)).eq.0) then
2066	capcal(ig)=capcalocean
2067	fluxgrd(ig)=0.
2068	!! fluxgrdocean(ig)=pctsrf_sic(ig)flux_g(ig)+(1-pctsrf_sic(ig))(dt_hdiff(ig,1)+dt_ekman(ig,1))
2069	! Dividing by cell area to have flux in W/m2
2070	fluxgrdocean(ig)=flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1))/cell_area(ig)
2071	do iq=1,nsoilmx
2072	tsoil(ig,iq)=tsurf(ig)
2073	enddo
2074	if (pctsrf_sic(ig).gt.0.5) then
2075	capcal(ig)=capcalseaice
2076	if (qsurf(ig,igcm_h2o_ice).gt.0.) then
2077	capcal(ig)=capcalsno
2078	endif
2079	endif
2080	endif
2081
2082	enddo
2083
2084	endif !end of 'ok_slab_ocean'
2085
2086
2087	! Test energy conservation
2088	if(enertest)then
2089	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf(:)/totarea_planet,dEtots)
2090	if (is_master) print*,'Surface energy change =',dEtots,' W m-2'
2091	endif
2092
2093
2094	!---------------------------------------------------
2095	! VIII. Perform diagnostics and write output files
2096	!---------------------------------------------------
2097
2098	! Note : For output only: the actual model integration is performed in the dynamics.
2099
2100
2101
2102	! Temperature, zonal and meridional winds.
2103	zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
2104	zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
2105	zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
2106
2107	! Recast thermal plume vertical velocity array for outputs
2108	IF (calltherm) THEN
2109	DO ig=1,ngrid
2110	DO l=1,nlayer
2111	zw2_bis(ig,l) = zw2(ig,l)
2112	fm_bis(ig,l) = fm(ig,l)
2113	ENDDO
2114	ENDDO
2115	ENDIF
2116
2117	! Diagnostic.
2118	zdtdyn(1:ngrid,1:nlayer) = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep
2119	ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
2120
2121	zdudyn(1:ngrid,1:nlayer) = (pu(1:ngrid,1:nlayer)-zuprevious(1:ngrid,1:nlayer)) / ptimestep
2122	zuprevious(1:ngrid,1:nlayer) = zu(1:ngrid,1:nlayer)
2123
2124	if(firstcall)then
2125	zdtdyn(1:ngrid,1:nlayer)=0.0
2126	zdudyn(1:ngrid,1:nlayer)=0.0
2127	endif
2128
2129	! Dynamical heating diagnostic.
2130	do ig=1,ngrid
2131	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp
2132	enddo
2133
2134	! Tracers.
2135	zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
2136
2137	! Surface pressure.
2138	ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
2139
2140
2141
2142	! Surface and soil temperature information
2143	call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
2144	call planetwide_minval(tsurf(:),Ts2)
2145	call planetwide_maxval(tsurf(:),Ts3)
2146	if(callsoil)then
2147	TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
2148	if (is_master) then
2149	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
2150	print*,Ts1,Ts2,Ts3,TsS
2151	end if
2152	else
2153	if (is_master) then
2154	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
2155	print*,Ts1,Ts2,Ts3
2156	endif
2157	end if
2158
2159
2160	! Check the energy balance of the simulation during the run
2161	if(corrk)then
2162
2163	call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
2164	call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
2165	call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
2166	call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
2167	call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
2168	do ig=1,ngrid
2169	if(fluxtop_dn(ig).lt.0.0)then
2170	print*,'fluxtop_dn has gone crazy'
2171	print*,'fluxtop_dn=',fluxtop_dn(ig)
2172	print*,'tau_col=',tau_col(ig)
2173	print*,'aerosol=',aerosol(ig,:,:)
2174	print*,'temp= ',pt(ig,:)
2175	print*,'pplay= ',pplay(ig,:)
2176	call abort
2177	endif
2178	end do
2179
2180	if(ngrid.eq.1)then
2181	DYN=0.0
2182	endif
2183
2184	if (is_master) then
2185	print*,' ISR ASR OLR GND DYN [W m^-2]'
2186	print*, ISR,ASR,OLR,GND,DYN
2187	endif
2188
2189	if(enertest .and. is_master)then
2190	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
2191	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
2192	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
2193	endif
2194
2195	if(meanOLR .and. is_master)then
2196	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
2197	! to record global radiative balance
2198	open(92,file="rad_bal.out",form='formatted',position='append')
2199	write(92,*) zday,ISR,ASR,OLR
2200	close(92)
2201	open(93,file="tem_bal.out",form='formatted',position='append')
2202	if(callsoil)then
2203	write(93,*) zday,Ts1,Ts2,Ts3,TsS
2204	else
2205	write(93,*) zday,Ts1,Ts2,Ts3
2206	endif
2207	close(93)
2208	endif
2209	endif
2210
2211	endif ! end of 'corrk'
2212
2213
2214	! Diagnostic to test radiative-convective timescales in code.
2215	if(testradtimes)then
2216	open(38,file="tau_phys.out",form='formatted',position='append')
2217	ig=1
2218	do l=1,nlayer
2219	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
2220	enddo
2221	close(38)
2222	print*,'As testradtimes enabled,'
2223	print*,'exiting physics on first call'
2224	call abort
2225	endif
2226
2227
2228	! Compute column amounts (kg m-2) if tracers are enabled.
2229	if(tracer)then
2230	qcol(1:ngrid,1:nq)=0.0
2231	do iq=1,nq
2232	do ig=1,ngrid
2233	qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
2234	enddo
2235	enddo
2236
2237	! Generalised for arbitrary aerosols now. By LK
2238	reffcol(1:ngrid,1:naerkind)=0.0
2239	if(co2cond.and.(iaero_co2.ne.0))then
2240	call co2_reffrad(ngrid,nlayer,nq,zq,reffrad(1,1,iaero_co2))
2241	do ig=1,ngrid
2242	reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayer,igcm_co2_ice)reffrad(ig,1:nlayer,iaero_co2)mass(ig,1:nlayer))
2243	enddo
2244	endif
2245	if(water.and.(iaero_h2o.ne.0))then
2246	call h2o_reffrad(ngrid,nlayer,zq(1,1,igcm_h2o_ice),zt, &
2247	reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
2248	do ig=1,ngrid
2249	reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayer,igcm_h2o_ice)reffrad(ig,1:nlayer,iaero_h2o)mass(ig,1:nlayer))
2250	enddo
2251	endif
2252
2253	endif ! end of 'tracer'
2254
2255
2256	! Test for water conservation.
2257	if(water)then
2258
2259	call planetwide_sumval(cell_area(:)*qsurf_hist(:,igcm_h2o_ice)/totarea_planet,icesrf)
2260	call planetwide_sumval(cell_area(:)*qsurf_hist(:,igcm_h2o_vap)/totarea_planet,liqsrf)
2261	call planetwide_sumval(cell_area(:)*qcol(:,igcm_h2o_ice)/totarea_planet,icecol)
2262	call planetwide_sumval(cell_area(:)*qcol(:,igcm_h2o_vap)/totarea_planet,vapcol)
2263
2264	h2otot = icesrf + liqsrf + icecol + vapcol
2265
2266	if (is_master) then
2267	print*,' Total water amount [kg m^-2]: ',h2otot
2268	print*,' Surface ice Surface liq. Atmos. con. Atmos. vap. [kg m^-2] '
2269	print*, icesrf,liqsrf,icecol,vapcol
2270	endif
2271
2272	if(meanOLR .and. is_master)then
2273	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
2274	! to record global water balance
2275	open(98,file="h2o_bal.out",form='formatted',position='append')
2276	write(98,*) zday,icesrf,liqsrf,icecol,vapcol
2277	close(98)
2278	endif
2279	endif
2280
2281	endif
2282
2283
2284	! Calculate RH (Relative Humidity) for diagnostic.
2285	if(water)then
2286
2287	do l = 1, nlayer
2288	do ig=1,ngrid
2289	call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l))
2290	RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
2291	enddo
2292	enddo
2293
2294	! Compute maximum possible H2O column amount (100% saturation).
2295	do ig=1,ngrid
2296	H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
2297	enddo
2298
2299	endif ! end of 'water'
2300
2301	! Calculate RH_generic (Generic Relative Humidity) for diagnostic.
2302	if(generic_condensation)then
2303	RH_generic(:,:,:)=0.0
2304	do iq=1,nq
2305
2306	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
2307
2308	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
2309
2310	do l = 1, nlayer
2311	do ig=1,ngrid
2312	call Psat_generic(zt(ig,l),pplay(ig,l),metallicity,psat_tmp_generic,qsat_generic(ig,l,iq))
2313	RH_generic(ig,l,iq) = zq(ig,l,igcm_generic_vap) / qsat_generic(ig,l,iq)
2314	enddo
2315	enddo
2316
2317	end if
2318
2319	end do ! iq=1,nq
2320
2321	endif ! end of 'generic_condensation'
2322
2323
2324	if (is_master) print,'--> Ls =',zls180./pi
2325
2326
2327	! -------------------------
2328	! IX. Thermosphere
2329	! -------------------------
2330	if (callthermos) then
2331	call conduction(ngrid,nlayer,nq,ptimestep,pplay,pplev, &
2332	pt,tsurf,zzlev,zzlay,muvar,pq,firstcall,zdtconduc)
2333	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtconduc(1:ngrid,1:nlayer)
2334	endif ! of if (callthermos)
2335
2336
2337	!----------------------------------------------------------------------
2338	! Writing NetCDF file "RESTARTFI" at the end of the run
2339	!----------------------------------------------------------------------
2340
2341	! Note: 'restartfi' is stored just before dynamics are stored
2342	! in 'restart'. Between now and the writting of 'restart',
2343	! there will have been the itau=itau+1 instruction and
2344	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
2345	! thus we store for time=time+dtvr
2346
2347
2348
2349	if(lastcall) then
2350	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
2351
2352	#ifndef MESOSCALE
2353
2354	if (ngrid.ne.1) then
2355	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
2356
2357	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
2358	ptimestep,ztime_fin, &
2359	tsurf,tsoil,emis,q2,qsurf_hist, &
2360	cloudfrac,totcloudfrac,hice, &
2361	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
2362	endif
2363	#endif
2364	if(ok_slab_ocean) then
2365	call ocean_slab_final!(tslab, seaice)
2366	end if
2367
2368	endif ! end of 'lastcall'
2369
2370
2371	! -----------------------------------------------------------------
2372	! WSTATS: Saving statistics
2373	! -----------------------------------------------------------------
2374	! ("stats" stores and accumulates key variables in file "stats.nc"
2375	! which can later be used to make the statistic files of the run:
2376	! if flag "callstats" from callphys.def is .true.)
2377
2378
2379	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
2380	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2381	call wstats(ngrid,"fluxsurf_lw", &
2382	"Thermal IR radiative flux to surface","W.m-2",2, &
2383	fluxsurf_lw)
2384	call wstats(ngrid,"fluxtop_lw", &
2385	"Thermal IR radiative flux to space","W.m-2",2, &
2386	fluxtop_lw)
2387
2388	! call wstats(ngrid,"fluxsurf_sw", &
2389	! "Solar radiative flux to surface","W.m-2",2, &
2390	! fluxsurf_sw_tot)
2391	! call wstats(ngrid,"fluxtop_sw", &
2392	! "Solar radiative flux to space","W.m-2",2, &
2393	! fluxtop_sw_tot)
2394
2395
2396	call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
2397	call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
2398	call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
2399	!call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
2400	!call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
2401	call wstats(ngrid,"p","Pressure","Pa",3,pplay)
2402	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
2403	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
2404	call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
2405	call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
2406	call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
2407
2408	if (tracer) then
2409	do iq=1,nq
2410	call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
2411	call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2412	'kg m^-2',2,qsurf(1,iq) )
2413	call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2414	'kg m^-2',2,qcol(1,iq) )
2415
2416	! call wstats(ngrid,trim(noms(iq))//'_reff', &
2417	! trim(noms(iq))//'_reff', &
2418	! 'm',3,reffrad(1,1,iq))
2419
2420	end do
2421
2422	if (water) then
2423	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
2424	call wstats(ngrid,"vmr_h2ovapor", &
2425	"H2O vapour volume mixing ratio","mol/mol", &
2426	3,vmr)
2427	endif
2428
2429	endif ! end of 'tracer'
2430
2431	if(watercond.and.CLFvarying)then
2432	call wstats(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
2433	call wstats(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
2434	call wstats(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
2435	call wstats(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
2436	call wstats(ngrid,"RH","relative humidity"," ",3,RH)
2437	endif
2438
2439	if (ok_slab_ocean) then
2440	call wstats(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
2441	call wstats(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
2442	call wstats(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1))
2443	call wstats(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2))
2444	call wstats(ngrid,"tslab1","tslab1","K",2,tslab(:,1))
2445	call wstats(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
2446	call wstats(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
2447	call wstats(ngrid,"tsea_ice","top layer temp, snow/ice","K",2,tsea_ice)
2448	! call wstats(ngrid,"tice","ice temp if snow on top","K",2,tice)
2449	call wstats(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
2450	call wstats(ngrid,"rnat","nature of the surface","",2,rnat)
2451	endif
2452
2453	if(lastcall.and.callstats) then
2454	write (,) "Writing stats..."
2455	call mkstats(ierr)
2456	endif
2457
2458
2459	#ifndef MESOSCALE
2460
2461	!-----------------------------------------------------------------------------------------------------
2462	! OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
2463	!
2464	! Note 1 : output with period "ecritphy", set in "run.def"
2465	!
2466	! Note 2 : writediagfi can also be called from any other subroutine for any variable,
2467	! but its preferable to keep all the calls in one place ...
2468	!-----------------------------------------------------------------------------------------------------
2469
2470	call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
2471	call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
2472	call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
2473	call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
2474	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2475	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
2476	call writediagfi(ngrid,"temp","temperature","K",3,zt)
2477	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
2478	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
2479	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
2480	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
2481	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
2482
2483	! Subsurface temperatures
2484	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2485	! call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
2486
2487	! Oceanic layers
2488	if(ok_slab_ocean) then
2489	call writediagfi(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
2490	call writediagfi(ngrid,"tsea_ice","top layer temp, snow/ice","K",2,tsea_ice)
2491	! call writediagfi(ngrid,"tice","ice temp if snow on top","K",2,tice)
2492	call writediagfi(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
2493	call writediagfi(ngrid,"tslab1","tslab1","K",2,tslab(:,1))
2494	call writediagfi(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
2495	call writediagfi(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
2496	call writediagfi(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
2497	call writediagfi(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1))
2498	call writediagfi(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2))
2499	call writediagfi(ngrid,"rnat","nature of the surface","",2,rnat)
2500	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
2501	call writediagfi(ngrid,"latentFlux","latent heat flux","w.m^-2",2,zdqsdif(:,igcm_h2o_vap)*RLVTT)
2502	endif
2503
2504	! Thermal plume model
2505	if (calltherm) then
2506	call writediagfi(ngrid,'entr','Entrainment','kg m$^{-2}$ s$^{-1}$', 3, entr)
2507	call writediagfi(ngrid,'detr','Detrainment','kg m$^{-2}$ s$^{-1}$', 3, detr)
2508	call writediagfi(ngrid,'fm','Mass flux','kg m$^{-2}$ s$^{-1}$', 3, fm_bis)
2509	call writediagfi(ngrid,'w_plm','Squared vertical velocity','m s$^{-1}$', 3, zw2_bis)
2510	call writediagfi(ngrid,'fraca','Updraft fraction','', 3, fraca)
2511	endif
2512
2513	! GW non-oro outputs
2514	if (calllott_nonoro) then
2515	call WRITEDIAGFI(ngrid,"dugwno","GW non-oro dU","m/s2", 3,d_u_hin)
2516	call WRITEDIAGFI(ngrid,"dvgwno","GW non-oro dV","m/s2", 3,d_v_hin)
2517	endif
2518
2519	! Total energy balance diagnostics
2520	if(callrad.and.(.not.newtonian))then
2521
2522	!call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
2523	!call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
2524	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
2525	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
2526	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
2527	call writediagfi(ngrid,"shad","rings"," ", 2, fract)
2528
2529	! call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
2530	! call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
2531	! call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
2532	! call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
2533	! call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
2534	! call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
2535	! call writediagfi(ngrid,"netfluxsurflw","lw net surface flux","W m-2",2,net_fluxsurf_lw)
2536
2537	if(ok_slab_ocean) then
2538	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrdocean)
2539	else
2540	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
2541	endif
2542
2543	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
2544
2545	endif ! end of 'callrad'
2546
2547	if(enertest) then
2548
2549	if (calldifv) then
2550
2551	call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
2552	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
2553
2554	! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
2555	! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
2556	! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
2557
2558	endif
2559
2560	if (corrk) then
2561	call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
2562	call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
2563	endif
2564
2565	if(watercond) then
2566
2567	call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
2568	if ((.not.calltherm).and.moistadjustment) then
2569	call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
2570	endif
2571	call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
2572
2573	! call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz)
2574	! call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)
2575	! call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
2576
2577	endif
2578
2579	if (generic_condensation) then
2580
2581	call writediagfi(ngrid,"genericconddE","heat from generic condensation","W m-2",2,genericconddE)
2582	call writediagfi(ngrid,"dt_generic_condensation","heating from generic condensation","K s-1",3,dt_generic_condensation)
2583
2584	endif
2585
2586	endif ! end of 'enertest'
2587
2588	! Diagnostics of optical thickness
2589	! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
2590	if (diagdtau) then
2591	do nw=1,L_NSPECTV
2592	write(str2,'(i2.2)') nw
2593	call writediagfi(ngrid,'dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',1,int_dtauv(:,nlayer:1:-1,nw))
2594	enddo
2595	do nw=1,L_NSPECTI
2596	write(str2,'(i2.2)') nw
2597	call writediagfi(ngrid,'dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',1,int_dtaui(:,nlayer:1:-1,nw))
2598	enddo
2599	endif
2600
2601
2602	! Temporary inclusions for heating diagnostics.
2603	call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
2604	call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
2605	call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
2606	call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
2607	call writediagfi(ngrid,"dtmoistadj","moist adj heating","K s-1",3,dtmoist)
2608	call writediagfi(ngrid,"dtdryadj","dry adj heating","K s-1",3,zdtadj)
2609
2610	! For Debugging.
2611	!call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
2612	!call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
2613
2614
2615	! Output aerosols.
2616	if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
2617	call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
2618	if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
2619	call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
2620	if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
2621	call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
2622	if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
2623	call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
2624
2625	! Output tracers.
2626	if (tracer) then
2627
2628	do iq=1,nq
2629	call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
2630	call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2631	'kg m^-2',2,qsurf_hist(1,iq) )
2632	call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', &
2633	'kg m^-2',2,qcol(1,iq) )
2634	! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
2635	! 'kg m^-2',2,qsurf(1,iq) )
2636
2637	if(watercond.or.CLFvarying)then
2638	call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
2639	call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
2640	call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
2641	call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
2642	call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
2643	endif
2644
2645	if(waterrain)then
2646	call writediagfi(ngrid,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
2647	call writediagfi(ngrid,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
2648	call writediagfi(ngrid,"reevap","reevaporation of precipitation","kg m-2 s-1",2,reevap_precip)
2649	endif
2650
2651	if(generic_condensation)then
2652	call writediagfi(ngrid,"rneb_generic","GCS cloud fraction (generic condensation)"," ",3,rneb_generic)
2653	call writediagfi(ngrid,"CLF","GCS cloud fraction"," ",3,cloudfrac)
2654	call writediagfi(ngrid,"RH_generic","GCS relative humidity"," ",3,RH_generic)
2655	endif
2656
2657	if(generic_rain)then
2658	call writediagfi(ngrid,"rain","generic rainfall","kg m-2 s-1",2,zdqsrain_generic)
2659	call writediagfi(ngrid,"snow","generic snowfall","kg m-2 s-1",2,zdqssnow_generic)
2660	call writediagfi(ngrid,"reevap","generic reevaporation of precipitation","kg m-2 s-1",2,reevap_precip_generic)
2661	endif
2662
2663	if((hydrology).and.(.not.ok_slab_ocean))then
2664	call writediagfi(ngrid,"hice","oceanic ice height","m",2,hice)
2665	endif
2666
2667	call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
2668
2669	enddo ! end of 'nq' loop
2670
2671	endif ! end of 'tracer'
2672
2673
2674	! Output spectrum.
2675	if(specOLR.and.corrk)then
2676	call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
2677	call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
2678	call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W/m^2/cm^-1",3,GSR_nu)
2679	endif
2680
2681	#else
2682	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
2683	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
2684	comm_ALBEQ(1:ngrid)=albedo_equivalent(1:ngrid)
2685	if (.not.calldifv) comm_LATENT_HF(:)=0.0
2686	if ((tracer).and.(water)) then
2687	comm_CLOUDFRAC(1:ngrid,1:nlayer)=cloudfrac(1:ngrid,1:nlayer)
2688	comm_TOTCLOUDFRAC(1:ngrid)=totcloudfrac(1:ngrid)
2689	comm_SURFRAIN(1:ngrid)=zdqsrain(1:ngrid)
2690	comm_DQVAP(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_h2o_vap)
2691	comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_h2o_ice)
2692	comm_H2OICE_REFF(1:ngrid,1:nlayer)=reffrad(1:ngrid,1:nlayer,iaero_h2o)
2693	comm_REEVAP(1:ngrid)=reevap_precip(1:ngrid)
2694	comm_DTRAIN(1:ngrid,1:nlayer)=zdtrain(1:ngrid,1:nlayer)
2695	comm_DTLSC(1:ngrid,1:nlayer)=dtlscale(1:ngrid,1:nlayer)
2696	comm_RH(1:ngrid,1:nlayer)=RH(1:ngrid,1:nlayer)
2697
2698	else if ((tracer).and.(generic_condensation).and.(.not. water)) then
2699
2700	! If you have set generic_condensation (and not water) and you have set several GCS
2701	! then the outputs given to WRF will be only the ones for the last generic tracer
2702	! (Because it is rewritten every tracer in the loop)
2703	! WRF can take only one moist tracer
2704
2705	do iq=1,nq
2706	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
2707
2708	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
2709
2710	reffrad_generic_zeros_for_wrf(:,:) = 1.
2711
2712	comm_CLOUDFRAC(1:ngrid,1:nlayer) = cloudfrac(1:ngrid,1:nlayer)
2713	comm_TOTCLOUDFRAC(1:ngrid) = totcloudfrac(1:ngrid) !??????
2714	comm_SURFRAIN(1:ngrid) = zdqsrain_generic(1:ngrid,iq)
2715	comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_vap)
2716	comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_generic_ice)
2717	comm_H2OICE_REFF(1:ngrid,1:nlayer) = reffrad_generic_zeros_for_wrf(1:ngrid,1:nlayer) ! for now zeros !
2718	!comm_H2OICE_REFF(1:ngrid,1:nlayer) = 0*zdtrain_generic(1:ngrid,1:nlayer) ! for now zeros !
2719	comm_REEVAP(1:ngrid) = reevap_precip_generic(1:ngrid,iq)
2720	comm_DTRAIN(1:ngrid,1:nlayer) = zdtrain_generic(1:ngrid,1:nlayer)
2721	comm_DTLSC(1:ngrid,1:nlayer) = dt_generic_condensation(1:ngrid,1:nlayer)
2722	comm_RH(1:ngrid,1:nlayer) = RH_generic(1:ngrid,1:nlayer,iq)
2723
2724	endif
2725	end do ! do iq=1,nq loop on tracers
2726
2727	else
2728	comm_CLOUDFRAC(1:ngrid,1:nlayer)=0.
2729	comm_TOTCLOUDFRAC(1:ngrid)=0.
2730	comm_SURFRAIN(1:ngrid)=0.
2731	comm_DQVAP(1:ngrid,1:nlayer)=0.
2732	comm_DQICE(1:ngrid,1:nlayer)=0.
2733	comm_H2OICE_REFF(1:ngrid,1:nlayer)=0.
2734	comm_REEVAP(1:ngrid)=0.
2735	comm_DTRAIN(1:ngrid,1:nlayer)=0.
2736	comm_DTLSC(1:ngrid,1:nlayer)=0.
2737	comm_RH(1:ngrid,1:nlayer)=0.
2738
2739	endif ! if water, if generic_condensation, else
2740
2741	comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid)
2742	comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid)
2743	comm_FLUXTOP_LW(1:ngrid)=fluxtop_lw(1:ngrid)
2744	comm_FLUXSURF_SW(1:ngrid)=fluxsurf_sw(1:ngrid)
2745	comm_FLUXSURF_LW(1:ngrid)=fluxsurf_lw(1:ngrid)
2746	comm_FLXGRD(1:ngrid)=fluxgrd(1:ngrid)
2747	sensibFlux(1:ngrid) = zflubid(1:ngrid) - capcal(1:ngrid)*zdtsdif(1:ngrid) !!! ????
2748	comm_HR_DYN(1:ngrid,1:nlayer) = zdtdyn(1:ngrid,1:nlayer)
2749
2750	! if (turb_resolved) then
2751	! open(17,file='lsf.txt',form='formatted',status='old')
2752	! rewind(17)
2753	! DO l=1,nlayer
2754	! read(17,*) lsf_dt(l),lsf_dq(l)
2755	! ENDDO
2756	! close(17)
2757	! do ig=1,ngrid
2758	! if ((tracer).and.(water)) then
2759	! pdq(ig,:,igcm_h2o_vap) = pdq(ig,:,igcm_h2o_vap) + lsf_dq(:)
2760	! endif
2761	! pdt(ig,:) = pdt(ig,:) + lsf_dt(:)
2762	! comm_HR_DYN(ig,:) = lsf_dt(:)
2763	! enddo
2764	! endif
2765	#endif
2766
2767	! XIOS outputs
2768	#ifdef CPP_XIOS
2769	! Send fields to XIOS: (NB these fields must also be defined as
2770	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
2771	CALL send_xios_field("ls",zls)
2772
2773	CALL send_xios_field("ps",ps)
2774	CALL send_xios_field("area",cell_area)
2775	CALL send_xios_field("p",pplay)
2776	CALL send_xios_field("temperature",zt)
2777	CALL send_xios_field("u",zu)
2778	CALL send_xios_field("v",zv)
2779	CALL send_xios_field("omega",omega)
2780
2781	IF (calltherm) THEN
2782	CALL send_xios_field('w_plm',zw2_bis)
2783	CALL send_xios_field('entr',entr)
2784	CALL send_xios_field('detr',detr)
2785	! CALL send_xios_field('fm',fm_bis)
2786	! CALL send_xios_field('fraca',fraca)
2787	ENDIF
2788
2789	IF (water) THEN
2790	CALL send_xios_field('h2o_vap',zq(:,:,igcm_h2o_vap))
2791	CALL send_xios_field('h2o_ice',zq(:,:,igcm_h2o_ice))
2792
2793	CALL send_xios_field('h2o_layer1',zq(:,1,igcm_h2o_vap))
2794	CALL send_xios_field('co2_layer1',zq(:,1,igcm_co2_ice))
2795	CALL send_xios_field('tsurf',tsurf)
2796	CALL send_xios_field('co2ice',qsurf(1:ngrid,igcm_co2_ice))
2797	CALL send_xios_field('h2o_ice_s',qsurf(1:ngrid,igcm_h2o_ice))
2798	ENDIF
2799
2800	CALL send_xios_field("ISR",fluxtop_dn)
2801	CALL send_xios_field("OLR",fluxtop_lw)
2802	CALL send_xios_field("ASR",fluxabs_sw)
2803
2804	if (specOLR .and. corrk) then
2805	call send_xios_field("OLR3D",OLR_nu)
2806	call send_xios_field("IR_Bandwidth",DWNI)
2807	call send_xios_field("VI_Bandwidth",DWNV)
2808	call send_xios_field("OSR3D",OSR_nu)
2809	call send_xios_field("GSR3D",GSR_nu)
2810	endif
2811
2812	if (lastcall.and.is_omp_master) then
2813	write(,) "physiq: call xios_context_finalize"
2814	call xios_context_finalize
2815	endif
2816	#endif
2817
2818	if (check_physics_outputs) then
2819	! Check the validity of updated fields at the end of the physics step
2820	call check_physics_fields("end of physiq:", zt, zu, zv, pplev, zq)
2821	endif
2822
2823	icount=icount+1
2824
2825	end subroutine physiq
2826
2827	end module physiq_mod

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

Download in other formats:

Original Format