Context Navigation

physiq_mod.F90 @ 3026

Last change on this file since 3026 was 2987, checked in by emillour, 18 months ago

Fixed the hardcoded value for the altitude of topmost level

In physiq_mod, the variable zzlev(:,nlayer+1) was hardcoded to 1e7. This is not enough when simulating the atmosphere of inflated hot Jupiter. Thus, to compute this value, we now use the thickness of the level below the top one to compute the altitude of topmost level. This variable is used only in the turbulent diffusion (we use 1/thickness) and the sedimentation (use of thickness).

Property svn:executable set to *

File size: 118.6 KB

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

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

Download in other formats:

Original Format