Context Navigation

physiq_mod.F90 @ 3441

Last change on this file since 3441 was 3437, checked in by emillour, 9 months ago

Generic PCM Integrate update from newton branch by LT:

remove call to old newtrelax routine and replace it with a new one, newtcool_MOCHA, tuned to the MOCHA intercomparison project

Property svn:executable set to *

File size: 125.3 KB

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

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

Download in other formats:

Original Format