Context Navigation

physiq_mod.F90 @ 3816

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

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

Download in other formats:

Original Format