Context Navigation

physiq_mod.F90 @ 3300

Last change on this file since 3300 was 3299, checked in by emillour, 8 months ago

Generic PCM:

add some auxiliary functions in module tracer_h.F90: is_known_tracer(tname) returns ".true." if "tname" is a known tracer tracer_index(tname) returns the tracer index of tracer "tname"
add the possibility of "volcanic outgasing" as point sources. controled by a "callvolcano" flag (default .false.) details about source location and injection details (frequency, length) need be provided in a "volcano.def" ASCII file read at runtime. (see routine "read_volcano" in volcano.F90)

AB+EM

Property svn:executable set to *

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

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

Download in other formats:

Original Format