Context Navigation

physiq_mod.F90 @ 3567

Last change on this file since 3567 was 3562, checked in by mmaurice, 4 weeks ago

Generic PCM

1D restart operational: a restart1D.nc file is created that contains
psurf, tracers, winds and temperature profiles. A retartfi.nc file is
also created. Move those to and start1D.nc and startfi.nc and set
"restart" flag to .true. in rcm1d.def to restart from the files (also
make sure that day0 corresponds to the value in startfi.nc).

Property svn:executable set to *

File size: 124.1 KB

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

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

Download in other formats:

Original Format