Context Navigation

physiq_mod.F90 @ 2810

Last change on this file since 2810 was 2803, checked in by aslmd, 3 years ago

Initialisation of Radiative Generic Condensable Aerosols

We can activate the scheme by putting aerogeneric = # of aerosols in callphys.def. This is the only needed thing for activating the radiative effects. They must be tracer in modern tracer.def

Commented out the abort if we use more than 4 aerosols

Added reading of optprop files for Radiative Generic Condensable Aerosols

We use the same file for IR and VI channel. For now, only MnS, Cr,Fe and Mg2SiO4 can be read. If you want to add another specie, check the code, it is explained how to quickly do that (right above the Initializations)

Added radii calculation for Radiative Generic Condensable aerosols

Changed the hardcoded size of the totalemit array

The hardcoded size is now 1900 instead of 100 so we don't exceed the array size when working at high spectral resolution (very rare case)

Added opacity computation for Radiative Generic Condensable aerosols

We do this computation in the same fashion as what's performed on water and dust.

switch iniaerosol and initracer order, to prepare for the RGCS scheme

Needed to switch the order of initialization so we can use the RGCS scheme without the assumption that ice and vap tracer of the same specie are following each other in the traceur.def file

Property svn:executable set to *

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

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

Download in other formats:

Original Format