Context Navigation

physiq_mod.F @ 2932

Last change on this file since 2932 was 2932, checked in by romain.vande, 21 months ago

Mars PCM:
Add a new routine to write output called write_output.
It needs to be imported (for example : use write_output_mod, only: write_output)
Then, it requires only 4 arguments : the name of the variable, its title, its units and the variable itself.
It detects the dimension of the variable and decide to output either in diagfi or diagsoil.
It is also compatible with XIOS (xml file needs to be adapted)
Writediagfi and writediagsoil routines are still available in case.
RV

File size: 156.8 KB

Line
1	MODULE physiq_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	SUBROUTINE physiq(
8	$ ngrid,nlayer,nq
9	$ ,firstcall,lastcall
10	$ ,pday,ptime,ptimestep
11	$ ,pplev,pplay,pphi
12	$ ,pu,pv,pt,pq
13	$ ,flxw
14	$ ,pdu,pdv,pdt,pdq,pdpsrf)
15
16	use watercloud_mod, only: watercloud, zdqcloud, zdqscloud
17	use calchim_mod, only: calchim, ichemistry, zdqchim, zdqschim
18	use watersat_mod, only: watersat
19	use co2condens_mod, only: co2condens
20	use co2cloud_mod, only: co2cloud
21	use callradite_mod, only: callradite
22	use callsedim_mod, only: callsedim
23	use rocketduststorm_mod, only: rocketduststorm, dustliftday
24	use calcstormfract_mod, only: calcstormfract
25	use topmons_mod, only: topmons,topmons_setup
26	use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2, igcm_co2_ice,
27	& igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice,
28	& igcm_hdo_vap, igcm_hdo_ice,
29	& igcm_ccn_mass, igcm_ccn_number,
30	& igcm_ccnco2_mass, igcm_ccnco2_number,
31	& igcm_ccnco2_h2o_mass_ice,
32	& igcm_ccnco2_h2o_mass_ccn,
33	& igcm_ccnco2_h2o_number,
34	& igcm_ccnco2_meteor_mass,
35	& igcm_ccnco2_meteor_number,
36	& rho_ice_co2,nuiceco2_sed,nuiceco2_ref,
37	& igcm_dust_mass, igcm_dust_number, igcm_h2o2,
38	& nuice_ref, rho_ice, rho_dust, ref_r0,
39	& igcm_he, igcm_stormdust_mass,
40	& igcm_stormdust_number, igcm_topdust_mass,
41	& igcm_topdust_number,
42	& qperemin
43	use comsoil_h, only: inertiedat, ! soil thermal inertia
44	& tsoil, nsoilmx,!number of subsurface layers
45	& mlayer,layer ! soil mid layer depths
46	use geometry_mod, only: longitude, latitude, cell_area,
47	& longitude_deg
48	use comgeomfi_h, only: sinlon, coslon, sinlat, coslat
49	use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam,
50	& zthe, z0, albedo_h2o_cap,albedo_h2o_frost,
51	& frost_albedo_threshold,frost_metam_threshold,
52	& tsurf, emis,
53	& capcal, fluxgrd, qsurf,
54	& hmons,summit,base,watercap,watercaptag
55	use comsaison_h, only: dist_sol, declin, zls,
56	& mu0, fract, local_time
57	use slope_mod, only: theta_sl, psi_sl
58	use conc_mod, only: rnew, cpnew, mmean
59	use time_phylmdz_mod, only: iphysiq, day_step, ecritstart, daysec
60	use dimradmars_mod, only: aerosol, totcloudfrac,
61	& dtrad, fluxrad_sky, fluxrad, albedo,
62	& naerkind, iaer_dust_doubleq,
63	& iaer_stormdust_doubleq, iaer_h2o_ice,
64	& flux_1AU
65	use dust_param_mod, only: doubleq, lifting, callddevil,
66	& tauscaling, odpref, dustbin,
67	& dustscaling_mode, dust_rad_adjust,
68	& freedust, reff_driven_IRtoVIS_scenario
69	use turb_mod, only: q2, wstar, ustar, sensibFlux,
70	& zmax_th, hfmax_th, turb_resolved
71	use planete_h, only: aphelie, periheli, year_day, peri_day,
72	& obliquit
73	USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad
74	USE calldrag_noro_mod, ONLY: calldrag_noro
75	USE vdifc_mod, ONLY: vdifc
76	use param_v4_h, only: nreact,n_avog,
77	& fill_data_thermos, allocate_param_thermos
78	use iono_h, only: allocate_param_iono
79	use compute_dtau_mod, only: compute_dtau
80	use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
81	use check_fields_mod, only: check_physics_fields
82	#ifdef MESOSCALE
83	use comsoil_h, only: mlayer,layer
84	use surfdat_h, only: z0_default
85	use comm_wrf
86	#else
87	USE planetwide_mod, ONLY: planetwide_maxval, planetwide_minval,
88	& planetwide_sumval
89	use phyredem, only: physdem0, physdem1
90	use phyetat0_mod, only: phyetat0, tab_cntrl_mod
91	use wstats_mod, only: callstats, wstats, mkstats
92	use eofdump_mod, only: eofdump
93	USE vertical_layers_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt
94	USE mod_phys_lmdz_omp_data, ONLY: is_omp_master
95	USE time_phylmdz_mod, ONLY: day_end
96	#endif
97
98	#ifdef CPP_XIOS
99	use xios_output_mod, only: initialize_xios_output,
100	& update_xios_timestep,
101	& send_xios_field
102	use wxios, only: wxios_context_init, xios_context_finalize
103	#endif
104	USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured
105	use ioipsl_getin_p_mod, only: getin_p
106	use comslope_mod, ONLY: nslope,def_slope,def_slope_mean,
107	& subslope_dist,iflat,sky_slope,
108	& major_slope,compute_meshgridavg,
109	& ini_comslope_h
110	USE ioipsl_getincom, only: getin
111	use write_output_mod, only: write_output
112	IMPLICIT NONE
113	c=======================================================================
114	c
115	c subject:
116	c --------
117	c
118	c Organisation of the physical parametrisations of the LMD
119	c martian atmospheric general circulation model.
120	c
121	c The GCM can be run without or with tracer transport
122	c depending on the value of Logical "tracer" in file "callphys.def"
123	c Tracers may be water vapor, ice OR chemical species OR dust particles
124	c
125	c SEE comments in initracer.F about numbering of tracer species...
126	c
127	c It includes:
128	c
129	c 1. Initialization:
130	c 1.1 First call initializations
131	c 1.2 Initialization for every call to physiq
132	c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
133	c 2. Compute radiative transfer tendencies
134	c (longwave and shortwave) for CO2 and aerosols.
135	c 3. Gravity wave and subgrid scale topography drag :
136	c 4. Vertical diffusion (turbulent mixing):
137	c 5. Convective adjustment
138	c 6. Condensation and sublimation of carbon dioxide.
139	c 7. TRACERS :
140	c 7a. water, water ice, co2 ice (clouds)
141	c 7b. call for photochemistry when tracers are chemical species
142	c 7c. other scheme for tracer (dust) transport (lifting, sedimentation)
143	c 7d. updates (CO2 pressure variations, surface budget)
144	c 8. Contribution to tendencies due to thermosphere
145	c 9. Surface and sub-surface temperature calculations
146	c 10. Write outputs :
147	c - "startfi", "histfi" (if it's time)
148	c - Saving statistics (if "callstats = .true.")
149	c - Dumping eof (if "calleofdump = .true.")
150	c - Output any needed variables in "diagfi"
151	c 11. Diagnostic: mass conservation of tracers
152	c
153	c author:
154	c -------
155	c Frederic Hourdin 15/10/93
156	c Francois Forget 1994
157	c Christophe Hourdin 02/1997
158	c Subroutine completly rewritten by F.Forget (01/2000)
159	c Introduction of the photochemical module: S. Lebonnois (11/2002)
160	c Introduction of the thermosphere module: M. Angelats i Coll (2002)
161	c Water ice clouds: Franck Montmessin (update 06/2003)
162	c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
163	c Nb: See callradite.F for more information.
164	c Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
165	c jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization
166	c
167	c 10/16 J. Audouard: modifications for CO2 clouds scheme
168
169	c arguments:
170	c ----------
171	c
172	c input:
173	c ------
174	c ecri period (in dynamical timestep) to write output
175	c ngrid Size of the horizontal grid.
176	c All internal loops are performed on that grid.
177	c nlayer Number of vertical layers.
178	c nq Number of advected fields
179	c firstcall True at the first call
180	c lastcall True at the last call
181	c pday Number of days counted from the North. Spring
182	c equinoxe.
183	c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
184	c ptimestep timestep (s)
185	c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
186	c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
187	c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
188	c pu(ngrid,nlayer) u component of the wind (ms-1)
189	c pv(ngrid,nlayer) v component of the wind (ms-1)
190	c pt(ngrid,nlayer) Temperature (K)
191	c pq(ngrid,nlayer,nq) Advected fields
192	c pudyn(ngrid,nlayer) \|
193	c pvdyn(ngrid,nlayer) \| Dynamical temporal derivative for the
194	c ptdyn(ngrid,nlayer) \| corresponding variables
195	c pqdyn(ngrid,nlayer,nq) \|
196	c flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
197	c
198	c output:
199	c -------
200	c
201	c pdu(ngrid,nlayer) \|
202	c pdv(ngrid,nlayer) \| Temporal derivative of the corresponding
203	c pdt(ngrid,nlayer) \| variables due to physical processes.
204	c pdq(ngrid,nlayer,nq) \|
205	c pdpsrf(ngrid) \|
206
207	c
208	c=======================================================================
209	c
210	c 0. Declarations :
211	c ------------------
212
213	include "callkeys.h"
214	include "comg1d.h"
215	include "nlteparams.h"
216	include "netcdf.inc"
217
218	c Arguments :
219	c -----------
220
221	c inputs:
222	c -------
223	INTEGER,INTENT(in) :: ngrid ! number of atmospheric columns
224	INTEGER,INTENT(in) :: nlayer ! number of atmospheric layers
225	INTEGER,INTENT(in) :: nq ! number of tracers
226	LOGICAL,INTENT(in) :: firstcall ! signals first call to physics
227	LOGICAL,INTENT(in) :: lastcall ! signals last call to physics
228	REAL,INTENT(in) :: pday ! number of elapsed sols since reference Ls=0
229	REAL,INTENT(in) :: ptime ! "universal time", given as fraction of sol (e.g.: 0.5 for noon)
230	REAL,INTENT(in) :: ptimestep ! physics timestep (s)
231	REAL,INTENT(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
232	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
233	REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential at mid-layer (m2s-2)
234	REAL,INTENT(in) :: pu(ngrid,nlayer) ! zonal wind component (m/s)
235	REAL,INTENT(in) :: pv(ngrid,nlayer) ! meridional wind component (m/s)
236	REAL,INTENT(in) :: pt(ngrid,nlayer) ! temperature (K)
237	REAL,INTENT(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
238	REAL,INTENT(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s)
239	! at lower boundary of layer
240
241	c outputs:
242	c --------
243	c physical tendencies
244	REAL,INTENT(out) :: pdu(ngrid,nlayer) ! zonal wind tendency (m/s/s)
245	REAL,INTENT(out) :: pdv(ngrid,nlayer) ! meridional wind tendency (m/s/s)
246	REAL,INTENT(out) :: pdt(ngrid,nlayer) ! temperature tendency (K/s)
247	REAL,INTENT(out) :: pdq(ngrid,nlayer,nq) ! tracer tendencies (../kg/s)
248	REAL,INTENT(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
249
250
251
252	c Local saved variables:
253	c ----------------------
254	INTEGER,SAVE :: day_ini ! Initial date of the run (sol since Ls=0)
255	INTEGER,SAVE :: icount ! counter of calls to physiq during the run.
256	REAL,SAVE :: time_phys
257
258	!$OMP THREADPRIVATE(day_ini,icount,time_phys)
259
260	#ifdef DUSTSTORM
261	REAL pq_tmp(ngrid, nlayer, 2) ! To compute tendencies due the dust bomb
262	#endif
263
264	c Variables used by the water ice microphysical scheme:
265	REAL rice(ngrid,nlayer) ! Water ice geometric mean radius (m)
266	REAL nuice(ngrid,nlayer) ! Estimated effective variance
267	! of the size distribution
268	real rsedcloud(ngrid,nlayer) ! Cloud sedimentation radius
269	real rhocloud(ngrid,nlayer) ! Cloud density (kg.m-3)
270	real rsedcloudco2(ngrid,nlayer) ! CO2 Cloud sedimentation radius
271	real rhocloudco2(ngrid,nlayer) ! CO2 Cloud density (kg.m-3)
272	real nuiceco2(ngrid,nlayer) ! Estimated effective variance of the
273	! size distribution
274	REAL inertiesoil(ngrid,nsoilmx,nslope) ! Time varying subsurface
275	! thermal inertia (J.s-1/2.m-2.K-1)
276	! (used only when tifeedback=.true.)
277	c Variables used by the CO2 clouds microphysical scheme:
278	DOUBLE PRECISION riceco2(ngrid,nlayer) ! co2 ice geometric mean radius (m)
279	real zdqssed_co2(ngrid) ! CO2 flux at the surface (kg.m-2.s-1)
280	real zdqssed_ccn(ngrid,nq) ! CCN flux at the surface (kg.m-2.s-1)
281	real zcondicea_co2microp(ngrid,nlayer)
282	c Variables used by the photochemistry
283	REAL surfdust(ngrid,nlayer) ! dust surface area (m2/m3, if photochemistry)
284	REAL surfice(ngrid,nlayer) ! ice surface area (m2/m3, if photochemistry)
285	c Variables used by the slope model
286	REAL sl_ls, sl_lct, sl_lat
287	REAL sl_tau, sl_alb, sl_the, sl_psi
288	REAL sl_fl0, sl_flu
289	REAL sl_ra, sl_di0
290	REAL sky
291	REAL fluxsurf_dir_dn_sw(ngrid) ! Incident direct solar flux on Mars at surface (W.m-2)
292
293	REAL,PARAMETER :: stephan = 5.67e-08 ! Stephan Boltzman constant
294
295	c Local variables :
296	c -----------------
297
298	REAL CBRT
299	EXTERNAL CBRT
300
301	! CHARACTER*80 fichier
302	INTEGER l,ig,ierr,igout,iq,tapphys,isoil
303
304	REAL fluxsurf_lw(ngrid,nslope) !incident LW (IR) surface flux (W.m-2)
305	REAL fluxsurf_dn_sw(ngrid,2,nslope) ! Incident SW (solar) surface flux (W.m-2)
306	REAL fluxsurf_up_sw(ngrid,2) ! Reflected SW (solar) surface flux (W.m-2)
307	REAL fluxtop_lw(ngrid) !Outgoing LW (IR) flux to space (W.m-2)
308	REAL fluxtop_dn_sw(ngrid,2) ! Incoming SW (solar) flux from space (W.m-2)
309	REAL fluxtop_up_sw(ngrid,2) ! Outgoing SW (solar) flux to space (W.m-2)
310	REAL tau_pref_scenario(ngrid) ! prescribed dust column visible opacity
311	! at odpref
312	REAL IRtoVIScoef(ngrid) ! conversion coefficient to apply on
313	! scenario absorption IR (9.3um) CDOD
314	! = tau_pref_gcm_VIS / tau_pref_gcm_IR
315	REAL tau_pref_gcm(ngrid) ! dust column visible opacity at odpref in the GCM
316	c rocket dust storm
317	REAL totstormfract(ngrid) ! fraction of the mesh where the dust storm is contained
318	logical clearatm ! clearatm used to calculate twice the radiative
319	! transfer when rdstorm is active :
320	! - in a mesh with stormdust and background dust (false)
321	! - in a mesh with background dust only (true)
322	c entrainment by mountain top dust flows
323	logical nohmons ! nohmons used to calculate twice the radiative
324	! transfer when topflows is active :
325	! - in a mesh with topdust and background dust (false)
326	! - in a mesh with background dust only (true)
327
328	REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
329	! AS: TBD: this one should be in a module !
330	REAL zday ! date (time since Ls=0, in martian days)
331	REAL zzlay(ngrid,nlayer) ! altitude at the middle of the layers
332	REAL zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
333	! REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
334
335	c Tendancies due to various processes:
336	REAL dqsurf(ngrid,nq,nslope) ! tendency for tracers on surface (Kg/m2/s)
337	REAL zdtlw(ngrid,nlayer) ! (K/s)
338	REAL zdtsw(ngrid,nlayer) ! (K/s)
339	REAL pdqrds(ngrid,nlayer,nq) ! tendency for dust after rocketduststorm
340
341	REAL zdtnirco2(ngrid,nlayer) ! (K/s)
342	REAL zdtnlte(ngrid,nlayer) ! (K/s)
343	REAL zdtsurf(ngrid,nslope) ! (K/s)
344	REAL zdtcloud(ngrid,nlayer),zdtcloudco2(ngrid,nlayer)
345	REAL zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! (m.s-2)
346	REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid,nslope) ! (K/s)
347	REAL zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! (m.s-2)
348	REAL zdhadj(ngrid,nlayer) ! (K/s)
349	REAL zdtgw(ngrid,nlayer) ! (K/s)
350	REAL zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer) ! (m.s-2)
351	REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid,nslope)
352	REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer)
353
354	REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq,nslope)
355	REAL zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq)
356	REAL zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq)
357	REAL zdqadj(ngrid,nlayer,nq)
358	REAL zdqc(ngrid,nlayer,nq)
359	REAL zdqcloudco2(ngrid,nlayer,nq)
360	REAL zdqsc(ngrid,nq,nslope)
361
362	REAL zdteuv(ngrid,nlayer) ! (K/s)
363	REAL zdtconduc(ngrid,nlayer) ! (K/s)
364	REAL zdumolvis(ngrid,nlayer)
365	REAL zdvmolvis(ngrid,nlayer)
366	real zdqmoldiff(ngrid,nlayer,nq)
367	real*8 PhiEscH,PhiEscH2,PhiEscD
368
369	REAL dwatercap(ngrid,nslope), dwatercap_dif(ngrid,nslope) ! (kg/m-2)
370
371	c Local variable for local intermediate calcul:
372	REAL zflubid(ngrid,nslope)
373	REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
374	REAL zdum1(ngrid,nlayer)
375	REAL zdum2(ngrid,nlayer)
376	REAL ztim1,ztim2,ztim3, z1,z2
377	REAL ztime_fin
378	REAL zdh(ngrid,nlayer)
379	REAL zh(ngrid,nlayer) ! potential temperature (K)
380	REAL pw(ngrid,nlayer) ! vertical velocity (m/s) (>0 when downwards)
381	INTEGER length
382	PARAMETER (length=100)
383
384	c Variables for the total dust for diagnostics
385	REAL qdusttotal(ngrid,nlayer) !it equals to dust + stormdust
386
387	INTEGER iaer
388
389	c local variables only used for diagnostic (output in file "diagfi" or "stats")
390	c -----------------------------------------------------------------------------
391	REAL ps(ngrid), zt(ngrid,nlayer)
392	REAL zu(ngrid,nlayer),zv(ngrid,nlayer)
393	REAL zq(ngrid,nlayer,nq)
394
395	REAL fluxtop_dn_sw_tot(ngrid), fluxtop_up_sw_tot(ngrid)
396	REAL fluxsurf_dn_sw_tot(ngrid,nslope), fluxsurf_up_sw_tot(ngrid)
397	character*2 str2
398	! character*5 str5
399	real zdtdif(ngrid,nlayer), zdtadj(ngrid,nlayer)
400	real rdust(ngrid,nlayer) ! dust geometric mean radius (m)
401	real rstormdust(ngrid,nlayer) ! stormdust geometric mean radius (m)
402	real rtopdust(ngrid,nlayer) ! topdust geometric mean radius (m)
403	integer igmin, lmin
404	logical tdiag
405
406	real co2col(ngrid) ! CO2 column
407	! pplev and pplay are dynamical inputs and must not be modified in the physics.
408	! instead, use zplay and zplev :
409	REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
410	! REAL zstress(ngrid),cd
411	real tmean, zlocal(nlayer)
412	real rho(ngrid,nlayer) ! density
413	real vmr(ngrid,nlayer) ! volume mixing ratio
414	real rhopart(ngrid,nlayer) ! number density of a given species
415	real colden(ngrid,nq) ! vertical column of tracers
416	real mass(nq) ! global mass of tracers (g)
417	REAL mtot(ngrid) ! Total mass of water vapor (kg/m2)
418	REAL mstormdtot(ngrid) ! Total mass of stormdust tracer (kg/m2)
419	REAL mdusttot(ngrid) ! Total mass of dust tracer (kg/m2)
420	REAL icetot(ngrid) ! Total mass of water ice (kg/m2)
421	REAL mtotco2(ngrid) ! Total mass of co2, including ice at the surface (kg/m2)
422	REAL vaptotco2(ngrid) ! Total mass of co2 vapor (kg/m2)
423	REAL icetotco2(ngrid) ! Total mass of co2 ice (kg/m2)
424	REAL Nccntot(ngrid) ! Total number of ccn (nbr/m2)
425	REAL NccnCO2tot(ngrid) ! Total number of ccnCO2 (nbr/m2)
426	REAL Mccntot(ngrid) ! Total mass of ccn (kg/m2)
427	REAL rave(ngrid) ! Mean water ice effective radius (m)
428	REAL opTES(ngrid,nlayer) ! abs optical depth at 825 cm-1
429	REAL tauTES(ngrid) ! column optical depth at 825 cm-1
430	REAL Qabsice ! Water ice absorption coefficient
431	REAL taucloudtes(ngrid) ! Cloud opacity at infrared
432	! reference wavelength using
433	! Qabs instead of Qext
434	! (direct comparison with TES)
435	REAL mtotD(ngrid) ! Total mass of HDO vapor (kg/m2)
436	REAL icetotD(ngrid) ! Total mass of HDO ice (kg/m2)
437	REAL DoH_vap(ngrid,nlayer) !D/H ratio
438	REAL DoH_ice(ngrid,nlayer) !D/H ratio
439	REAL DoH_surf(ngrid) !D/H ratio surface
440
441	REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s)
442	REAL dndustsurf(ngrid) ! surface n dust flux (number/m2/s)
443	REAL ndust(ngrid,nlayer) ! true n dust (kg/kg)
444	REAL qdust(ngrid,nlayer) ! true q dust (kg/kg)
445	REAL nccn(ngrid,nlayer) ! true n ccn (kg/kg)
446	REAL qccn(ngrid,nlayer) ! true q ccn (kg/kg)
447	c definition tendancies of stormdust tracers
448	REAL rdsdqdustsurf(ngrid) ! surface q stormdust flux (kg/m2/s)
449	REAL rdsdndustsurf(ngrid) ! surface n stormdust flux (number/m2/s)
450	REAL rdsndust(ngrid,nlayer) ! true n stormdust (kg/kg)
451	REAL rdsqdust(ngrid,nlayer) ! true q stormdust (kg/kg)
452	REAL wspeed(ngrid,nlayer+1) ! vertical velocity stormdust tracer
453	REAL wtop(ngrid,nlayer+1) ! vertical velocity topdust tracer
454
455	REAL dsodust(ngrid,nlayer) ! density scaled opacity for background dust
456	REAL dsords(ngrid,nlayer) ! density scaled opacity for stormdust
457	REAL dsotop(ngrid,nlayer) ! density scaled opacity for topdust
458
459	c Test 1d/3d scavenging
460	real h2otot(ngrid)
461	real hdotot(ngrid)
462	REAL satu(ngrid,nlayer) ! satu ratio for output
463	REAL zqsat(ngrid,nlayer) ! saturation
464	REAL satuco2(ngrid,nlayer) ! co2 satu ratio for output
465	REAL zqsatco2(ngrid,nlayer) ! saturation co2
466
467
468	! Added for new NLTE scheme
469
470	real co2vmr_gcm(ngrid,nlayer)
471	real n2vmr_gcm(ngrid,nlayer)
472	real ovmr_gcm(ngrid,nlayer)
473	real covmr_gcm(ngrid,nlayer)
474	integer ierr_nlte
475	real*8 varerr
476
477	C Non-oro GW drag & Calcul of Brunt-Vaisala freq. (BV2)
478	REAL ztetalev(ngrid,nlayer)
479	real zdtetalev(ngrid,nlayer), zdzlev(ngrid,nlayer)
480	REAL bv2(ngrid,nlayer) ! BV2 at zlev
481	c Non-oro GW tendencies
482	REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
483	REAL d_t_hin(ngrid,nlayer)
484	c Diagnostics 2D of gw_nonoro
485	REAL zustrhi(ngrid), zvstrhi(ngrid)
486	c Variables for PBL
487	REAL zz1(ngrid)
488	REAL lmax_th_out(ngrid)
489	REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
490	REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
491	INTEGER lmax_th(ngrid),dimout,n_out,n
492	CHARACTER(50) zstring
493	REAL dtke_th(ngrid,nlayer+1)
494	REAL zcdv(ngrid), zcdh(ngrid)
495	REAL, ALLOCATABLE, DIMENSION(:,:) :: T_out
496	REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out
497	REAL u_out1(ngrid)
498	REAL T_out1(ngrid)
499	REAL, ALLOCATABLE, DIMENSION(:) :: z_out ! height of interpolation between z0 and z1 [meters]
500	REAL tstar(ngrid) ! friction velocity and friction potential temp
501	REAL L_mo(ngrid),vhf(ngrid),vvv(ngrid)
502	real qdustrds0(ngrid,nlayer),qdustrds1(ngrid,nlayer)
503	real qstormrds0(ngrid,nlayer),qstormrds1(ngrid,nlayer)
504	real qdusttotal0(ngrid),qdusttotal1(ngrid)
505
506	c sub-grid scale water ice clouds (A. Pottier 2013)
507	logical clearsky
508	! flux for the part without clouds
509	real zdtswclf(ngrid,nlayer)
510	real zdtlwclf(ngrid,nlayer)
511	real fluxsurf_lwclf(ngrid)
512	real fluxsurf_dn_swclf(ngrid,2),fluxsurf_up_swclf(ngrid,2)
513	real fluxtop_lwclf(ngrid)
514	real fluxtop_dn_swclf(ngrid,2),fluxtop_up_swclf(ngrid,2)
515	real taucloudtesclf(ngrid)
516	real tf_clf, ntf_clf ! tf: fraction of clouds, ntf: fraction without clouds
517	real rave2(ngrid), totrave2(ngrid) ! Mean water ice mean radius (m)
518	C test de conservation de la masse de CO2
519	REAL co2totA
520	REAL co2totB
521	REAL co2conservation
522
523	c entrainment by mountain top dust flows above sub-grid scale topography
524	REAL pdqtop(ngrid,nlayer,nq) ! tendency for dust after topmons
525
526	c when no startfi file is asked for init
527	real alpha,lay1 ! coefficients for building layers
528	integer iloop
529
530	! flags to trigger extra sanity checks
531	logical,save :: check_physics_inputs=.false.
532	logical,save :: check_physics_outputs=.false.
533
534	!$OMP THREADPRIVATE(check_physics_inputs,check_physics_outputs)
535
536	c Sub-grid scale slopes
537	real :: tsurf_meshavg(ngrid) ! Surface temperature grid box averaged [K]
538	real :: albedo_meshavg(ngrid,2) ! albedo temperature grid box averaged [1]
539	real :: emis_meshavg(ngrid,2) ! emis temperature grid box averaged [1]
540	real :: qsurf_meshavg(ngrid,nq) ! surface tracer mesh averaged [kg/m^2]
541	real :: qsurf_tmp(ngrid,nq) ! temporary qsurf for chimie
542	real :: inertiedat_tmp(ngrid,nsoilmx,nslope) ! temporary inertiedat for soil.F
543	integer :: islope
544	real :: zdqsdif_meshavg_tmp(ngrid,nq) ! temporary for dust lifting
545
546	logical :: write_restart
547
548	c=======================================================================
549	pdq(:,:,:) = 0.
550
551	c 1. Initialisation:
552	c -----------------
553	c 1.1 Initialisation only at first call
554	c ---------------------------------------
555
556	IF (firstcall) THEN
557
558	call getin_p("check_physics_inputs",check_physics_inputs)
559	call getin_p("check_physics_outputs",check_physics_outputs)
560
561	c variables set to 0
562	c ~~~~~~~~~~~~~~~~~~
563	aerosol(:,:,:)=0
564	dtrad(:,:)=0
565
566	#ifndef MESOSCALE
567	fluxrad(:,:)=0
568	wstar(:)=0.
569	#endif
570
571	#ifdef CPP_XIOS
572	! Initialize XIOS context
573	write(,) "physiq: call wxios_context_init"
574	CALL wxios_context_init
575	#endif
576
577	c read startfi
578	c ~~~~~~~~~~~~
579	#ifndef MESOSCALE
580
581	! GCM. Read netcdf initial physical parameters.
582	CALL phyetat0 ("startfi.nc",0,0,
583	& nsoilmx,ngrid,nlayer,nq,
584	& day_ini,time_phys,
585	& tsurf,tsoil,albedo,emis,
586	& q2,qsurf,tauscaling,totcloudfrac,wstar,
587	& watercap,def_slope,def_slope_mean,subslope_dist)
588
589	DO islope=1,nslope
590	sky_slope(islope) = (1.+cos(pi*def_slope_mean(islope)/180.))/2.
591	END DO
592
593	! Sky view:
594	DO islope=1,nslope
595	sky_slope(islope) = (1.+cos(pi*def_slope_mean(islope)/180.))/2.
596	END DO
597	! Determine the 'flatest' slopes
598	iflat = 1
599	DO islope=2,nslope
600	IF(abs(def_slope_mean(islope)).lt.
601	& abs(def_slope_mean(iflat)))THEN
602	iflat = islope
603	ENDIF
604	ENDDO
605	PRINT*,'Flat slope for islope = ',iflat
606	PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
607
608	#else
609	! MESOSCALE. Supposedly everything is already set in modules.
610	! So we just check. And we fill day_ini
611	print*,"check: --- in physiq.F"
612	print*,"check: rad,cpp,g,r,rcp,daysec"
613	print*,rad,cpp,g,r,rcp,daysec
614	PRINT*,'check: tsurf ',tsurf(1,:),tsurf(ngrid,:)
615	PRINT*,'check: tsoil ',tsoil(1,1,:),tsoil(ngrid,nsoilmx,:)
616	PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)
617	PRINT*,'check: midlayer,layer ', mlayer(:),layer(:)
618	PRINT*,'check: tracernames ', noms
619	PRINT*,'check: emis ',emis(1,:),emis(ngrid,:)
620	PRINT*,'check: q2 ',q2(1,1),q2(ngrid,nlayer+1)
621	PRINT*,'check: qsurf ',qsurf(1,1,:),qsurf(ngrid,nq,:)
622	PRINT*,'check: co2ice ',qsurf(1,igcm_co2,:),qsurf(ngrid,igcm_co2,:)
623	!!!
624	day_ini = pday
625	!!! a couple initializations (dummy for mesoscale) done in phyetat0
626	!!! --- maybe this should be done in update_inputs_physiq_mod
627
628	tauscaling(:)=1.0 !! probably important
629	totcloudfrac(:)=1.0
630	albedo(:,1)=albedodat(:)
631	albedo(:,2)=albedo(:,1)
632	watercap(:)=0.0
633	#endif
634	#ifndef MESOSCALE
635	if (.not.startphy_file) then
636	! starting without startfi.nc and with callsoil
637	! is not yet possible as soildepth default is not defined
638	if (callsoil) then
639	! default mlayer distribution, following a power law:
640	! mlayer(k)=lay1alpha*(k-1/2)
641	lay1=2.e-4
642	alpha=2
643	do iloop=0,nsoilmx-1
644	mlayer(iloop)=lay1(alpha*(iloop-0.5))
645	enddo
646	lay1=sqrt(mlayer(0)*mlayer(1))
647	alpha=mlayer(1)/mlayer(0)
648	do iloop=1,nsoilmx
649	layer(iloop)=lay1(alpha*(iloop-1))
650	enddo
651	endif
652	! additionnal "academic" initialization of physics
653	do islope = 1,nslope
654	tsurf(:,islope)=pt(:,1)
655	enddo
656	write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
657	do isoil=1,nsoilmx
658	tsoil(1:ngrid,isoil,:)=tsurf(1:ngrid,:)
659	enddo
660	write(,) "Physiq: initializing inertiedat !!"
661	inertiedat(:,:)=400.
662	write(,) "Physiq: initializing day_ini to pdat !"
663	day_ini=pday
664	endif
665	DO islope = 1,nslope
666	inertiedat_tmp(:,:,islope) = inertiedat(:,:)
667	ENDDO
668	#endif
669	if (pday.ne.day_ini) then
670	write(,) "PHYSIQ: ERROR: bad synchronization between ",
671	& "physics and dynamics"
672	write(,) "dynamics day [pday]: ",pday
673	write(,) "physics day [day_ini]: ",day_ini
674	call abort_physic("physiq","dynamics day /= physics day",1)
675	endif
676
677	write (,) 'In physiq day_ini =', day_ini
678
679	c initialize tracers
680	c ~~~~~~~~~~~~~~~~~~
681
682	CALL initracer(ngrid,nq,qsurf)
683
684	c Initialize albedo and orbital calculation
685	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
686	CALL surfini(ngrid,qsurf)
687
688	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
689
690	c initialize soil
691	c ~~~~~~~~~~~~~~~
692	IF (callsoil) THEN
693	c Thermal inertia feedback:
694	IF (tifeedback) THEN
695	DO islope = 1,nslope
696	CALL soil_tifeedback(ngrid,nsoilmx,
697	s qsurf(:,:,islope),inertiesoil(:,:,islope))
698	ENDDO
699	CALL soil(ngrid,nsoilmx,firstcall,inertiesoil,
700	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
701	ELSE
702	CALL soil(ngrid,nsoilmx,firstcall,inertiedat_tmp,
703	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
704	ENDIF ! of IF (tifeedback)
705	ELSE
706	PRINT*,
707	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
708	DO ig=1,ngrid
709	capcal(ig,:)=1.e5
710	fluxgrd(ig,:)=0.
711	ENDDO
712	ENDIF
713	icount=1
714
715	#ifndef MESOSCALE
716	c Initialize thermospheric parameters
717	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
718
719	if (callthermos) then
720	call fill_data_thermos
721	call allocate_param_thermos(nlayer)
722	call allocate_param_iono(nlayer,nreact)
723	call param_read_e107
724	endif
725	#endif
726	c Initialize R and Cp as constant
727
728	if (.not.callthermos .and. .not.photochem) then
729	do l=1,nlayer
730	do ig=1,ngrid
731	rnew(ig,l)=r
732	cpnew(ig,l)=cpp
733	mmean(ig,l)=mugaz
734	enddo
735	enddo
736	endif
737
738	if(callnlte.and.nltemodel.eq.2) call nlte_setup
739	if(callnirco2.and.nircorr.eq.1) call NIR_leedat
740
741
742	IF (water.AND.(ngrid.NE.1)) THEN
743	write(,)"physiq: water_param Surface water frost albedo:",
744	. albedo_h2o_frost
745	write(,)"physiq: water_param Surface watercap albedo:",
746	. albedo_h2o_cap
747	ENDIF
748
749	#ifndef MESOSCALE
750	if(ngrid.ne.1) then
751	! no need to compute slopes when in 1D; it is an input
752	if (callslope) call getslopes(ngrid,phisfi)
753	endif !1D
754	if (ecritstart.GT.0) then
755	call physdem0("restartfi.nc",longitude,latitude,
756	& nsoilmx,ngrid,nlayer,nq,
757	& ptimestep,pday,0.,cell_area,
758	& albedodat,inertiedat,def_slope,
759	& subslope_dist)
760	else
761	call physdem0("restartfi.nc",longitude,latitude,
762	& nsoilmx,ngrid,nlayer,nq,
763	& ptimestep,float(day_end),0.,cell_area,
764	& albedodat,inertiedat,def_slope,
765	& subslope_dist)
766	endif
767
768	c Initialize mountain mesh fraction for the entrainment by top flows param.
769	c ~~~~~~~~~~~~~~~
770	if (topflows) call topmons_setup(ngrid)
771
772	#endif
773
774	#ifdef CPP_XIOS
775	! XIOS outputs
776	write(,) "physiq firstcall: call initialize_xios_output"
777	call initialize_xios_output(pday,ptime,ptimestep,daysec,
778	& presnivs,pseudoalt,mlayer)
779	#endif
780	ENDIF ! (end of "if firstcall")
781
782	if (check_physics_inputs) then
783	! Check the validity of input fields coming from the dynamics
784	call check_physics_fields("begin physiq:",pt,pu,pv,pplev,pq)
785	endif
786
787	c ---------------------------------------------------
788	c 1.2 Initializations done at every physical timestep:
789	c ---------------------------------------------------
790	c
791
792	#ifdef CPP_XIOS
793	! update XIOS time/calendar
794	call update_xios_timestep
795	#endif
796
797
798
799	c Initialize various variables
800	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
801	pdv(:,:)=0
802	pdu(:,:)=0
803	pdt(:,:)=0
804	pdq(:,:,:)=0
805	pdpsrf(:)=0
806	zflubid(:,:)=0
807	zdtsurf(:,:)=0
808	dqsurf(:,:,:)=0
809	dsodust(:,:)=0.
810	dsords(:,:)=0.
811	dsotop(:,:)=0.
812	dwatercap(:,:)=0
813
814	call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf,
815	& albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg)
816
817	! Dust scenario conversion coefficient from IRabs to VISext
818	IRtoVIScoef(1:ngrid)=2.6 ! initialized with former value from Montabone et al 2015
819	! recomputed in aeropacity if reff_driven_IRtoVIS_scenario=.true.
820
821	#ifdef DUSTSTORM
822	pq_tmp(:,:,:)=0
823	#endif
824	igout=ngrid/2+1
825
826
827	zday=pday+ptime ! compute time, in sols (and fraction thereof)
828	! Compute local time at each grid point
829	DO ig=1,ngrid
830	local_time(ig)=modulo(1.+(zday-INT(zday))
831	& +(longitude_deg(ig)/15)/24,1.)
832	ENDDO
833	c Compute Solar Longitude (Ls) :
834	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
835	if (season) then
836	call solarlong(zday,zls)
837	else
838	call solarlong(float(day_ini),zls)
839	end if
840
841	c Initialize pressure levels
842	c ~~~~~~~~~~~~~~~~~~
843	zplev(:,:) = pplev(:,:)
844	zplay(:,:) = pplay(:,:)
845	ps(:) = pplev(:,1)
846
847	c Compute geopotential at interlayers
848	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
849	c ponderation des altitudes au niveau des couches en dp/p
850
851	DO l=1,nlayer
852	DO ig=1,ngrid
853	zzlay(ig,l)=pphi(ig,l)/g
854	ENDDO
855	ENDDO
856	DO ig=1,ngrid
857	zzlev(ig,1)=0.
858	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
859	ENDDO
860	DO l=2,nlayer
861	DO ig=1,ngrid
862	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
863	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
864	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
865	ENDDO
866	ENDDO
867
868
869	! Potential temperature calculation not the same in physiq and dynamic
870
871	c Compute potential temperature
872	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
873	DO l=1,nlayer
874	DO ig=1,ngrid
875	zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp
876	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
877	ENDDO
878	ENDDO
879
880	#ifndef MESOSCALE
881	c-----------------------------------------------------------------------
882	c 1.2.5 Compute mean mass, cp, and R
883	c --------------------------------
884
885	if(photochem.or.callthermos) then
886	call concentrations(ngrid,nlayer,nq,
887	& zplay,pt,pdt,pq,pdq,ptimestep)
888	endif
889	#endif
890
891	! Compute vertical velocity (m/s) from vertical mass flux
892	! w = F / (rhoarea) and rho = P/(rT)
893	! but first linearly interpolate mass flux to mid-layers
894	do l=1,nlayer-1
895	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
896	enddo
897	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
898	do l=1,nlayer
899	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) /
900	& (pplay(1:ngrid,l)*cell_area(1:ngrid))
901	! NB: here we use r and not rnew since this diagnostic comes
902	! from the dynamics
903	enddo
904
905	! test for co2 conservation with co2 microphysics
906	if (igcm_co2_ice.ne.0) then
907	! calculates the amount of co2 at the beginning of physics
908	co2totA = 0.
909	do ig=1,ngrid
910	do l=1,nlayer
911	co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g*
912	& (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice)
913	& +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
914	end do
915	do islope = 1,nslope
916	co2totA = co2totA + qsurf(ig,igcm_co2,islope)*
917	& subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
918	enddo
919	end do
920	else
921	co2totA = 0.
922	do ig=1,ngrid
923	do l=1,nlayer
924	co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g*
925	& (pq(ig,l,igcm_co2)
926	& +pdq(ig,l,igcm_co2)*ptimestep)
927	end do
928	do islope = 1,nslope
929	co2totA = co2totA + qsurf(ig,igcm_co2,islope)*
930	& subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
931	enddo
932	end do
933	endif ! of if (igcm_co2_ice.ne.0)
934	c-----------------------------------------------------------------------
935	c 2. Compute radiative tendencies :
936	c------------------------------------
937
938	IF (callrad) THEN
939
940	c Local Solar zenith angle
941	c ~~~~~~~~~~~~~~~~~~~~~~~~
942
943	CALL orbite(zls,dist_sol,declin)
944
945	IF (diurnal) THEN
946	ztim1=SIN(declin)
947	ztim2=COS(declin)COS(2.pi*(zday-.5))
948	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
949
950	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
951	& ztim1,ztim2,ztim3, mu0,fract)
952
953	ELSE
954	CALL mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad)
955	ENDIF ! of IF (diurnal)
956
957	IF( MOD(icount-1,iradia).EQ.0) THEN
958
959	c NLTE cooling from CO2 emission
960	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
961	IF(callnlte) then
962	if(nltemodel.eq.0.or.nltemodel.eq.1) then
963	CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte)
964	else if(nltemodel.eq.2) then
965	co2vmr_gcm(1:ngrid,1:nlayer)=
966	& pq(1:ngrid,1:nlayer,igcm_co2)*
967	& mmean(1:ngrid,1:nlayer)/mmol(igcm_co2)
968	n2vmr_gcm(1:ngrid,1:nlayer)=
969	& pq(1:ngrid,1:nlayer,igcm_n2)*
970	& mmean(1:ngrid,1:nlayer)/mmol(igcm_n2)
971	covmr_gcm(1:ngrid,1:nlayer)=
972	& pq(1:ngrid,1:nlayer,igcm_co)*
973	& mmean(1:ngrid,1:nlayer)/mmol(igcm_co)
974	ovmr_gcm(1:ngrid,1:nlayer)=
975	& pq(1:ngrid,1:nlayer,igcm_o)*
976	& mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
977
978	CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
979	$ pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
980	$ ovmr_gcm, zdtnlte,ierr_nlte,varerr )
981	if(ierr_nlte.gt.0) then
982	write(,)
983	$ 'WARNING: nlte_tcool output with error message',
984	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
985	write(,)'I will continue anyway'
986	endif
987
988	zdtnlte(1:ngrid,1:nlayer)=
989	& zdtnlte(1:ngrid,1:nlayer)/86400.
990	endif
991	ELSE
992	zdtnlte(:,:)=0.
993	ENDIF !end callnlte
994
995	c Find number of layers for LTE radiation calculations
996	IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
997	& CALL nlthermeq(ngrid,nlayer,zplev,zplay)
998
999	c rocketstorm : compute dust storm mesh fraction
1000	IF (rdstorm) THEN
1001	CALL calcstormfract(ngrid,nlayer,nq,pq,
1002	& totstormfract)
1003	ENDIF
1004
1005	c Note: Dustopacity.F has been transferred to callradite.F
1006
1007	#ifdef DUSTSTORM
1008	!! specific case: save the quantity of dust before adding perturbation
1009
1010	if (firstcall) then
1011	pq_tmp(1:ngrid,1:nlayer,1)=pq(1:ngrid,1:nlayer,igcm_dust_mass)
1012	pq_tmp(1:ngrid,1:nlayer,2)=pq(1:ngrid,1:nlayer,igcm_dust_number)
1013	endif
1014	#endif
1015
1016	c Call main radiative transfer scheme
1017	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1018	c Transfer through CO2 (except NIR CO2 absorption)
1019	c and aerosols (dust and water ice)
1020	! callradite for background dust (out of the rdstorm fraction)
1021	clearatm=.true.
1022	!! callradite for background dust (out of the topflows fraction)
1023	nohmons=.true.
1024
1025	c Radiative transfer
1026	c ------------------
1027	! callradite for the part with clouds
1028	clearsky=.false. ! part with clouds for both cases CLFvarying true/false
1029	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
1030	& albedo_meshavg,emis_meshavg,
1031	& mu0,zplev,zplay,pt,tsurf_meshavg,fract,dist_sol,igout,
1032	& zdtlw,zdtsw,fluxsurf_lw(:,iflat),fluxsurf_dn_sw(:,:,iflat),
1033	& fluxsurf_up_sw,
1034	& fluxtop_lw,fluxtop_dn_sw,fluxtop_up_sw,
1035	& tau_pref_scenario,tau_pref_gcm,
1036	& tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
1037	& taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,
1038	& qsurf_meshavg(:,igcm_co2),rstormdust,rtopdust,totstormfract,
1039	& clearatm,dsords,dsotop,nohmons,clearsky,totcloudfrac)
1040
1041	DO islope=1,nslope
1042	fluxsurf_lw(:,islope) =fluxsurf_lw(:,iflat)
1043	fluxsurf_dn_sw(:,:,islope) =fluxsurf_dn_sw(:,:,iflat)
1044	ENDDO
1045
1046	! case of sub-grid water ice clouds: callradite for the clear case
1047	IF (CLFvarying) THEN
1048	! ---> PROBLEMS WITH ALLOCATED ARRAYS
1049	! (temporary solution in callcorrk: do not deallocate
1050	! if
1051	! CLFvarying ...) ?? AP ??
1052	clearsky=.true.
1053	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
1054	& albedo_meshavg,emis_meshavg,mu0,zplev,zplay,pt,
1055	& tsurf_meshavg,fract,
1056	& dist_sol,igout,zdtlwclf,zdtswclf,
1057	& fluxsurf_lwclf,fluxsurf_dn_swclf,fluxsurf_up_swclf,
1058	& fluxtop_lwclf,fluxtop_dn_swclf,fluxtop_up_swclf,
1059	& tau_pref_scenario,tau_pref_gcm,tau,aerosol,
1060	& dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
1061	& taucloudtesclf,rdust,
1062	& rice,nuice,riceco2, nuiceco2,
1063	& qsurf_meshavg(:,igcm_co2),
1064	& rstormdust,rtopdust,totstormfract,
1065	& clearatm,dsords,dsotop,
1066	& nohmons,clearsky,totcloudfrac)
1067	clearsky = .false. ! just in case.
1068	! Sum the fluxes and heating rates from cloudy/clear
1069	! cases
1070	DO ig=1,ngrid
1071	tf_clf=totcloudfrac(ig)
1072	ntf_clf=1.-tf_clf
1073	DO islope=1,nslope
1074	fluxsurf_lw(ig,islope) = ntf_clf*fluxsurf_lwclf(ig)
1075	& + tf_clf*fluxsurf_lw(ig,islope)
1076	fluxsurf_dn_sw(ig,1:2,islope) =
1077	& ntf_clf*fluxsurf_dn_swclf(ig,1:2)
1078	& + tf_clf*fluxsurf_dn_sw(ig,1:2,islope)
1079	ENDDO
1080	fluxsurf_up_sw(ig,1:2) =
1081	& ntf_clf*fluxsurf_up_swclf(ig,1:2)
1082	& + tf_clf*fluxsurf_up_sw(ig,1:2)
1083	fluxtop_lw(ig) = ntf_clf*fluxtop_lwclf(ig)
1084	& + tf_clf*fluxtop_lw(ig)
1085	fluxtop_dn_sw(ig,1:2)=ntf_clf*fluxtop_dn_swclf(ig,1:2)
1086	& + tf_clf*fluxtop_dn_sw(ig,1:2)
1087	fluxtop_up_sw(ig,1:2)=ntf_clf*fluxtop_up_swclf(ig,1:2)
1088	& + tf_clf*fluxtop_up_sw(ig,1:2)
1089	taucloudtes(ig) = ntf_clf*taucloudtesclf(ig)
1090	& + tf_clf*taucloudtes(ig)
1091	zdtlw(ig,1:nlayer) = ntf_clf*zdtlwclf(ig,1:nlayer)
1092	& + tf_clf*zdtlw(ig,1:nlayer)
1093	zdtsw(ig,1:nlayer) = ntf_clf*zdtswclf(ig,1:nlayer)
1094	& + tf_clf*zdtsw(ig,1:nlayer)
1095	ENDDO
1096
1097	ENDIF ! (CLFvarying)
1098
1099	!============================================================================
1100
1101	#ifdef DUSTSTORM
1102	!! specific case: compute the added quantity of dust for perturbation
1103	if (firstcall) then
1104	pdq(1:ngrid,1:nlayer,igcm_dust_mass)=
1105	& pdq(1:ngrid,1:nlayer,igcm_dust_mass)
1106	& - pq_tmp(1:ngrid,1:nlayer,1)
1107	& + pq(1:ngrid,1:nlayer,igcm_dust_mass)
1108	pdq(1:ngrid,1:nlayer,igcm_dust_number)=
1109	& pdq(1:ngrid,1:nlayer,igcm_dust_number)
1110	& - pq_tmp(1:ngrid,1:nlayer,2)
1111	& + pq(1:ngrid,1:nlayer,igcm_dust_number)
1112	endif
1113	#endif
1114
1115	c Outputs for basic check (middle of domain)
1116	c ------------------------------------------
1117	write(*,'("Ls =",f11.6," check lat =",f10.6,
1118	& " lon =",f11.6)')
1119	& zls180./pi,latitude(igout)180/pi,
1120	& longitude(igout)*180/pi
1121
1122	write(*,'(" tau_pref_gcm(",f4.0," Pa) =",f9.6,
1123	& " tau(",f4.0," Pa) =",f9.6)')
1124	& odpref,tau_pref_gcm(igout),
1125	& odpref,tau(igout,1)*odpref/zplev(igout,1)
1126
1127
1128	c ---------------------------------------------------------
1129	c Call slope parameterization for direct and scattered flux
1130	c ---------------------------------------------------------
1131	IF(callslope) THEN
1132	! assume that in this case, nslope = 1
1133	if(nslope.ne.1) then
1134	call abort_physic(
1135	& "physiq","callslope=true but nslope.ne.1",1)
1136	endif
1137	print *, 'Slope scheme is on and computing...'
1138	DO ig=1,ngrid
1139	sl_the = theta_sl(ig)
1140	IF (sl_the .ne. 0.) THEN
1141	ztim1=fluxsurf_dn_sw(ig,1,1)+fluxsurf_dn_sw(ig,2,1)
1142	DO l=1,2
1143	sl_lct = ptime24. + 180.longitude(ig)/pi/15.
1144	sl_ra = pi*(1.0-sl_lct/12.)
1145	sl_lat = 180.*latitude(ig)/pi
1146	sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
1147	sl_alb = albedo(ig,l,1)
1148	sl_psi = psi_sl(ig)
1149	sl_fl0 = fluxsurf_dn_sw(ig,l,1)
1150	sl_di0 = 0.
1151	if ((mu0(ig) .gt. 0.).and.(ztim1.gt.0.)) then
1152	sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
1153	sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
1154	sl_di0 = sl_di0/ztim1
1155	sl_di0 = fluxsurf_dn_sw(ig,l,1)*sl_di0
1156	endif
1157	! you never know (roundup concern...)
1158	if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
1159	!!!!!!!!!!!!!!!!!!!!!!!!!!
1160	CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
1161	& sl_tau, sl_alb, sl_the, sl_psi,
1162	& sl_di0, sl_fl0, sl_flu )
1163	!!!!!!!!!!!!!!!!!!!!!!!!!!
1164	fluxsurf_dn_sw(ig,l,1) = sl_flu
1165	ENDDO
1166	!!! compute correction on IR flux as well
1167	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
1168	fluxsurf_lw(ig,:)= fluxsurf_lw(ig,:)*sky
1169	ENDIF
1170	ENDDO
1171	c ELSE ! not calling subslope, nslope might be > 1
1172	c DO islope = 1,nslope
1173	c IF(def_slope_mean(islope).ge.0.) THEN
1174	c sl_psi = 0. !Northward slope
1175	c ELSE
1176	c sl_psi = 180. !Southward slope
1177	c ENDIF
1178	c sl_the=abs(def_slope_mean(islope))
1179	c IF (sl_the .gt. 1e-6) THEN
1180	c DO ig=1,ngrid
1181	c ztim1=fluxsurf_dn_sw(ig,1,islope)
1182	c s +fluxsurf_dn_sw(ig,2,islope)
1183	c DO l=1,2
1184	c sl_lct = ptime24. + 180.longitude(ig)/pi/15.
1185	c sl_ra = pi*(1.0-sl_lct/12.)
1186	c sl_lat = 180.*latitude(ig)/pi
1187	c sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
1188	c sl_alb = albedo(ig,l,islope)
1189	c sl_psi = psi_sl(ig)
1190	c sl_fl0 = fluxsurf_dn_sw(ig,l,islope)
1191	c sl_di0 = 0.
1192	c if (mu0(ig) .gt. 0.) then
1193	c sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
1194	c sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
1195	c sl_di0 = sl_di0/ztim1
1196	c sl_di0 = fluxsurf_dn_sw(ig,l,islope)*sl_di0
1197	c endif
1198	c ! you never know (roundup concern...)
1199	c if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
1200	c !!!!!!!!!!!!!!!!!!!!!!!!!!
1201	c CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
1202	c & sl_tau, sl_alb, sl_the, sl_psi,
1203	c & sl_di0, sl_fl0, sl_flu )
1204	c !!!!!!!!!!!!!!!!!!!!!!!!!!
1205	c fluxsurf_dn_sw(ig,l,islope) = sl_flu
1206	c ENDDO
1207	c !!! compute correction on IR flux as well
1208	c
1209	c fluxsurf_lw(ig,islope)= fluxsurf_lw(ig,islope)
1210	c & *sky_slope(islope)
1211	c ENDDO
1212	c ENDIF
1213	c ENDDO ! islope = 1,nslope
1214	ENDIF ! callslope
1215
1216	c CO2 near infrared absorption
1217	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1218	zdtnirco2(:,:)=0
1219	if (callnirco2) then
1220	call nirco2abs (ngrid,nlayer,zplay,dist_sol,nq,pq,
1221	. mu0,fract,declin, zdtnirco2)
1222	endif
1223
1224	c Radiative flux from the sky absorbed by the surface (W.m-2)
1225	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1226	DO ig=1,ngrid
1227	DO islope = 1,nslope
1228	fluxrad_sky(ig,islope) =
1229	$ emis(ig,islope)*fluxsurf_lw(ig,islope)
1230	$ +fluxsurf_dn_sw(ig,1,islope)*(1.-albedo(ig,1,islope))
1231	$ +fluxsurf_dn_sw(ig,2,islope)*(1.-albedo(ig,2,islope))
1232	ENDDO
1233	ENDDO
1234
1235
1236	c Net atmospheric radiative heating rate (K.s-1)
1237	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1238	IF(callnlte) THEN
1239	CALL blendrad(ngrid, nlayer, zplay,
1240	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
1241	ELSE
1242	DO l=1,nlayer
1243	DO ig=1,ngrid
1244	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
1245	& +zdtnirco2(ig,l)
1246	ENDDO
1247	ENDDO
1248	ENDIF
1249
1250	ENDIF ! of if(mod(icount-1,iradia).eq.0)
1251
1252	c Transformation of the radiative tendencies:
1253	c -------------------------------------------
1254
1255	c Net radiative surface flux (W.m-2)
1256	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
1257	c
1258	DO ig=1,ngrid
1259	DO islope = 1,nslope
1260	zplanck(ig)=tsurf(ig,islope)*tsurf(ig,islope)
1261	zplanck(ig)=emis(ig,islope)*
1262	$ stephanzplanck(ig)zplanck(ig)
1263	fluxrad(ig,islope)=fluxrad_sky(ig,islope)-zplanck(ig)
1264	IF(callslope) THEN
1265	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
1266	fluxrad(ig,nslope)=fluxrad(ig,nslope)+
1267	$ (1.-sky)*zplanck(ig)
1268	ELSE
1269	fluxrad(ig,islope)=fluxrad(ig,islope) +
1270	$ (1.-sky_slope(iflat))emis(ig,iflat)
1271	$ stephantsurf(ig,iflat)*4
1272	ENDIF
1273	ENDDO
1274	ENDDO
1275
1276	DO l=1,nlayer
1277	DO ig=1,ngrid
1278	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
1279	ENDDO
1280	ENDDO
1281
1282	ENDIF ! of IF (callrad)
1283
1284	c 3.1 Rocket dust storm
1285	c -------------------------------------------
1286	IF (rdstorm) THEN
1287	clearatm=.false.
1288	pdqrds(:,:,:)=0.
1289	qdusttotal0(:)=0.
1290	qdusttotal1(:)=0.
1291	do ig=1,ngrid
1292	do l=1,nlayer
1293	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) ! updated potential
1294	! temperature tendency
1295	! for diagnostics
1296	qdustrds0(ig,l)=pq(ig,l,igcm_dust_mass)+
1297	& pdq(ig,l,igcm_dust_mass)*ptimestep
1298	qstormrds0(ig,l)=pq(ig,l,igcm_stormdust_mass)+
1299	& pdq(ig,l,igcm_stormdust_mass)*ptimestep
1300	qdusttotal0(ig)=qdusttotal0(ig)+(qdustrds0(ig,l)+
1301	& qstormrds0(ig,l))*(zplev(ig,l)-
1302	& zplev(ig,l+1))/g
1303	enddo
1304	enddo
1305	call write_output('qdustrds0','qdust before rds',
1306	& 'kg/kg ',qdustrds0(:,:))
1307	call write_output('qstormrds0','qstorm before rds',
1308	& 'kg/kg ',qstormrds0(:,:))
1309
1310	CALL rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep,
1311	& pq,pdq,pt,pdt,zplev,zplay,zzlev,
1312	& zzlay,zdtsw,zdtlw,
1313	c for radiative transfer
1314	& clearatm,icount,zday,zls,
1315	& tsurf,qsurf_meshavg(:,igcm_co2),igout,
1316	& totstormfract,tauscaling,
1317	& dust_rad_adjust,IRtoVIScoef,
1318	c input sub-grid scale cloud
1319	& clearsky,totcloudfrac,
1320	c input sub-grid scale topography
1321	& nohmons,
1322	c output
1323	& pdqrds,wspeed,dsodust,dsords,dsotop,
1324	& tau_pref_scenario,tau_pref_gcm)
1325
1326	c update the tendencies of both dust after vertical transport
1327	DO l=1,nlayer
1328	DO ig=1,ngrid
1329	pdq(ig,l,igcm_stormdust_mass)=
1330	& pdq(ig,l,igcm_stormdust_mass)+
1331	& pdqrds(ig,l,igcm_stormdust_mass)
1332	pdq(ig,l,igcm_stormdust_number)=
1333	& pdq(ig,l,igcm_stormdust_number)+
1334	& pdqrds(ig,l,igcm_stormdust_number)
1335
1336	pdq(ig,l,igcm_dust_mass)=
1337	& pdq(ig,l,igcm_dust_mass)+ pdqrds(ig,l,igcm_dust_mass)
1338	pdq(ig,l,igcm_dust_number)=
1339	& pdq(ig,l,igcm_dust_number)+
1340	& pdqrds(ig,l,igcm_dust_number)
1341
1342	ENDDO
1343	ENDDO
1344	do l=1,nlayer
1345	do ig=1,ngrid
1346	qdustrds1(ig,l)=pq(ig,l,igcm_dust_mass)+
1347	& pdq(ig,l,igcm_dust_mass)*ptimestep
1348	qstormrds1(ig,l)=pq(ig,l,igcm_stormdust_mass)+
1349	& pdq(ig,l,igcm_stormdust_mass)*ptimestep
1350	qdusttotal1(ig)=qdusttotal1(ig)+(qdustrds1(ig,l)+
1351	& qstormrds1(ig,l))*(zplev(ig,l)-
1352	& zplev(ig,l+1))/g
1353	enddo
1354	enddo
1355
1356	c for diagnostics
1357	call write_output('qdustrds1','qdust after rds',
1358	& 'kg/kg ',qdustrds1(:,:))
1359	call write_output('qstormrds1','qstorm after rds',
1360	& 'kg/kg ',qstormrds1(:,:))
1361
1362	call write_output('qdusttotal0','q sum before rds',
1363	& 'kg/m2 ',qdusttotal0(:))
1364	call write_output('qdusttotal1','q sum after rds',
1365	& 'kg/m2 ',qdusttotal1(:))
1366
1367	ENDIF ! end of if(rdstorm)
1368
1369	c 3.2 Dust entrained from the PBL up to the top of sub-grid scale topography
1370	c -------------------------------------------
1371	IF (topflows) THEN
1372	clearatm=.true. ! stormdust is not accounted in the extra heating on top of the mountains
1373	nohmons=.false.
1374	pdqtop(:,:,:)=0.
1375	CALL topmons(ngrid,nlayer,nq,ptime,ptimestep,
1376	& pq,pdq,pt,pdt,zplev,zplay,zzlev,
1377	& zzlay,zdtsw,zdtlw,
1378	& icount,zday,zls,tsurf(:,iflat),
1379	& qsurf_meshavg(:,igcm_co2),
1380	& igout,aerosol,tauscaling,dust_rad_adjust,
1381	& IRtoVIScoef,totstormfract,clearatm,
1382	& clearsky,totcloudfrac,
1383	& nohmons,
1384	& pdqtop,wtop,dsodust,dsords,dsotop,
1385	& tau_pref_scenario,tau_pref_gcm)
1386
1387	c update the tendencies of both dust after vertical transport
1388	DO l=1,nlayer
1389	DO ig=1,ngrid
1390	pdq(ig,l,igcm_topdust_mass)=
1391	& pdq(ig,l,igcm_topdust_mass)+
1392	& pdqtop(ig,l,igcm_topdust_mass)
1393	pdq(ig,l,igcm_topdust_number)=
1394	& pdq(ig,l,igcm_topdust_number)+
1395	& pdqtop(ig,l,igcm_topdust_number)
1396	pdq(ig,l,igcm_dust_mass)=
1397	& pdq(ig,l,igcm_dust_mass)+ pdqtop(ig,l,igcm_dust_mass)
1398	pdq(ig,l,igcm_dust_number)=
1399	& pdq(ig,l,igcm_dust_number)+pdqtop(ig,l,igcm_dust_number)
1400
1401	ENDDO
1402	ENDDO
1403
1404	ENDIF ! end of if (topflows)
1405
1406	c 3.3 Dust injection from the surface
1407	c -------------------------------------------
1408	if (dustinjection.gt.0) then
1409
1410	CALL compute_dtau(ngrid,nlayer,
1411	& zday,pplev,tau_pref_gcm,
1412	& ptimestep,local_time,IRtoVIScoef,
1413	& dustliftday)
1414	endif ! end of if (dustinjection.gt.0)
1415
1416
1417	c-----------------------------------------------------------------------
1418	c 4. Gravity wave and subgrid scale topography drag :
1419	c -------------------------------------------------
1420
1421
1422	IF(calllott)THEN
1423	CALL calldrag_noro(ngrid,nlayer,ptimestep,
1424	& zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
1425
1426	DO l=1,nlayer
1427	DO ig=1,ngrid
1428	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
1429	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
1430	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
1431	ENDDO
1432	ENDDO
1433	ENDIF
1434
1435	c-----------------------------------------------------------------------
1436	c 5. Vertical diffusion (turbulent mixing):
1437	c -----------------------------------------
1438
1439	IF (calldifv) THEN
1440	DO ig=1,ngrid
1441	DO islope = 1,nslope
1442	zflubid(ig,islope)=fluxrad(ig,islope)
1443	& +fluxgrd(ig,islope)
1444	ENDDO
1445	ENDDO
1446	zdum1(:,:)=0
1447	zdum2(:,:)=0
1448	do l=1,nlayer
1449	do ig=1,ngrid
1450	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
1451	enddo
1452	enddo
1453
1454	c ----------------------
1455	c Treatment of a special case : using new surface layer (Richardson based)
1456	c without using the thermals in gcm and mesoscale can yield problems in
1457	c weakly unstable situations when winds are near to 0. For those cases, we add
1458	c a unit subgrid gustiness. Remember that thermals should be used we using the
1459	c Richardson based surface layer model.
1460	IF ( .not.calltherm
1461	. .and. callrichsl
1462	. .and. .not.turb_resolved) THEN
1463
1464	DO ig=1, ngrid
1465	IF (zh(ig,1) .lt. tsurf_meshavg(ig)) THEN
1466	wstar(ig)=1.
1467	hfmax_th(ig)=0.2
1468	ELSE
1469	wstar(ig)=0.
1470	hfmax_th(ig)=0.
1471	ENDIF
1472	ENDDO
1473	ENDIF
1474	c ----------------------
1475
1476	IF (tke_heat_flux .ne. 0.) THEN
1477
1478	zz1(:)=(pt(:,1)+pdt(:,1)ptimestep)(r/g)*
1479	& (-alog(zplay(:,1)/zplev(:,1)))
1480	pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
1481	ENDIF
1482
1483	c Calling vdif (Martian version WITH CO2 condensation)
1484	dwatercap_dif(:,:) = 0.
1485	zcdh(:) = 0.
1486	zcdv(:) = 0.
1487
1488	CALL vdifc(ngrid,nlayer,nq,zpopsk,
1489	$ ptimestep,capcal,lwrite,
1490	$ zplay,zplev,zzlay,zzlev,z0,
1491	$ pu,pv,zh,pq,tsurf,emis,qsurf,
1492	$ zdum1,zdum2,zdh,pdq,zflubid,
1493	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
1494	& zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,
1495	& zcondicea_co2microp,sensibFlux,
1496	& dustliftday,local_time,watercap,dwatercap_dif)
1497
1498	DO ig=1,ngrid
1499	zdtsurf(ig,:)=zdtsurf(ig,:)+zdtsdif(ig,:)
1500	dwatercap(ig,:)=dwatercap(ig,:)+dwatercap_dif(ig,:)
1501	ENDDO
1502	call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,zdqsdif,
1503	& albedo_meshavg,emis_meshavg,tsurf_meshavg,zdqsdif_meshavg_tmp)
1504	IF (.not.turb_resolved) THEN
1505	DO l=1,nlayer
1506	DO ig=1,ngrid
1507	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
1508	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
1509	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
1510
1511	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
1512	ENDDO
1513	ENDDO
1514
1515	DO iq=1, nq
1516	DO l=1,nlayer
1517	DO ig=1,ngrid
1518	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
1519	ENDDO
1520	ENDDO
1521	ENDDO
1522	DO iq=1, nq
1523	DO ig=1,ngrid
1524	dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:)
1525	ENDDO
1526	ENDDO
1527
1528	ELSE
1529	write (,) '******************************************'
1530	write (,) '** LES mode: the difv part is only used to'
1531	write (,) '** - provide HFX and UST to the dynamics'
1532	write (,) '** - update TSURF'
1533	write (,) '******************************************'
1534	!! Specific treatment for lifting in turbulent-resolving mode (AC)
1535	IF (lifting .and. doubleq) THEN
1536	!! lifted dust is injected in the first layer.
1537	!! Sedimentation must be called after turbulent mixing, i.e. on next step, after WRF.
1538	!! => lifted dust is not incremented before the sedimentation step.
1539	zdqdif(1:ngrid,1,1:nq)=0.
1540	zdqdif(1:ngrid,1,igcm_dust_number) =
1541	. -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_number)
1542	zdqdif(1:ngrid,1,igcm_dust_mass) =
1543	. -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_mass)
1544	zdqdif(1:ngrid,2:nlayer,1:nq) = 0.
1545	DO iq=1, nq
1546	IF ((iq .ne. igcm_dust_mass)
1547	& .and. (iq .ne. igcm_dust_number)) THEN
1548	zdqsdif(:,iq,:)=0.
1549	ENDIF
1550	ENDDO
1551	ELSE
1552	zdqdif(1:ngrid,1:nlayer,1:nq) = 0.
1553	zdqsdif(1:ngrid,1:nq,1:nslope) = 0.
1554	ENDIF
1555	ENDIF
1556	ELSE
1557	DO ig=1,ngrid
1558	DO islope=1,nslope
1559	zdtsurf(ig,islope)=zdtsurf(ig,islope)+
1560	s (fluxrad(ig,islope)+fluxgrd(ig,islope))/capcal(ig,islope)
1561	ENDDO
1562	ENDDO
1563	IF (turb_resolved) THEN
1564	write(,) 'Turbulent-resolving mode !'
1565	write(,) 'Please set calldifv to T in callphys.def'
1566	call abort_physic("physiq","turbulent-resolving mode",1)
1567	ENDIF
1568	ENDIF ! of IF (calldifv)
1569
1570	c-----------------------------------------------------------------------
1571	c 6. Thermals :
1572	c -----------------------------
1573
1574	if(calltherm .and. .not.turb_resolved) then
1575
1576	call calltherm_interface(ngrid,nlayer,nq,igcm_co2,
1577	$ zzlev,zzlay,
1578	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
1579	$ zplay,zplev,pphi,zpopsk,
1580	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
1581	$ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
1582
1583	DO l=1,nlayer
1584	DO ig=1,ngrid
1585	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
1586	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
1587	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
1588	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
1589	ENDDO
1590	ENDDO
1591
1592	DO ig=1,ngrid
1593	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
1594	ENDDO
1595
1596	DO iq=1,nq
1597	DO l=1,nlayer
1598	DO ig=1,ngrid
1599	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
1600	ENDDO
1601	ENDDO
1602	ENDDO
1603
1604	lmax_th_out(:)=real(lmax_th(:))
1605
1606	else !of if calltherm
1607	lmax_th(:)=0
1608	wstar(:)=0.
1609	hfmax_th(:)=0.
1610	lmax_th_out(:)=0.
1611	end if
1612
1613	c-----------------------------------------------------------------------
1614	c 7. Dry convective adjustment:
1615	c -----------------------------
1616
1617	IF(calladj) THEN
1618
1619	DO l=1,nlayer
1620	DO ig=1,ngrid
1621	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
1622	ENDDO
1623	ENDDO
1624	zduadj(:,:)=0
1625	zdvadj(:,:)=0
1626	zdhadj(:,:)=0
1627
1628	CALL convadj(ngrid,nlayer,nq,ptimestep,
1629	$ zplay,zplev,zpopsk,lmax_th,
1630	$ pu,pv,zh,pq,
1631	$ pdu,pdv,zdh,pdq,
1632	$ zduadj,zdvadj,zdhadj,
1633	$ zdqadj)
1634
1635	DO l=1,nlayer
1636	DO ig=1,ngrid
1637	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
1638	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
1639	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
1640
1641	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
1642	ENDDO
1643	ENDDO
1644
1645	DO iq=1, nq
1646	DO l=1,nlayer
1647	DO ig=1,ngrid
1648	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
1649	ENDDO
1650	ENDDO
1651	ENDDO
1652	ENDIF ! of IF(calladj)
1653
1654	c-----------------------------------------------------
1655	c 8. Non orographic Gravity waves :
1656	c -------------------------------------------------
1657
1658	IF (calllott_nonoro) THEN
1659
1660	CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep,
1661	& cpnew,rnew,
1662	& zplay,
1663	& zmax_th, ! max altitude reached by thermals (m)
1664	& pt, pu, pv,
1665	& pdt, pdu, pdv,
1666	& zustrhi,zvstrhi,
1667	& d_t_hin, d_u_hin, d_v_hin)
1668
1669	! Update tendencies
1670	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)
1671	& +d_t_hin(1:ngrid,1:nlayer)
1672	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)
1673	& +d_u_hin(1:ngrid,1:nlayer)
1674	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)
1675	& +d_v_hin(1:ngrid,1:nlayer)
1676
1677	ENDIF ! of IF (calllott_nonoro)
1678
1679	c-----------------------------------------------------------------------
1680	c 9. Specific parameterizations for tracers
1681	c: -----------------------------------------
1682
1683
1684	c 9a. Water and ice
1685	c ---------------
1686
1687	c ---------------------------------------
1688	c Water ice condensation in the atmosphere
1689	c ----------------------------------------
1690	IF (water) THEN
1691
1692	call watercloud(ngrid,nlayer,ptimestep,
1693	& zplev,zplay,pdpsrf,zzlay, pt,pdt,
1694	& pq,pdq,zdqcloud,zdtcloud,
1695	& nq,tau,tauscaling,rdust,rice,nuice,
1696	& rsedcloud,rhocloud,totcloudfrac)
1697	c Temperature variation due to latent heat release
1698	if (activice) then
1699	pdt(1:ngrid,1:nlayer) =
1700	& pdt(1:ngrid,1:nlayer) +
1701	& zdtcloud(1:ngrid,1:nlayer)
1702	endif
1703
1704	! increment water vapour and ice atmospheric tracers tendencies
1705	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) =
1706	& pdq(1:ngrid,1:nlayer,igcm_h2o_vap) +
1707	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_vap)
1708	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
1709	& pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
1710	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice)
1711
1712	if (hdo) then
1713	! increment HDO vapour and ice atmospheric tracers tendencies
1714	pdq(1:ngrid,1:nlayer,igcm_hdo_vap) =
1715	& pdq(1:ngrid,1:nlayer,igcm_hdo_vap) +
1716	& zdqcloud(1:ngrid,1:nlayer,igcm_hdo_vap)
1717	pdq(1:ngrid,1:nlayer,igcm_hdo_ice) =
1718	& pdq(1:ngrid,1:nlayer,igcm_hdo_ice) +
1719	& zdqcloud(1:ngrid,1:nlayer,igcm_hdo_ice)
1720	endif !hdo
1721
1722	! increment dust and ccn masses and numbers
1723	! We need to check that we have Nccn & Ndust > 0
1724	! This is due to single precision rounding problems
1725	if (microphys) then
1726	pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
1727	& pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
1728	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
1729	pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
1730	& pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
1731	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
1732	where (pq(:,:,igcm_ccn_mass) +
1733	& ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
1734	pdq(:,:,igcm_ccn_mass) =
1735	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1736	pdq(:,:,igcm_ccn_number) =
1737	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1738	end where
1739	where (pq(:,:,igcm_ccn_number) +
1740	& ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
1741	pdq(:,:,igcm_ccn_mass) =
1742	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1743	pdq(:,:,igcm_ccn_number) =
1744	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1745	end where
1746	endif
1747
1748	if (scavenging) then
1749	pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
1750	& pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
1751	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_mass)
1752	pdq(1:ngrid,1:nlayer,igcm_dust_number) =
1753	& pdq(1:ngrid,1:nlayer,igcm_dust_number) +
1754	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_number)
1755	where (pq(:,:,igcm_dust_mass) +
1756	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
1757	pdq(:,:,igcm_dust_mass) =
1758	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1759	pdq(:,:,igcm_dust_number) =
1760	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1761	end where
1762	where (pq(:,:,igcm_dust_number) +
1763	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
1764	pdq(:,:,igcm_dust_mass) =
1765	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1766	pdq(:,:,igcm_dust_number) =
1767	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1768	end where
1769	endif ! of if scavenging
1770
1771	END IF ! of IF (water)
1772
1773	c 9a bis. CO2 clouds (CL & JA)
1774	c ---------------------------------------
1775	c CO2 ice cloud condensation in the atmosphere
1776	c ----------------------------------------
1777	c flag needed in callphys.def:
1778	c co2clouds=.true. is mandatory (default is .false.)
1779	c co2useh2o=.true. if you want to allow co2 condensation
1780	c on water ice particles
1781	c meteo_flux=.true. if you want to add a meteoritic
1782	c supply of CCN
1783	c CLFvaryingCO2=.true. if you want to have a sub-grid
1784	c temperature distribution
1785	c spantCO2=integer (i.e. 3) amplitude of the sub-grid T disti
1786	c nuiceco2_sed=0.2 variance of the size distribution for the
1787	c sedimentation
1788	c nuiceco2_ref=0.2 variance of the size distribution for the
1789	c nucleation
1790	c imicroco2=50 micro-timestep is 1/50 of physical timestep
1791	zdqssed_co2(:) = 0.
1792	zdqssed_ccn(:,:) = 0.
1793
1794	IF (co2clouds) THEN
1795	call co2cloud(ngrid,nlayer,ptimestep,
1796	& zplev,zplay,pdpsrf,zzlay,pt,pdt,
1797	& pq,pdq,zdqcloudco2,zdtcloudco2,
1798	& nq,tau,tauscaling,rdust,rice,riceco2,nuice,
1799	& rhocloud, rsedcloudco2,rhocloudco2,zzlev,zdqssed_co2,
1800	& zdqssed_ccn,pdu,pu,zcondicea_co2microp)
1801
1802	DO iq=1, nq
1803	DO ig=1,ngrid
1804	DO islope = 1,nslope
1805	dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope)+
1806	& zdqssed_ccn(ig,iq)cos(pidef_slope_mean(islope)/180.)
1807	ENDDO !(islope)
1808	ENDDO ! (ig)
1809	ENDDO ! (iq)q)
1810	c Temperature variation due to latent heat release
1811	pdt(1:ngrid,1:nlayer) =
1812	& pdt(1:ngrid,1:nlayer) +
1813	& zdtcloudco2(1:ngrid,1:nlayer)
1814
1815	! increment dust and ccn masses and numbers
1816	! We need to check that we have Nccn & Ndust > 0
1817	! This is due to single precision rounding problems
1818	! increment dust tracers tendancies
1819	pdq(:,:,igcm_dust_mass) = pdq(:,:,igcm_dust_mass)
1820	& + zdqcloudco2(:,:,igcm_dust_mass)
1821
1822	pdq(:,:,igcm_dust_number) = pdq(:,:,igcm_dust_number)
1823	& + zdqcloudco2(:,:,igcm_dust_number)
1824
1825	pdq(:,:,igcm_co2) = pdq(:,:,igcm_co2)
1826	& + zdqcloudco2(:,:,igcm_co2)
1827
1828	pdq(:,:,igcm_co2_ice) = pdq(:,:,igcm_co2_ice)
1829	& + zdqcloudco2(:,:,igcm_co2_ice)
1830
1831	pdq(:,:,igcm_ccnco2_mass) = pdq(:,:,igcm_ccnco2_mass)
1832	& + zdqcloudco2(:,:,igcm_ccnco2_mass)
1833
1834	pdq(:,:,igcm_ccnco2_number) = pdq(:,:,igcm_ccnco2_number)
1835	& + zdqcloudco2(:,:,igcm_ccnco2_number)
1836
1837	if (meteo_flux) then
1838	pdq(:,:,igcm_ccnco2_meteor_mass) =
1839	& pdq(:,:,igcm_ccnco2_meteor_mass) +
1840	& zdqcloudco2(:,:,igcm_ccnco2_meteor_mass)
1841
1842	pdq(:,:,igcm_ccnco2_meteor_number) =
1843	& pdq(:,:,igcm_ccnco2_meteor_number)
1844	& + zdqcloudco2(:,:,igcm_ccnco2_meteor_number)
1845	end if
1846	!Update water ice clouds values as well
1847	if (co2useh2o) then
1848	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
1849	& pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
1850	& zdqcloudco2(1:ngrid,1:nlayer,igcm_h2o_ice)
1851	pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
1852	& pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
1853	& zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_mass)
1854	pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
1855	& pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
1856	& zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
1857
1858	pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
1859	& pdq(:,:,igcm_ccnco2_h2o_mass_ice) +
1860	& zdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ice)
1861
1862	pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
1863	& pdq(:,:,igcm_ccnco2_h2o_mass_ccn) +
1864	& zdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ccn)
1865
1866	pdq(:,:,igcm_ccnco2_h2o_number) =
1867	& pdq(:,:,igcm_ccnco2_h2o_number) +
1868	& zdqcloudco2(:,:,igcm_ccnco2_h2o_number)
1869
1870	c Negative values?
1871	where (pq(:,:,igcm_ccn_mass) +
1872	& ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
1873	pdq(:,:,igcm_ccn_mass) =
1874	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1875	pdq(:,:,igcm_ccn_number) =
1876	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1877	end where
1878	c Negative values?
1879	where (pq(:,:,igcm_ccn_number) +
1880	& ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
1881	pdq(:,:,igcm_ccn_mass) =
1882	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1883	pdq(:,:,igcm_ccn_number) =
1884	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1885	end where
1886	where (pq(:,:,igcm_ccnco2_h2o_mass_ice) +
1887	& pq(:,:,igcm_ccnco2_h2o_mass_ccn) +
1888	& (pdq(:,:,igcm_ccnco2_h2o_mass_ice) +
1889	& pdq(:,:,igcm_ccnco2_h2o_mass_ccn)
1890	& )*ptimestep < 0.)
1891	pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
1892	& - pq(:,:,igcm_ccnco2_h2o_mass_ice)
1893	& /ptimestep + 1.e-30
1894	pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
1895	& - pq(:,:,igcm_ccnco2_h2o_mass_ccn)
1896	& /ptimestep + 1.e-30
1897	pdq(:,:,igcm_ccnco2_h2o_number) =
1898	& - pq(:,:,igcm_ccnco2_h2o_number)
1899	& /ptimestep + 1.e-30
1900	end where
1901
1902	where (pq(:,:,igcm_ccnco2_h2o_number) +
1903	& (pdq(:,:,igcm_ccnco2_h2o_number)
1904	& )*ptimestep < 0.)
1905	pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
1906	& - pq(:,:,igcm_ccnco2_h2o_mass_ice)
1907	& /ptimestep + 1.e-30
1908	pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
1909	& - pq(:,:,igcm_ccnco2_h2o_mass_ccn)
1910	& /ptimestep + 1.e-30
1911	pdq(:,:,igcm_ccnco2_h2o_number) =
1912	& - pq(:,:,igcm_ccnco2_h2o_number)
1913	& /ptimestep + 1.e-30
1914	end where
1915	endif ! of if (co2useh2o)
1916	c Negative values?
1917	where (pq(:,:,igcm_ccnco2_mass) +
1918	& ptimestep*pdq(:,:,igcm_ccnco2_mass) < 0.)
1919	pdq(:,:,igcm_ccnco2_mass) =
1920	& - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
1921	pdq(:,:,igcm_ccnco2_number) =
1922	& - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
1923	end where
1924	where (pq(:,:,igcm_ccnco2_number) +
1925	& ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.)
1926	pdq(:,:,igcm_ccnco2_mass) =
1927	& - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
1928	pdq(:,:,igcm_ccnco2_number) =
1929	& - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
1930	end where
1931
1932	c Negative values?
1933	where (pq(:,:,igcm_dust_mass) +
1934	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
1935	pdq(:,:,igcm_dust_mass) =
1936	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1937	pdq(:,:,igcm_dust_number) =
1938	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1939	end where
1940	where (pq(:,:,igcm_dust_number) +
1941	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
1942	pdq(:,:,igcm_dust_mass) =
1943	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1944	pdq(:,:,igcm_dust_number) =
1945	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1946	end where
1947	if (meteo_flux) then
1948	where (pq(:,:,igcm_ccnco2_meteor_mass) +
1949	& ptimestep*pdq(:,:,igcm_ccnco2_meteor_mass) < 0.)
1950	pdq(:,:,igcm_ccnco2_meteor_mass) =
1951	& - pq(:,:,igcm_ccnco2_meteor_mass)/ptimestep + 1.e-30
1952	pdq(:,:,igcm_ccnco2_meteor_number) =
1953	& - pq(:,:,igcm_ccnco2_meteor_number)/ptimestep + 1.e-30
1954	end where
1955	where (pq(:,:,igcm_ccnco2_meteor_number) +
1956	& ptimestep*pdq(:,:,igcm_ccnco2_meteor_number) < 0.)
1957	pdq(:,:,igcm_ccnco2_meteor_mass) =
1958	& - pq(:,:,igcm_ccnco2_meteor_mass)/ptimestep + 1.e-30
1959	pdq(:,:,igcm_ccnco2_meteor_number) =
1960	& - pq(:,:,igcm_ccnco2_meteor_number)/ptimestep + 1.e-30
1961	end where
1962	end if
1963	END IF ! of IF (co2clouds)
1964
1965	c 9b. Aerosol particles
1966	c -------------------
1967	c ----------
1968	c Dust devil :
1969	c ----------
1970	IF(callddevil) then
1971	call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
1972	& zdqdev,zdqsdev)
1973
1974	if (dustbin.ge.1) then
1975	do iq=1,nq
1976	DO l=1,nlayer
1977	DO ig=1,ngrid
1978	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
1979	ENDDO
1980	ENDDO
1981	enddo
1982	do iq=1,nq
1983	DO ig=1,ngrid
1984	DO islope = 1,nslope
1985	dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) +
1986	& zdqsdev(ig,iq)cos(pidef_slope_mean(islope)/180.)
1987	ENDDO
1988	ENDDO
1989	enddo
1990	endif ! of if (dustbin.ge.1)
1991
1992	END IF ! of IF (callddevil)
1993
1994	c -------------
1995	c Sedimentation : acts also on water ice
1996	c -------------
1997	IF (sedimentation) THEN
1998	zdqsed(1:ngrid,1:nlayer,1:nq)=0
1999	zdqssed(1:ngrid,1:nq)=0
2000
2001	c Sedimentation for co2 clouds tracers are inside co2cloud microtimestep
2002	c Zdqssed isn't
2003	call callsedim(ngrid,nlayer,ptimestep,
2004	& zplev,zzlev,zzlay,pt,pdt,
2005	& rdust,rstormdust,rtopdust,
2006	& rice,rsedcloud,rhocloud,
2007	& pq,pdq,zdqsed,zdqssed,nq,
2008	& tau,tauscaling)
2009
2010	c Flux at the surface of co2 ice computed in co2cloud microtimestep
2011	IF (rdstorm) THEN
2012	c Storm dust cannot sediment to the surface
2013	DO ig=1,ngrid
2014	zdqsed(ig,1,igcm_stormdust_mass)=
2015	& zdqsed(ig,1,igcm_stormdust_mass)+
2016	& zdqssed(ig,igcm_stormdust_mass) /
2017	& ((pplev(ig,1)-pplev(ig,2))/g)
2018	zdqsed(ig,1,igcm_stormdust_number)=
2019	& zdqsed(ig,1,igcm_stormdust_number)+
2020	& zdqssed(ig,igcm_stormdust_number) /
2021	& ((pplev(ig,1)-pplev(ig,2))/g)
2022	zdqssed(ig,igcm_stormdust_mass)=0.
2023	zdqssed(ig,igcm_stormdust_number)=0.
2024	ENDDO
2025	ENDIF !rdstorm
2026
2027	DO iq=1, nq
2028	DO l=1,nlayer
2029	DO ig=1,ngrid
2030	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
2031	ENDDO
2032	ENDDO
2033	ENDDO
2034	DO iq=1, nq
2035	DO ig=1,ngrid
2036	DO islope = 1,nslope
2037	dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) +
2038	& zdqssed(ig,iq)cos(pidef_slope_mean(islope)/180.)
2039	ENDDO
2040	ENDDO
2041	ENDDO
2042
2043	END IF ! of IF (sedimentation)
2044
2045	c Add lifted dust to tendancies after sedimentation in the LES (AC)
2046	IF (turb_resolved) THEN
2047	DO iq=1, nq
2048	DO l=1,nlayer
2049	DO ig=1,ngrid
2050	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
2051	ENDDO
2052	ENDDO
2053	ENDDO
2054	DO iq=1, nq
2055	DO ig=1,ngrid
2056	dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:)
2057	ENDDO
2058	ENDDO
2059	ENDIF
2060	c
2061	c 9c. Chemical species
2062	c ------------------
2063
2064	#ifndef MESOSCALE
2065	c --------------
2066	c photochemistry :
2067	c --------------
2068	IF (photochem) then
2069
2070	if (modulo(icount-1,ichemistry).eq.0) then
2071	! compute chemistry every ichemistry physics step
2072
2073	! dust and ice surface area
2074	call surfacearea(ngrid, nlayer, naerkind,
2075	$ ptimestep, zplay, zzlay,
2076	$ pt, pq, pdq, nq,
2077	$ rdust, rice, tau, tauscaling,
2078	$ surfdust, surfice)
2079	! call photochemistry
2080	DO ig = 1,ngrid
2081	qsurf_tmp(ig,:) = qsurf(ig,:,major_slope(ig))
2082	ENDDO
2083	call calchim(ngrid,nlayer,nq,
2084	& ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
2085	$ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
2086	$ zdqcloud,zdqscloud,tau(:,1),
2087	$ qsurf_tmp(:,igcm_co2),
2088	$ pu,pdu,pv,pdv,surfdust,surfice)
2089
2090	endif ! of if (modulo(icount-1,ichemistry).eq.0)
2091
2092	! increment values of tracers:
2093	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
2094	! tracers is zero anyways
2095	DO l=1,nlayer
2096	DO ig=1,ngrid
2097	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
2098	ENDDO
2099	ENDDO
2100	ENDDO ! of DO iq=1,nq
2101
2102	! add condensation tendency for H2O2
2103	if (igcm_h2o2.ne.0) then
2104	DO l=1,nlayer
2105	DO ig=1,ngrid
2106	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
2107	& +zdqcloud(ig,l,igcm_h2o2)
2108	ENDDO
2109	ENDDO
2110	endif
2111
2112	! increment surface values of tracers:
2113	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
2114	! tracers is zero anyways
2115	DO ig=1,ngrid
2116	DO islope = 1,nslope
2117	dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope) +
2118	& zdqschim(ig,iq)cos(pidef_slope_mean(islope)/180.)
2119	ENDDO
2120	ENDDO
2121	ENDDO ! of DO iq=1,nq
2122
2123	! add condensation tendency for H2O2
2124	if (igcm_h2o2.ne.0) then
2125	DO ig=1,ngrid
2126	DO islope = 1,nslope
2127	dqsurf(ig,igcm_h2o2,islope)=dqsurf(ig,igcm_h2o2,islope)+
2128	& zdqscloud(ig,igcm_h2o2)cos(pidef_slope_mean(islope)/180.)
2129	ENDDO
2130	ENDDO
2131	endif
2132
2133	END IF ! of IF (photochem)
2134	#endif
2135
2136
2137	#ifndef MESOSCALE
2138	c-----------------------------------------------------------------------
2139	c 10. THERMOSPHERE CALCULATION
2140	c-----------------------------------------------------------------------
2141
2142	if (callthermos) then
2143	call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol,
2144	$ mu0,ptimestep,ptime,zday,tsurf_meshavg,zzlev,zzlay,
2145	& pt,pq,pu,pv,pdt,pdq,
2146	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff,
2147	$ PhiEscH,PhiEscH2,PhiEscD)
2148
2149	DO l=1,nlayer
2150	DO ig=1,ngrid
2151	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
2152	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)+zdteuv(ig,l)
2153	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
2154	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
2155	DO iq=1, nq
2156	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
2157	ENDDO
2158	ENDDO
2159	ENDDO
2160
2161	endif ! of if (callthermos)
2162	#endif
2163	c-----------------------------------------------------------------------
2164	c 11. Carbon dioxide condensation-sublimation:
2165	c (should be the last atmospherical physical process to be computed)
2166	c -------------------------------------------
2167	IF (tituscap) THEN
2168	!!! get the actual co2 seasonal cap from Titus observations
2169	CALL geticecover(ngrid, 180.*zls/pi,
2170	. 180.longitude/pi, 180.latitude/pi,
2171	. qsurf_tmp(:,igcm_co2) )
2172	qsurf_tmp(:,igcm_co2) = qsurf_tmp(:,igcm_co2) * 10000.
2173	ENDIF
2174
2175
2176	IF (callcond) THEN
2177	zdtc(:,:) = 0.
2178	zdtsurfc(:,:) = 0.
2179	zduc(:,:) = 0.
2180	zdvc(:,:) = 0.
2181	zdqc(:,:,:) = 0.
2182	zdqsc(:,:,:) = 0.
2183	CALL co2condens(ngrid,nlayer,nq,ptimestep,
2184	$ capcal,zplay,zplev,tsurf,pt,
2185	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
2186	$ qsurf(:,igcm_co2,:),albedo,emis,rdust,
2187	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
2188	$ fluxsurf_dn_sw,zls,
2189	$ zdqssed_co2,zcondicea_co2microp,
2190	& zdqsc)
2191
2192	DO iq=1, nq
2193	DO ig=1,ngrid
2194	dqsurf(ig,iq,:)=dqsurf(ig,iq,:)+zdqsc(ig,iq,:)
2195	ENDDO ! (ig)
2196	ENDDO ! (iq)
2197	DO l=1,nlayer
2198	DO ig=1,ngrid
2199	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
2200	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
2201	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
2202	ENDDO
2203	ENDDO
2204	DO ig=1,ngrid
2205	zdtsurf(ig,:) = zdtsurf(ig,:) + zdtsurfc(ig,:)
2206	ENDDO
2207
2208	DO iq=1, nq
2209	DO l=1,nlayer
2210	DO ig=1,ngrid
2211	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
2212	ENDDO
2213	ENDDO
2214	ENDDO
2215
2216	#ifndef MESOSCALE
2217	! update surface pressure
2218	DO ig=1,ngrid
2219	ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
2220	ENDDO
2221	! update pressure levels
2222	DO l=1,nlayer
2223	DO ig=1,ngrid
2224	zplay(ig,l) = aps(l) + bps(l)*ps(ig)
2225	zplev(ig,l) = ap(l) + bp(l)*ps(ig)
2226	ENDDO
2227	ENDDO
2228	zplev(:,nlayer+1) = 0.
2229	! update layers altitude
2230	DO l=2,nlayer
2231	DO ig=1,ngrid
2232	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
2233	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
2234	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
2235	ENDDO
2236	ENDDO
2237	#endif
2238	ENDIF ! of IF (callcond)
2239
2240	c-----------------------------------------------------------------------
2241	c Updating tracer budget on surface
2242	c-----------------------------------------------------------------------
2243	DO iq=1, nq
2244	DO ig=1,ngrid
2245	DO islope = 1,nslope
2246	qsurf(ig,iq,islope)=qsurf(ig,iq,islope)+
2247	& ptimestep*dqsurf(ig,iq,islope)
2248	ENDDO
2249	ENDDO ! (ig)
2250	ENDDO ! (iq)
2251	c-----------------------------------------------------------------------
2252	c 12. Surface and sub-surface soil temperature
2253	c-----------------------------------------------------------------------
2254	c
2255	c
2256	c 12.1 Increment Surface temperature:
2257	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2258
2259	DO ig=1,ngrid
2260	DO islope = 1,nslope
2261	tsurf(ig,islope)=tsurf(ig,islope)+
2262	& ptimestep*zdtsurf(ig,islope)
2263	ENDDO
2264	ENDDO
2265
2266	c Prescribe a cold trap at south pole (except at high obliquity !!)
2267	c Temperature at the surface is set there to be the temperature
2268	c corresponding to equilibrium temperature between phases of CO2
2269
2270
2271	IF (water) THEN
2272	!#ifndef MESOSCALE
2273	! if (caps.and.(obliquit.lt.27.)) then => now done in co2condens
2274	! NB: Updated surface pressure, at grid point 'ngrid', is
2275	! ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
2276	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
2277	! & (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
2278	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095ps(ngrid)))
2279	! endif
2280	!#endif
2281	c -------------------------------------------------------------
2282	c Change of surface albedo in case of ground frost
2283	c everywhere except on the north permanent cap and in regions
2284	c covered by dry ice.
2285	c ALWAYS PLACE these lines after co2condens !!!
2286	c -------------------------------------------------------------
2287	do ig=1,ngrid
2288	do islope = 1,nslope
2289	if ((qsurf(ig,igcm_co2,islope).eq.0).and.
2290	& (qsurf(ig,igcm_h2o_ice,islope)
2291	& .gt.frost_albedo_threshold)) then
2292	if ((watercaptag(ig)).and.(cst_cap_albedo)) then
2293	albedo(ig,1,islope) = albedo_h2o_cap
2294	albedo(ig,2,islope) = albedo_h2o_cap
2295	else
2296	albedo(ig,1,islope) = albedo_h2o_frost
2297	albedo(ig,2,islope) = albedo_h2o_frost
2298	endif !((watercaptag(ig)).and.(cst_cap_albedo)) then
2299	c write(,) "frost thickness", qsurf(ig,igcm_h2o_ice)
2300	c write(,) "physiq.F frost :"
2301	c & ,latitude(ig)180./pi, longitude(ig)180./pi
2302	endif
2303	enddo ! islope
2304	enddo ! of do ig=1,ngrid
2305	ENDIF ! of IF (water)
2306
2307	c
2308	c 12.2 Compute soil temperatures and subsurface heat flux:
2309	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2310	IF (callsoil) THEN
2311	c Thermal inertia feedback
2312	IF (tifeedback) THEN
2313	DO islope = 1,nslope
2314	CALL soil_tifeedback(ngrid,nsoilmx,
2315	s qsurf(:,:,islope),inertiesoil(:,:,islope))
2316	ENDDO
2317	CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
2318	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
2319	ELSE
2320	CALL soil(ngrid,nsoilmx,.false.,inertiedat_tmp,
2321	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
2322	ENDIF
2323	ENDIF
2324
2325	c To avoid negative values
2326	IF (rdstorm) THEN
2327	where (pq(:,:,igcm_stormdust_mass) +
2328	& ptimestep*pdq(:,:,igcm_stormdust_mass) < 0.)
2329	pdq(:,:,igcm_stormdust_mass) =
2330	& - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30
2331	pdq(:,:,igcm_stormdust_number) =
2332	& - pq(:,:,igcm_stormdust_number)/ptimestep + 1.e-30
2333	end where
2334	where (pq(:,:,igcm_stormdust_number) +
2335	& ptimestep*pdq(:,:,igcm_stormdust_number) < 0.)
2336	pdq(:,:,igcm_stormdust_mass) =
2337	& - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30
2338	pdq(:,:,igcm_stormdust_number) =
2339	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
2340	end where
2341
2342	where (pq(:,:,igcm_dust_mass) +
2343	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
2344	pdq(:,:,igcm_dust_mass) =
2345	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
2346	pdq(:,:,igcm_dust_number) =
2347	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
2348	end where
2349	where (pq(:,:,igcm_dust_number) +
2350	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
2351	pdq(:,:,igcm_dust_mass) =
2352	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
2353	pdq(:,:,igcm_dust_number) =
2354	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
2355	end where
2356	ENDIF !(rdstorm)
2357
2358	c-----------------------------------------------------------------------
2359	c J. Naar : Surface and sub-surface water ice
2360	c-----------------------------------------------------------------------
2361	c
2362	c
2363	c Increment Watercap (surface h2o reservoirs):
2364	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2365
2366	DO ig=1,ngrid
2367	DO islope = 1,nslope
2368	watercap(ig,islope)=watercap(ig,islope)+
2369	s ptimestep*dwatercap(ig,islope)
2370	ENDDO
2371	ENDDO
2372
2373	IF (refill_watercap) THEN
2374
2375	DO ig=1,ngrid
2376	DO islope = 1,nslope
2377	if (watercaptag(ig).and. (qsurf(ig,igcm_h2o_ice,islope)
2378	& .gt.frost_metam_threshold)) then
2379
2380	watercap(ig,islope)=watercap(ig,islope)
2381	& +qsurf(ig,igcm_h2o_ice,islope)
2382	& - frost_metam_threshold
2383	qsurf(ig,igcm_h2o_ice,islope) = frost_metam_threshold
2384	endif ! (watercaptag(ig).and.
2385	ENDDO
2386	ENDDO
2387
2388	ENDIF ! (refill_watercap) THEN
2389
2390
2391	c-----------------------------------------------------------------------
2392	c 13. Write output files
2393	c ----------------------
2394
2395	call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf,
2396	& albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg)
2397
2398	c -------------------------------
2399	c Dynamical fields incrementation
2400	c -------------------------------
2401	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
2402	! temperature, zonal and meridional wind
2403	DO l=1,nlayer
2404	DO ig=1,ngrid
2405	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
2406	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
2407	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
2408	ENDDO
2409	ENDDO
2410
2411	! tracers
2412	DO iq=1, nq
2413	DO l=1,nlayer
2414	DO ig=1,ngrid
2415	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
2416	ENDDO
2417	ENDDO
2418	ENDDO
2419
2420	! Density
2421	DO l=1,nlayer
2422	DO ig=1,ngrid
2423	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
2424	ENDDO
2425	ENDDO
2426
2427	! Potential Temperature
2428
2429	DO ig=1,ngrid
2430	DO l=1,nlayer
2431	zh(ig,l) = zt(ig,l)(zplev(ig,1)/zplay(ig,l))*rcp
2432	ENDDO
2433	ENDDO
2434
2435	c Compute surface stress : (NB: z0 is a common in surfdat.h)
2436	c DO ig=1,ngrid
2437	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
2438	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
2439	c ENDDO
2440
2441	c Sum of fluxes in solar spectral bands (for output only)
2442	fluxtop_dn_sw_tot(1:ngrid)=fluxtop_dn_sw(1:ngrid,1) +
2443	& fluxtop_dn_sw(1:ngrid,2)
2444	fluxtop_up_sw_tot(1:ngrid)=fluxtop_up_sw(1:ngrid,1) +
2445	& fluxtop_up_sw(1:ngrid,2)
2446	fluxsurf_dn_sw_tot(1:ngrid,1:nslope)=
2447	& fluxsurf_dn_sw(1:ngrid,1,1:nslope) +
2448	& fluxsurf_dn_sw(1:ngrid,2,1:nslope)
2449	fluxsurf_up_sw_tot(1:ngrid)=fluxsurf_up_sw(1:ngrid,1) +
2450	& fluxsurf_up_sw(1:ngrid,2)
2451
2452	c ***** TEST ****************************************************
2453	ztim1 = 999
2454	DO l=1,nlayer
2455	DO ig=1,ngrid
2456	if (pt(ig,l).lt.ztim1) then
2457	ztim1 = pt(ig,l)
2458	igmin = ig
2459	lmin = l
2460	end if
2461	ENDDO
2462	ENDDO
2463	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
2464	write(,) 'PHYSIQ: stability WARNING :'
2465	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
2466	& 'ig l =', igmin, lmin
2467	end if
2468	c *******************************************************************
2469
2470	c ---------------------
2471	c Outputs to the screen
2472	c ---------------------
2473
2474	IF (lwrite) THEN
2475	PRINT*,'Global diagnostics for the physics'
2476	PRINT,'Variables and their increments x and dx/dt dt'
2477	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
2478	WRITE(*,'(2f10.5,2f15.5)')
2479	s tsurf(igout,:),zdtsurf(igout,:)*ptimestep,
2480	s zplev(igout,1),pdpsrf(igout)*ptimestep
2481	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
2482	WRITE(*,'(i4,6f10.5)') (l,
2483	s pu(igout,l),pdu(igout,l)*ptimestep,
2484	s pv(igout,l),pdv(igout,l)*ptimestep,
2485	s pt(igout,l),pdt(igout,l)*ptimestep,
2486	s l=1,nlayer)
2487	ENDIF ! of IF (lwrite)
2488
2489	c ----------------------------------------------------------
2490	c ----------------------------------------------------------
2491	c INTERPOLATIONS IN THE SURFACE-LAYER
2492	c ----------------------------------------------------------
2493	c ----------------------------------------------------------
2494
2495	n_out=0 ! number of elements in the z_out array.
2496	! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set
2497	! to 5
2498	IF (n_out .ne. 0) THEN
2499
2500	IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out))
2501	IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out))
2502	IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out))
2503
2504	z_out(:)=[3.,2.,1.,0.5,0.1]
2505	u_out(:,:)=0.
2506	T_out(:,:)=0.
2507
2508	call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
2509	& g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf(:,iflat),
2510	& zh,z_out,n_out,T_out,u_out,ustar,tstar,L_mo,vhf,vvv)
2511	! pourquoi ustar recalcule ici? fait dans vdifc.
2512
2513	#ifndef MESOSCALE
2514	DO n=1,n_out
2515	write(zstring, '(F8.6)') z_out(n)
2516	call write_output('T_out_'//trim(zstring),
2517	& 'potential temperature at z_out','K',T_out(:,n))
2518	call write_output('u_out_'//trim(zstring),
2519	& 'horizontal velocity norm at z_out','m/s',u_out(:,n))
2520	ENDDO
2521	call write_output('u_star',
2522	& 'friction velocity','m/s',ustar)
2523	call write_output('teta_star',
2524	& 'friction potential temperature','K',tstar)
2525	call write_output('vvv',
2526	& 'Vertical velocity variance at zout','m',vvv)
2527	call write_output('vhf',
2528	& 'Vertical heat flux at zout','m',vhf)
2529	#else
2530	T_out1(:)=T_out(:,1)
2531	u_out1(:)=u_out(:,1)
2532	#endif
2533
2534	ENDIF
2535
2536	c ----------------------------------------------------------
2537	c ----------------------------------------------------------
2538	c END OF SURFACE LAYER INTERPOLATIONS
2539	c ----------------------------------------------------------
2540	c ----------------------------------------------------------
2541
2542	#ifndef MESOSCALE
2543	c -------------------------------------------------------------------
2544	c Writing NetCDF file "RESTARTFI" at the end of the run
2545	c -------------------------------------------------------------------
2546	c Note: 'restartfi' is stored just before dynamics are stored
2547	c in 'restart'. Between now and the writting of 'restart',
2548	c there will have been the itau=itau+1 instruction and
2549	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
2550	c thus we store for time=time+dtvr
2551
2552	! default: not writing a restart file at this time step
2553	write_restart=.false.
2554	IF (ecritstart.GT.0) THEN
2555	! For when we store multiple time steps in the restart file
2556	IF (MODULO(icount*iphysiq,ecritstart).EQ.0) THEN
2557	write_restart=.true.
2558	ENDIF
2559	ENDIF
2560	IF (lastcall) THEN
2561	! Always write a restart at the end of the simulation
2562	write_restart=.true.
2563	ENDIF
2564
2565	IF (write_restart) THEN
2566	IF (grid_type==unstructured) THEN !IF DYNAMICO
2567
2568	! When running Dynamico, no need to add a dynamics time step to ztime_fin
2569	IF (ptime.LE. 1.E-10) THEN
2570	! Residual ptime occurs with Dynamico
2571	ztime_fin = pday !+ ptime + ptimestep/(float(iphysiq)*daysec)
2572	. - day_ini - time_phys
2573	ELSE
2574	ztime_fin = pday + ptime !+ ptimestep/(float(iphysiq)*daysec)
2575	. - day_ini - time_phys
2576	ENDIF
2577	if (ecritstart==0) then
2578	ztime_fin = ztime_fin-(day_end-day_ini)
2579	endif
2580
2581	ELSE ! IF LMDZ
2582
2583	if (ecritstart.GT.0) then !IF MULTIPLE RESTARTS nothing change
2584	ztime_fin = pday - day_ini + ptime
2585	. + ptimestep/(float(iphysiq)*daysec)
2586	else !IF ONE RESTART final time in top of day_end
2587	ztime_fin = pday - day_ini-(day_end-day_ini)
2588	. + ptime + ptimestep/(float(iphysiq)*daysec)
2589	endif
2590
2591	ENDIF ! of IF (grid_type==unstructured)
2592	write(*,'(A,I7,A,F12.5)')
2593	. 'PHYSIQ: writing in restartfi ; icount=',
2594	. icount,' date=',ztime_fin
2595
2596	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
2597	. ptimestep,ztime_fin,
2598	. tsurf,tsoil,albedo,
2599	. emis,q2,qsurf,tauscaling,totcloudfrac,wstar,
2600	. watercap)
2601	ENDIF ! of IF (write_restart)
2602
2603	#endif
2604
2605	c IF (ngrid.NE.1) then
2606
2607	c -------------------------------------------------------------------
2608	c Calculation of diagnostic variables written in both stats and
2609	c diagfi files
2610	c -------------------------------------------------------------------
2611
2612
2613	do ig=1,ngrid
2614	if(mu0(ig).le.0.01) then
2615	fluxsurf_dir_dn_sw(ig) = 0.
2616	else
2617	fluxsurf_dir_dn_sw(ig) = flux_1AU/dist_sol/dist_solmu0(ig)
2618	& exp(-(tau(ig,iaer_dust_doubleq)+
2619	& tau(ig,iaer_h2o_ice))/mu0(ig))
2620	endif
2621	enddo
2622
2623	! Density-scaled opacities
2624	do ig=1,ngrid
2625	dsodust(ig,:) =
2626	& dsodust(ig,:)*tauscaling(ig)
2627	dsords(ig,:) =
2628	& dsords(ig,:)*tauscaling(ig)
2629	dsotop(ig,:) =
2630	& dsotop(ig,:)*tauscaling(ig)
2631	enddo
2632
2633	if(doubleq) then
2634	do ig=1,ngrid
2635	IF (sedimentation) THEN
2636	dqdustsurf(ig) =
2637	& zdqssed(ig,igcm_dust_mass)*tauscaling(ig)
2638	dndustsurf(ig) =
2639	& zdqssed(ig,igcm_dust_number)*tauscaling(ig)
2640	ENDIF
2641	ndust(ig,:) =
2642	& zq(ig,:,igcm_dust_number)*tauscaling(ig)
2643	qdust(ig,:) =
2644	& zq(ig,:,igcm_dust_mass)*tauscaling(ig)
2645	enddo
2646	if (scavenging) then
2647	do ig=1,ngrid
2648	IF (sedimentation) THEN
2649	dqdustsurf(ig) = dqdustsurf(ig) +
2650	& zdqssed(ig,igcm_ccn_mass)*tauscaling(ig)
2651	dndustsurf(ig) = dndustsurf(ig) +
2652	& zdqssed(ig,igcm_ccn_number)*tauscaling(ig)
2653	ENDIF
2654	nccn(ig,:) =
2655	& zq(ig,:,igcm_ccn_number)*tauscaling(ig)
2656	qccn(ig,:) =
2657	& zq(ig,:,igcm_ccn_mass)*tauscaling(ig)
2658	enddo
2659	endif
2660	endif ! of (doubleq)
2661
2662	if (rdstorm) then ! diagnostics of stormdust tendancies for 1D and 3D
2663	mstormdtot(:)=0
2664	mdusttot(:)=0
2665	qdusttotal(:,:)=0
2666	do ig=1,ngrid
2667	rdsdqdustsurf(ig) =
2668	& zdqssed(ig,igcm_stormdust_mass)*tauscaling(ig)
2669	rdsdndustsurf(ig) =
2670	& zdqssed(ig,igcm_stormdust_number)*tauscaling(ig)
2671	rdsndust(ig,:) =
2672	& pq(ig,:,igcm_stormdust_number)*tauscaling(ig)
2673	rdsqdust(ig,:) =
2674	& pq(ig,:,igcm_stormdust_mass)*tauscaling(ig)
2675	do l=1,nlayer
2676	mstormdtot(ig) = mstormdtot(ig) +
2677	& zq(ig,l,igcm_stormdust_mass) *
2678	& (zplev(ig,l) - zplev(ig,l+1)) / g
2679	mdusttot(ig) = mdusttot(ig) +
2680	& zq(ig,l,igcm_dust_mass) *
2681	& (zplev(ig,l) - zplev(ig,l+1)) / g
2682	qdusttotal(ig,l) = qdust(ig,l)+rdsqdust(ig,l) !calculate total dust
2683	enddo
2684	enddo
2685	endif !(rdstorm)
2686
2687	if (water) then
2688	mtot(:)=0
2689	icetot(:)=0
2690	rave(:)=0
2691	tauTES(:)=0
2692
2693	IF (hdo) then
2694	mtotD(:)=0
2695	icetotD(:)=0
2696	ENDIF !hdo
2697
2698	do ig=1,ngrid
2699	do l=1,nlayer
2700	mtot(ig) = mtot(ig) +
2701	& zq(ig,l,igcm_h2o_vap) *
2702	& (zplev(ig,l) - zplev(ig,l+1)) / g
2703	icetot(ig) = icetot(ig) +
2704	& zq(ig,l,igcm_h2o_ice) *
2705	& (zplev(ig,l) - zplev(ig,l+1)) / g
2706	IF (hdo) then
2707	mtotD(ig) = mtotD(ig) +
2708	& zq(ig,l,igcm_hdo_vap) *
2709	& (zplev(ig,l) - zplev(ig,l+1)) / g
2710	icetotD(ig) = icetotD(ig) +
2711	& zq(ig,l,igcm_hdo_ice) *
2712	& (zplev(ig,l) - zplev(ig,l+1)) / g
2713	ENDIF !hdo
2714
2715	c Computing abs optical depth at 825 cm-1 in each
2716	c layer to simulate NEW TES retrieval
2717	Qabsice = min(
2718	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
2719	& )
2720	opTES(ig,l)= 0.75 * Qabsice *
2721	& zq(ig,l,igcm_h2o_ice) *
2722	& (zplev(ig,l) - zplev(ig,l+1)) / g
2723	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
2724	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
2725	enddo
2726	c rave(ig)=rave(ig)/max(icetot(ig),1.e-30) ! mass weight
2727	c if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
2728	enddo
2729	call watersat(ngrid*nlayer,zt,zplay,zqsat)
2730	satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
2731
2732	if (scavenging) then
2733	Nccntot(:)= 0
2734	Mccntot(:)= 0
2735	rave(:)=0
2736	do ig=1,ngrid
2737	do l=1,nlayer
2738	Nccntot(ig) = Nccntot(ig) +
2739	& zq(ig,l,igcm_ccn_number)*tauscaling(ig)
2740	& *(zplev(ig,l) - zplev(ig,l+1)) / g
2741	Mccntot(ig) = Mccntot(ig) +
2742	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
2743	& *(zplev(ig,l) - zplev(ig,l+1)) / g
2744	cccc Column integrated effective ice radius
2745	cccc is weighted by total ice surface area (BETTER than total ice mass)
2746	rave(ig) = rave(ig) +
2747	& tauscaling(ig) *
2748	& zq(ig,l,igcm_ccn_number) *
2749	& (zplev(ig,l) - zplev(ig,l+1)) / g *
2750	& rice(ig,l) * rice(ig,l)* (1.+nuice_ref)
2751	enddo
2752	rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75
2753	& /max(pi*rave(ig),1.e-30) ! surface weight
2754	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
2755	enddo
2756	else ! of if (scavenging)
2757	rave(:)=0
2758	do ig=1,ngrid
2759	do l=1,nlayer
2760	rave(ig) = rave(ig) +
2761	& zq(ig,l,igcm_h2o_ice) *
2762	& (zplev(ig,l) - zplev(ig,l+1)) / g *
2763	& rice(ig,l) * (1.+nuice_ref)
2764	enddo
2765	rave(ig) = max(rave(ig) /
2766	& max(icetot(ig),1.e-30),1.e-30) ! mass weight
2767	enddo
2768	endif ! of if (scavenging)
2769
2770	!Alternative A. Pottier weighting
2771	rave2(:) = 0.
2772	totrave2(:) = 0.
2773	do ig=1,ngrid
2774	do l=1,nlayer
2775	rave2(ig) =rave2(ig)+ zq(ig,l,igcm_h2o_ice)*rice(ig,l)
2776	totrave2(ig) = totrave2(ig) + zq(ig,l,igcm_h2o_ice)
2777	end do
2778	rave2(ig)=max(rave2(ig)/max(totrave2(ig),1.e-30),1.e-30)
2779	end do
2780
2781	endif ! of if (water)
2782
2783	if (co2clouds) then
2784	mtotco2(1:ngrid) = 0.
2785	icetotco2(1:ngrid) = 0.
2786	vaptotco2(1:ngrid) = 0.
2787	do ig=1,ngrid
2788	do l=1,nlayer
2789	vaptotco2(ig) = vaptotco2(ig) +
2790	& zq(ig,l,igcm_co2) *
2791	& (zplev(ig,l) - zplev(ig,l+1)) / g
2792	icetotco2(ig) = icetot(ig) +
2793	& zq(ig,l,igcm_co2_ice) *
2794	& (zplev(ig,l) - zplev(ig,l+1)) / g
2795	end do
2796	mtotco2(ig) = icetotco2(ig) + vaptotco2(ig)
2797	end do
2798	end if
2799
2800	#ifndef MESOSCALE
2801	c -----------------------------------------------------------------
2802	c WSTATS: Saving statistics
2803	c -----------------------------------------------------------------
2804	c ("stats" stores and accumulates key variables in file "stats.nc"
2805	c which can later be used to make the statistic files of the run:
2806	c if flag "callstats" from callphys.def is .true.)
2807
2808	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
2809	call wstats(ngrid,"tsurf","Surface temperature","K",2
2810	& ,tsurf(:,iflat))
2811	call wstats(ngrid,"co2ice","CO2 ice cover",
2812	& "kg.m-2",2,qsurf(:,igcm_co2,iflat))
2813	call wstats(ngrid,"watercap","H2O ice cover",
2814	& "kg.m-2",2,watercap(:,iflat))
2815	call wstats(ngrid,"tau_pref_scenario",
2816	& "prescribed visible dod at 610 Pa","NU",
2817	& 2,tau_pref_scenario)
2818	call wstats(ngrid,"tau_pref_gcm",
2819	& "visible dod at 610 Pa in the GCM","NU",
2820	& 2,tau_pref_gcm)
2821	call wstats(ngrid,"fluxsurf_lw",
2822	& "Thermal IR radiative flux to surface","W.m-2",2,
2823	& fluxsurf_lw(:,iflat))
2824	call wstats(ngrid,"fluxsurf_dn_sw",
2825	& "Incoming Solar radiative flux to surface","W.m-2",2,
2826	& fluxsurf_dn_sw_tot(:,iflat))
2827	call wstats(ngrid,"fluxsurf_up_sw",
2828	& "Reflected Solar radiative flux from surface","W.m-2",2,
2829	& fluxsurf_up_sw_tot)
2830	call wstats(ngrid,"fluxtop_lw",
2831	& "Thermal IR radiative flux to space","W.m-2",2,
2832	& fluxtop_lw)
2833	call wstats(ngrid,"fluxtop_dn_sw",
2834	& "Incoming Solar radiative flux from space","W.m-2",2,
2835	& fluxtop_dn_sw_tot)
2836	call wstats(ngrid,"fluxtop_up_sw",
2837	& "Outgoing Solar radiative flux to space","W.m-2",2,
2838	& fluxtop_up_sw_tot)
2839	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
2840	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
2841	call wstats(ngrid,"v","Meridional (North-South) wind",
2842	& "m.s-1",3,zv)
2843	call wstats(ngrid,"w","Vertical (down-up) wind",
2844	& "m.s-1",3,pw)
2845	call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
2846	call wstats(ngrid,"pressure","Pressure","Pa",3,zplay)
2847	call wstats(ngrid,"q2",
2848	& "Boundary layer eddy kinetic energy",
2849	& "m2.s-2",3,q2)
2850	call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
2851	& emis(:,iflat))
2852	call wstats(ngrid,"fluxsurf_dir_dn_sw",
2853	& "Direct incoming SW flux at surface",
2854	& "W.m-2",2,fluxsurf_dir_dn_sw)
2855
2856	if (calltherm) then
2857	call wstats(ngrid,"zmax_th","Height of thermals",
2858	& "m",2,zmax_th)
2859	call wstats(ngrid,"hfmax_th","Max thermals heat flux",
2860	& "K.m/s",2,hfmax_th)
2861	call wstats(ngrid,"wstar",
2862	& "Max vertical velocity in thermals",
2863	& "m/s",2,wstar)
2864	endif
2865
2866	if (water) then
2867	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
2868	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2869	call wstats(ngrid,"vmr_h2ovap",
2870	& "H2O vapor volume mixing ratio","mol/mol",
2871	& 3,vmr)
2872	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
2873	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
2874	call wstats(ngrid,"vmr_h2oice",
2875	& "H2O ice volume mixing ratio","mol/mol",
2876	& 3,vmr)
2877	! also store vmr_ice*rice for better diagnostics of rice
2878	vmr(1:ngrid,1:nlayer)=vmr(1:ngrid,1:nlayer)*
2879	& rice(1:ngrid,1:nlayer)
2880	call wstats(ngrid,"vmr_h2oice_rice",
2881	& "H2O ice mixing ratio times ice particule size",
2882	& "(mol/mol)*m",
2883	& 3,vmr)
2884	vmr=zqsat(1:ngrid,1:nlayer)
2885	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2886	call wstats(ngrid,"vmr_h2osat",
2887	& "saturation volume mixing ratio","mol/mol",
2888	& 3,vmr)
2889	call wstats(ngrid,"h2o_ice_s",
2890	& "surface h2o_ice","kg/m2",
2891	& 2,qsurf(1,igcm_h2o_ice,iflat))
2892	call wstats(ngrid,'albedo',
2893	& 'albedo',
2894	& '',2,albedo(1,1,iflat))
2895	call wstats(ngrid,"mtot",
2896	& "total mass of water vapor","kg/m2",
2897	& 2,mtot)
2898	call wstats(ngrid,"icetot",
2899	& "total mass of water ice","kg/m2",
2900	& 2,icetot)
2901	call wstats(ngrid,"reffice",
2902	& "Mean reff","m",
2903	& 2,rave)
2904	call wstats(ngrid,"Nccntot",
2905	& "condensation nuclei","Nbr/m2",
2906	& 2,Nccntot)
2907	call wstats(ngrid,"Mccntot",
2908	& "condensation nuclei mass","kg/m2",
2909	& 2,Mccntot)
2910	call wstats(ngrid,"rice",
2911	& "Ice particle size","m",
2912	& 3,rice)
2913	if (.not.activice) then
2914	call wstats(ngrid,"tauTESap",
2915	& "tau abs 825 cm-1","",
2916	& 2,tauTES)
2917	else
2918	call wstats(ngrid,'tauTES',
2919	& 'tau abs 825 cm-1',
2920	& '',2,taucloudtes)
2921	endif
2922
2923	endif ! of if (water)
2924
2925	if (co2clouds) then
2926	call wstats(ngrid,"mtotco2",
2927	& "total mass atm of co2","kg/m2",
2928	& 2,mtotco2)
2929	call wstats(ngrid,"icetotco2",
2930	& "total mass atm of co2 ice","kg/m2",
2931	& 2,icetotco2)
2932	call wstats(ngrid,"vaptotco2",
2933	& "total mass atm of co2 vapor","kg/m2",
2934	& 2,icetotco2)
2935	end if
2936
2937
2938	if (dustbin.ne.0) then
2939
2940	call wstats(ngrid,'tau','taudust','SI',2,tau(1,1))
2941
2942	if (doubleq) then
2943	call wstats(ngrid,'dqsdust',
2944	& 'deposited surface dust mass',
2945	& 'kg.m-2.s-1',2,dqdustsurf)
2946	call wstats(ngrid,'dqndust',
2947	& 'deposited surface dust number',
2948	& 'number.m-2.s-1',2,dndustsurf)
2949	call wstats(ngrid,'reffdust','reffdust',
2950	& 'm',3,rdust*ref_r0)
2951	call wstats(ngrid,'dustq','Dust mass mr',
2952	& 'kg/kg',3,qdust)
2953	call wstats(ngrid,'dustN','Dust number',
2954	& 'part/kg',3,ndust)
2955	if (rdstorm) then
2956	call wstats(ngrid,'reffstormdust','reffdust',
2957	& 'm',3,rstormdust*ref_r0)
2958	call wstats(ngrid,'rdsdustq','Dust mass mr',
2959	& 'kg/kg',3,rdsqdust)
2960	call wstats(ngrid,'rdsdustN','Dust number',
2961	& 'part/kg',3,rdsndust)
2962	end if
2963	else
2964	do iq=1,dustbin
2965	write(str2(1:2),'(i2.2)') iq
2966	call wstats(ngrid,'q'//str2,'mix. ratio',
2967	& 'kg/kg',3,zq(1,1,iq))
2968	call wstats(ngrid,'qsurf'//str2,'qsurf',
2969	& 'kg.m-2',2,qsurf(1,iq,iflat))
2970	end do
2971	endif ! (doubleq)
2972
2973	if (scavenging) then
2974	call wstats(ngrid,'ccnq','CCN mass mr',
2975	& 'kg/kg',3,qccn)
2976	call wstats(ngrid,'ccnN','CCN number',
2977	& 'part/kg',3,nccn)
2978	endif ! (scavenging)
2979
2980	endif ! (dustbin.ne.0)
2981
2982	if (photochem) then
2983	do iq=1,nq
2984	if (noms(iq) .ne. "dust_mass" .and.
2985	$ noms(iq) .ne. "dust_number" .and.
2986	$ noms(iq) .ne. "ccn_mass" .and.
2987	$ noms(iq) .ne. "ccn_number" .and.
2988	$ noms(iq) .ne. "ccnco2_mass" .and.
2989	$ noms(iq) .ne. "ccnco2_number" .and.
2990	$ noms(iq) .ne. "stormdust_mass" .and.
2991	$ noms(iq) .ne. "stormdust_number" .and.
2992	$ noms(iq) .ne. "topdust_mass" .and.
2993	$ noms(iq) .ne. "topdust_number") then
2994	! volume mixing ratio
2995
2996	vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
2997	& *mmean(1:ngrid,1:nlayer)/mmol(iq)
2998
2999	call wstats(ngrid,"vmr_"//trim(noms(iq)),
3000	$ "Volume mixing ratio","mol/mol",3,vmr)
3001	if ((noms(iq).eq."o")
3002	$ .or. (noms(iq).eq."co2")
3003	$ .or. (noms(iq).eq."o3")
3004	$ .or. (noms(iq).eq."ar")
3005	$ .or. (noms(iq).eq."o2")
3006	$ .or. (noms(iq).eq."h2o_vap") ) then
3007	call write_output("vmr_"//trim(noms(iq)),
3008	$ "Volume mixing ratio","mol/mol",vmr(:,:))
3009	end if
3010
3011	! number density (molecule.cm-3)
3012
3013	rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
3014	& rho(1:ngrid,1:nlayer)n_avog/
3015	& (1000*mmol(iq))
3016
3017	call wstats(ngrid,"num_"//trim(noms(iq)),
3018	$ "Number density","cm-3",3,rhopart)
3019	call write_output("num_"//trim(noms(iq)),
3020	$ "Number density","cm-3",rhopart(:,:))
3021
3022	! vertical column (molecule.cm-2)
3023
3024	do ig = 1,ngrid
3025	colden(ig,iq) = 0.
3026	end do
3027	do l=1,nlayer
3028	do ig=1,ngrid
3029	colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
3030	$ *(zplev(ig,l)-zplev(ig,l+1))
3031	$ 6.022e22/(mmol(iq)g)
3032	end do
3033	end do
3034
3035	call wstats(ngrid,"c_"//trim(noms(iq)),
3036	$ "column","mol cm-2",2,colden(1,iq))
3037	call write_output("c_"//trim(noms(iq)),
3038	$ "column","mol cm-2",colden(:,iq))
3039
3040	! global mass (g)
3041
3042	call planetwide_sumval(colden(:,iq)/6.022e23
3043	$ mmol(iq)1.e4*cell_area(:),mass(iq))
3044
3045	call write_output("mass_"//trim(noms(iq)),
3046	$ "global mass","g",mass(iq))
3047
3048	end if ! of if (noms(iq) .ne. "dust_mass" ...)
3049	end do ! of do iq=1,nq
3050	end if ! of if (photochem)
3051
3052
3053	IF(lastcall.and.callstats) THEN
3054	write (,) "Writing stats..."
3055	call mkstats(ierr)
3056	ENDIF
3057
3058
3059	c (Store EOF for Mars Climate database software)
3060	IF (calleofdump) THEN
3061	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
3062	ENDIF
3063	#endif
3064	!endif of ifndef MESOSCALE
3065
3066	#ifdef MESOSCALE
3067
3068	!! see comm_wrf.
3069	!! not needed when an array is already in a shared module.
3070	!! --> example : hfmax_th, zmax_th
3071
3072	CALL allocate_comm_wrf(ngrid,nlayer)
3073
3074	!state real HR_SW ikj misc 1 - h "HR_SW" "HEATING RATE SW" "K/s"
3075	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
3076	!state real HR_LW ikj misc 1 - h "HR_LW" "HEATING RATE LW" "K/s"
3077	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
3078	!state real SWDOWNZ ij misc 1 - h "SWDOWNZ" "DOWNWARD SW FLUX AT SURFACE" "W m-2"
3079	comm_SWDOWNZ(1:ngrid) = fluxsurf_dn_sw_tot(1:ngrid)
3080	!state real TAU_DUST ij misc 1 - h "TAU_DUST" "REFERENCE VISIBLE DUST OPACITY" ""
3081	comm_TAU_DUST(1:ngrid) = tau_pref_gcm(1:ngrid)
3082	!state real RDUST ikj misc 1 - h "RDUST" "DUST RADIUS" "m"
3083	comm_RDUST(1:ngrid,1:nlayer) = rdust(1:ngrid,1:nlayer)
3084	!state real QSURFDUST ij misc 1 - h "QSURFDUST" "DUST MASS AT SURFACE" "kg m-2"
3085	IF (igcm_dust_mass .ne. 0) THEN
3086	comm_QSURFDUST(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
3087	ELSE
3088	comm_QSURFDUST(1:ngrid) = 0.
3089	ENDIF
3090	!state real MTOT ij misc 1 - h "MTOT" "TOTAL MASS WATER VAPOR in pmic" "pmic"
3091	comm_MTOT(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
3092	!state real ICETOT ij misc 1 - h "ICETOT" "TOTAL MASS WATER ICE" "kg m-2"
3093	comm_ICETOT(1:ngrid) = icetot(1:ngrid) * 1.e6 / rho_ice
3094	!state real VMR_ICE ikj misc 1 - h "VMR_ICE" "VOL. MIXING RATIO ICE" "ppm"
3095	IF (igcm_h2o_ice .ne. 0) THEN
3096	comm_VMR_ICE(1:ngrid,1:nlayer) = 1.e6
3097	. * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
3098	. * mmean(1:ngrid,1:nlayer) / mmol(igcm_h2o_ice)
3099	ELSE
3100	comm_VMR_ICE(1:ngrid,1:nlayer) = 0.
3101	ENDIF
3102	!state real TAU_ICE ij misc 1 - h "TAU_ICE" "CLOUD OD at 825 cm-1 TES" ""
3103	if (activice) then
3104	comm_TAU_ICE(1:ngrid) = taucloudtes(1:ngrid)
3105	else
3106	comm_TAU_ICE(1:ngrid) = tauTES(1:ngrid)
3107	endif
3108	!state real RICE ikj misc 1 - h "RICE" "ICE RADIUS" "m"
3109	comm_RICE(1:ngrid,1:nlayer) = rice(1:ngrid,1:nlayer)
3110
3111	!! calculate sensible heat flux in W/m2 for outputs
3112	!! -- the one computed in vdifc is not the real one
3113	!! -- vdifc must have been called
3114	if (.not.callrichsl) then
3115	sensibFlux(1:ngrid) = zflubid(1:ngrid)
3116	. - capcal(1:ngrid)*zdtsdif(1:ngrid)
3117	else
3118	sensibFlux(1:ngrid) =
3119	& (pplay(1:ngrid,1)/(rpt(1:ngrid,1)))cpp
3120	& sqrt(pu(1:ngrid,1)pu(1:ngrid,1)+pv(1:ngrid,1)*pv(1:ngrid,1)
3121	& +(log(1.+0.7wstar(1:ngrid) + 2.3wstar(1:ngrid)2))2)
3122	& zcdh(1:ngrid)(tsurf(1:ngrid)-zh(1:ngrid,1))
3123	endif
3124
3125	#else
3126	#ifndef MESOINI
3127
3128	c ==========================================================
3129	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
3130	c any variable for diagnostic (output with period
3131	c "ecritphy", set in "run.def")
3132	c ==========================================================
3133	c WRITEDIAGFI can ALSO be called from any other subroutines
3134	c for any variables !!
3135	call write_output("emis","Surface emissivity","w.m-1",
3136	& emis(:,iflat))
3137	do islope=1,nslope
3138	write(str2(1:2),'(i2.2)') islope
3139	call write_output("emis_slope"//str2,
3140	& "Surface emissivity","w.m-1",emis(:,islope))
3141	ENDDO
3142	call write_output("zzlev","Interlayer altitude",
3143	& "m",zzlev(:,1:nlayer))
3144	call write_output("pphi","Geopotential","m2s-2",
3145	& pphi(:,:))
3146	call write_output("phisfi","Surface geopotential",
3147	& "m2s-2",phisfi(:))
3148	call write_output("tsurf","Surface temperature","K",
3149	& tsurf(:,iflat))
3150	do islope=1,nslope
3151	write(str2(1:2),'(i2.2)') islope
3152	call write_output("tsurf_slope"//str2,
3153	& "Surface temperature","K",
3154	& tsurf(:,islope))
3155	ENDDO
3156	call write_output("ps","surface pressure","Pa",ps(:))
3157	call write_output("co2ice","co2 ice thickness"
3158	& ,"kg.m-2",qsurf(:,igcm_co2,iflat))
3159	do islope=1,nslope
3160	write(str2(1:2),'(i2.2)') islope
3161	call write_output("co2ice_slope"//str2,"co2 ice thickness"
3162	& ,"kg.m-2",qsurf(:,igcm_co2,islope))
3163	ENDDO
3164	call write_output("watercap","Water ice thickness"
3165	& ,"kg.m-2",watercap(:,iflat))
3166	do islope=1,nslope
3167	write(str2(1:2),'(i2.2)') islope
3168	call write_output("watercap_slope"//str2,
3169	& "Water ice thickness"
3170	& ,"kg.m-2",watercap(:,islope))
3171	ENDDO
3172	call write_output("temp_layer1","temperature in layer 1",
3173	& "K",zt(:,1))
3174	call write_output("temp7","temperature in layer 7",
3175	& "K",zt(:,7))
3176	call write_output("fluxsurf_lw","fluxsurf_lw","W.m-2",
3177	& fluxsurf_lw(:,iflat))
3178	do islope=1,nslope
3179	write(str2(1:2),'(i2.2)') islope
3180	call write_output("fluxsurf_lw_slope"//str2,
3181	& "fluxsurf_lw","W.m-2",
3182	& fluxsurf_lw(:,islope))
3183	ENDDO
3184	call write_output("fluxsurf_dn_sw","fluxsurf_dn_sw",
3185	& "W.m-2",fluxsurf_dn_sw_tot(:,iflat))
3186	do islope=1,nslope
3187	write(str2(1:2),'(i2.2)') islope
3188	call write_output("fluxsurf_dn_sw_slope"//str2,
3189	& "fluxsurf_dn_sw",
3190	& "W.m-2",fluxsurf_dn_sw_tot(:,islope))
3191	ENDDO
3192	call write_output("fluxtop_lw","fluxtop_lw","W.m-2",
3193	& fluxtop_lw(:))
3194	call write_output("fluxtop_up_sw","fluxtop_up_sw","W.m-2",
3195	& fluxtop_up_sw_tot(:))
3196	call write_output("temp","temperature","K",zt(:,:))
3197	call write_output("Sols","Time","sols",zday)
3198	call write_output("Ls","Solar longitude","deg",
3199	& zls*180./pi)
3200	call write_output("u","Zonal wind","m.s-1",zu(:,:))
3201	call write_output("v","Meridional wind","m.s-1",zv(:,:))
3202	call write_output("w","Vertical wind","m.s-1",pw(:,:))
3203	call write_output("rho","density","kg.m-3",rho(:,:))
3204	call write_output("pressure","Pressure","Pa",zplay(:,:))
3205	call write_output("pplev","Pressure","Pa",zplev(:,:))
3206	call write_output('sw_htrt','sw heat. rate',
3207	& 'K/s',zdtsw(:,:))
3208	call write_output('lw_htrt','lw heat. rate',
3209	& 'K/s',zdtlw(:,:))
3210	call write_output("local_time","Local time",
3211	& 'sol',local_time(:))
3212	if (.not.activice) then
3213	CALL write_output('tauTESap',
3214	& 'tau abs 825 cm-1',
3215	& '',tauTES(:))
3216	else
3217	CALL write_output('tauTES',
3218	& 'tau abs 825 cm-1',
3219	& '',taucloudtes(:))
3220	endif
3221	#else
3222	!!! this is to ensure correct initialisation of mesoscale model
3223	call write_output("tsurf","Surface temperature","K",
3224	& tsurf(:,iflat))
3225	call write_output("ps","surface pressure","Pa",ps(:))
3226	call write_output("co2ice","co2 ice thickness","kg.m-2",
3227	& qsurf(:,igcm_co2,iflat))
3228	call write_output("temp","temperature","K",zt(:,:))
3229	call write_output("u","Zonal wind","m.s-1",zu(:,:))
3230	call write_output("v","Meridional wind","m.s-1",zv(:,:))
3231	call write_output("emis","Surface emissivity","w.m-1",
3232	& emis(:,iflat))
3233	call write_output("tsoil","Soil temperature",
3234	& "K",tsoil(:,:,iflat))
3235	call write_output("inertiedat","Soil inertia",
3236	& "K",inertiedat(:,:))
3237	#endif
3238
3239	c ----------------------------------------------------------
3240	c Outputs of the CO2 cycle
3241	c ----------------------------------------------------------
3242
3243	if (igcm_co2.ne.0) then
3244	call write_output("co2","co2 mass mixing ratio",
3245	& "kg.kg-1",zq(:,:,igcm_co2))
3246
3247	if (co2clouds) then
3248	call write_output('ccnqco2','CCNco2 mmr',
3249	& 'kg.kg-1',zq(:,:,igcm_ccnco2_mass))
3250
3251	call write_output('ccnNco2','CCNco2 number',
3252	& 'part.kg-1',zq(:,:,igcm_ccnco2_number))
3253
3254	call write_output('co2_ice','co2_ice mmr in atm',
3255	& 'kg.kg-1',zq(:,:,igcm_co2_ice))
3256
3257	call write_output("mtotco2","total mass atm of co2",
3258	& "kg.m-2",mtotco2(:))
3259	call write_output("icetotco2","total mass atm of co2 ice",
3260	& "kg.m-2", icetotco2(:))
3261	call write_output("vaptotco2","total mass atm of co2 "//
3262	& "vapor","kg.m-2", vaptotco2(:))
3263	if (co2useh2o) then
3264	call write_output('ccnqco2_h2o_m_ice',
3265	& 'CCNco2_h2o_mass_ice mmr',
3266	& 'kg.kg-1',zq(:,:,igcm_ccnco2_h2o_mass_ice))
3267
3268	call write_output('ccnqco2_h2o_m_ccn',
3269	& 'CCNco2_h2o_mass_ccn mmr',
3270	& 'kg.kg-1',zq(:,:,igcm_ccnco2_h2o_mass_ccn))
3271
3272	call write_output('ccnNco2_h2o','CCNco2_h2o number',
3273	& 'part.kg-1',zq(:,:,igcm_ccnco2_h2o_number))
3274	end if
3275
3276	if (meteo_flux) then
3277	call write_output('ccnqco2_meteor','CCNco2_meteor mmr',
3278	& 'kg.kg-1',zq(:,:,igcm_ccnco2_meteor_mass))
3279
3280	call write_output('ccnNco2_meteor','CCNco2_meteor number',
3281	& 'part.kg-1',zq(:,:,igcm_ccnco2_meteor_number))
3282	end if
3283
3284	end if ! of if (co2clouds)
3285	end if ! of if (igcm_co2.ne.0)
3286
3287	! Output He tracer, if there is one
3288	if (igcm_he.ne.0) then
3289	call write_output("he","helium mass mixing ratio",
3290	& "kg/kg",zq(:,:,igcm_he))
3291	vmr = zq(1:ngrid,1:nlayer,igcm_he)
3292	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_he)
3293	call write_output('vmr_he','helium vol. mixing ratio',
3294	& 'mol/mol',vmr(:,:))
3295	end if
3296
3297	c ----------------------------------------------------------
3298	c Outputs of the water cycle
3299	c ----------------------------------------------------------
3300	if (water) then
3301	#ifdef MESOINI
3302	!!!! waterice = q01, voir readmeteo.F90
3303	call write_output('q01',noms(igcm_h2o_ice),
3304	& 'kg/kg',
3305	& zq(:,:,igcm_h2o_ice))
3306	!!!! watervapor = q02, voir readmeteo.F90
3307	call write_output('q02',noms(igcm_h2o_vap),
3308	& 'kg/kg',
3309	& zq(:,:,igcm_h2o_vap))
3310	!!!! surface waterice qsurf02 (voir readmeteo)
3311	call write_output('qsurf02','surface tracer',
3312	& 'kg.m-2',
3313	& qsurf(:,igcm_h2o_ice,iflat))
3314	do islope=1,nslope
3315	write(str2(1:2),'(i2.2)') islope
3316	call write_output('qsurf02_slope'//str2,
3317	& 'surface tracer','kg.m-2',
3318	& qsurf(:,igcm_h2o_ice,islope))
3319	ENDDO
3320	#endif
3321	call write_output('mtot',
3322	& 'total mass of water vapor',
3323	& 'kg/m2',mtot(:))
3324	call write_output('icetot',
3325	& 'total mass of water ice',
3326	& 'kg/m2',icetot(:))
3327	vmr = zq(1:ngrid,1:nlayer,igcm_h2o_ice)
3328	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
3329	call write_output('vmr_h2oice','h2o ice vmr',
3330	& 'mol/mol',vmr(:,:))
3331	vmr = zq(1:ngrid,1:nlayer,igcm_h2o_vap)
3332	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
3333	call write_output('vmr_h2ovap','h2o vap vmr',
3334	& 'mol/mol',vmr(:,:))
3335	call write_output('reffice',
3336	& 'Mean reff',
3337	& 'm',rave(:))
3338	call write_output('h2o_ice','h2o_ice','kg/kg',
3339	& zq(:,:,igcm_h2o_ice))
3340	call write_output('h2o_vap','h2o_vap','kg/kg',
3341	& zq(:,:,igcm_h2o_vap))
3342
3343	if (hdo) then
3344	vmr=zq(1:ngrid,1:nlayer,igcm_hdo_ice)
3345	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_ice)
3346	CALL write_output('vmr_hdoice','hdo ice vmr',
3347	& 'mol/mol',vmr(:,:))
3348	vmr=zq(1:ngrid,1:nlayer,igcm_hdo_vap)
3349	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_vap)
3350	CALL write_output('vmr_hdovap','hdo vap vmr',
3351	& 'mol/mol',vmr(:,:))
3352	call write_output('hdo_ice','hdo_ice','kg/kg',
3353	& zq(:,:,igcm_hdo_ice))
3354	call write_output('hdo_vap','hdo_vap','kg/kg',
3355	& zq(:,:,igcm_hdo_vap))
3356
3357	CALL write_output('mtotD',
3358	& 'total mass of HDO vapor',
3359	& 'kg/m2',mtotD(:))
3360	CALL write_output('icetotD',
3361	& 'total mass of HDO ice',
3362	& 'kg/m2',icetotD(:))
3363
3364	C Calculation of the D/H ratio
3365	do l=1,nlayer
3366	do ig=1,ngrid
3367	if (zq(ig,l,igcm_h2o_vap).gt.qperemin) then
3368	DoH_vap(ig,l) = ( zq(ig,l,igcm_hdo_vap)/
3369	& zq(ig,l,igcm_h2o_vap) )1./(2.155.76e-6)
3370	else
3371	DoH_vap(ig,l) = 0.
3372	endif
3373	enddo
3374	enddo
3375
3376	do l=1,nlayer
3377	do ig=1,ngrid
3378	if (zq(ig,l,igcm_h2o_ice).gt.qperemin) then
3379	DoH_ice(ig,l) = ( zq(ig,l,igcm_hdo_ice)/
3380	& zq(ig,l,igcm_h2o_ice) )/(2.*155.76e-6)
3381	else
3382	DoH_ice(ig,l) = 0.
3383	endif
3384	enddo
3385	enddo
3386
3387	CALL write_output('DoH_vap',
3388	& 'D/H ratio in vapor',
3389	& ' ',DoH_vap(:,:))
3390	CALL write_output('DoH_ice',
3391	& 'D/H ratio in ice',
3392	& '',DoH_ice(:,:))
3393
3394	endif !hdo
3395
3396	!A. Pottier
3397	CALL write_output('rmoym',
3398	& 'alternative reffice',
3399	& 'm',rave2(:))
3400	call write_output('saturation',
3401	& 'h2o vap saturation ratio','dimless',satu(:,:))
3402	if (scavenging) then
3403	CALL write_output("Nccntot",
3404	& "condensation nuclei","Nbr/m2",
3405	& Nccntot(:))
3406	CALL write_output("Mccntot",
3407	& "mass condensation nuclei","kg/m2",
3408	& Mccntot(:))
3409	endif
3410	call write_output('rice','Ice particle size',
3411	& 'm',rice(:,:))
3412	call write_output('h2o_ice_s',
3413	& 'surface h2o_ice',
3414	& 'kg.m-2',qsurf(:,igcm_h2o_ice,iflat))
3415	do islope=1,nslope
3416	write(str2(1:2),'(i2.2)') islope
3417	call write_output('h2o_ice_s_slope'//str2,
3418	& 'surface h2o_ice',
3419	& 'kg.m-2',qsurf(:,igcm_h2o_ice,islope))
3420	ENDDO
3421	if (hdo) then
3422	call write_output('hdo_ice_s',
3423	& 'surface hdo_ice',
3424	& 'kg.m-2',qsurf(:,igcm_hdo_ice,iflat))
3425	do islope=1,nslope
3426	write(str2(1:2),'(i2.2)') islope
3427	call write_output('hdo_ice_s_slope'//str2,
3428	& 'surface hdo_ice',
3429	& 'kg.m-2',qsurf(:,igcm_hdo_ice,islope))
3430	ENDDO
3431
3432	do ig=1,ngrid
3433	if (qsurf_meshavg(ig,igcm_h2o_ice).gt.qperemin) then
3434	DoH_surf(ig) = 0.5*( qsurf_meshavg(ig,igcm_hdo_ice)/
3435	& qsurf_meshavg(ig,igcm_h2o_ice) )/155.76e-6
3436	else
3437	DoH_surf(ig) = 0.
3438	endif
3439	enddo
3440
3441	call write_output('DoH_surf',
3442	& 'surface D/H',
3443	& '',DoH_surf(:))
3444	endif ! hdo
3445
3446	CALL write_output('albedo',
3447	& 'albedo',
3448	& '',albedo(:,1,iflat))
3449	do islope=1,nslope
3450	write(str2(1:2),'(i2.2)') islope
3451	CALL write_output('albedo_slope'//str2,
3452	& 'albedo',
3453	& '',albedo(:,1,islope))
3454	ENDDO
3455	if (tifeedback) then
3456	call write_output("soiltemp",
3457	& "Soil temperature","K",
3458	& tsoil(:,:,iflat))
3459	call write_output('soilti',
3460	& 'Soil Thermal Inertia',
3461	& 'J.s-1/2.m-2.K-1',inertiesoil(:,:,iflat))
3462	do islope=1,nslope
3463	write(str2(1:2),'(i2.2)') islope
3464	call write_output('soilti_slope'//str2,
3465	& 'Soil Thermal Inertia',
3466	& 'J.s-1/2.m-2.K-1',inertiesoil(:,:,islope))
3467	ENDDO
3468	endif
3469	!A. Pottier
3470	if (CLFvarying) then !AP14 nebulosity
3471	call write_output('totcloudfrac',
3472	& 'Total cloud fraction',
3473	& ' ',totcloudfrac(:))
3474	end if !clf varying
3475	end if !(water)
3476
3477	c ----------------------------------------------------------
3478	c Outputs of the dust cycle
3479	c ----------------------------------------------------------
3480
3481	call write_output('tau_pref_scenario',
3482	& 'Prescribed visible dust optical depth at 610Pa',
3483	& 'NU',tau_pref_scenario(:))
3484
3485	call write_output('tau_pref_gcm',
3486	& 'Visible dust optical depth at 610Pa in the GCM',
3487	& 'NU',tau_pref_gcm(:))
3488
3489	if (reff_driven_IRtoVIS_scenario) then
3490	call write_output('IRtoVIScoef',
3491	& 'conversion coeff for dust tau from abs9.3um to ext0.67um',
3492	& '/',IRtoVIScoef(:))
3493	endif
3494
3495	if (dustbin.ne.0) then
3496
3497	call write_output('tau','taudust','SI',tau(:,1))
3498
3499	#ifndef MESOINI
3500	if (doubleq) then
3501	call write_output('dqsdust',
3502	& 'deposited surface dust mass',
3503	& 'kg.m-2.s-1',dqdustsurf(:))
3504	call write_output('dqndust',
3505	& 'deposited surface dust number',
3506	& 'number.m-2.s-1',dndustsurf(:))
3507	call write_output('reffdust','reffdust',
3508	& 'm',rdust(:,:)*ref_r0)
3509	call write_output('dustq','Dust mass mr',
3510	& 'kg/kg',qdust(:,:))
3511	call write_output('dustN','Dust number',
3512	& 'part/kg',ndust(:,:))
3513
3514	select case (trim(dustiropacity))
3515	case ("tes")
3516	call write_output('dsodust',
3517	& 'density scaled extinction opacity of std dust at 9.3um(TES)',
3518	& 'm2.kg-1',dsodust(:,:))
3519	call write_output('dso',
3520	& 'density scaled extinction opacity of all dust at 9.3um(TES)',
3521	& 'm2.kg-1',dsodust(:,:)+dsords(:,:)+dsotop(:,:))
3522	case ("mcs")
3523	call write_output('dsodust',
3524	& 'density scaled extinction opacity of std dust at 21.6um(MCS)',
3525	& 'm2.kg-1',dsodust(:,:))
3526	call write_output('dso',
3527	& 'density scaled extinction opacity of all dust at 21.6um(MCS)',
3528	& 'm2.kg-1',dsodust(:,:)+dsords(:,:)+dsotop(:,:))
3529	end select
3530	else ! (doubleq=.false.)
3531	do iq=1,dustbin
3532	write(str2(1:2),'(i2.2)') iq
3533	call write_output('q'//str2,'mix. ratio',
3534	& 'kg/kg',zq(:,:,iq))
3535	call write_output('qsurf'//str2,'qsurf',
3536	& 'kg.m-2',qsurf(:,iq,iflat))
3537	do islope=1,nslope
3538	write(str2(1:2),'(i2.2)') islope
3539	call write_output('qsurf_slope'//str2,'qsurf',
3540	& 'kg.m-2',qsurf(:,iq,islope))
3541	ENDDO
3542	end do
3543	endif ! (doubleq)
3544
3545	if (rdstorm) then ! writediagfi tendencies stormdust tracers
3546	call write_output('reffstormdust','reffstormdust',
3547	& 'm',rstormdust(:,:)*ref_r0)
3548	call write_output('mstormdtot',
3549	& 'total mass of stormdust only',
3550	& 'kg.m-2',mstormdtot(:))
3551	call write_output('mdusttot',
3552	& 'total mass of dust only',
3553	& 'kg.m-2',mdusttot(:))
3554	call write_output('rdsdqsdust',
3555	& 'deposited surface stormdust mass',
3556	& 'kg.m-2.s-1',rdsdqdustsurf(:))
3557	call write_output('rdsdustq','storm Dust mass mr',
3558	& 'kg/kg',rdsqdust(:,:))
3559	call write_output('rdsdustqmodel','storm Dust massmr',
3560	& 'kg/kg',pq(:,:,igcm_stormdust_mass))
3561	call write_output('rdsdustN','storm Dust number',
3562	& 'part/kg',rdsndust(:,:))
3563	call write_output("stormfract",
3564	& "fraction of the mesh, with stormdust","none",
3565	& totstormfract(:))
3566	call write_output('qsurf',
3567	& 'stormdust injection',
3568	& 'kg.m-2',qsurf(:,igcm_stormdust_mass,iflat))
3569	do islope=1,nslope
3570	write(str2(1:2),'(i2.2)') islope
3571	call write_output('qsurf_slope'//str2,
3572	& 'stormdust injection',
3573	& 'kg.m-2',qsurf(:,igcm_stormdust_mass,islope))
3574	ENDDO
3575	call write_output('pdqsurf',
3576	& 'tendancy stormdust mass at surface',
3577	& 'kg.m-2',dqsurf(:,igcm_stormdust_mass,iflat))
3578	do islope=1,nslope
3579	write(str2(1:2),'(i2.2)') islope
3580	call write_output('pdqsurf_slope'//str2,
3581	& 'tendancy stormdust mass at surface',
3582	& 'kg.m-2',dqsurf(:,igcm_stormdust_mass,islope))
3583	ENDDO
3584	call write_output('wspeed','vertical speed stormdust',
3585	& 'm/s',wspeed(:,:))
3586	call write_output('zdqsed_dustq'
3587	& ,'sedimentation q','kg.m-2.s-1',
3588	& zdqsed(:,:,igcm_dust_mass))
3589	call write_output('zdqssed_dustq'
3590	& ,'sedimentation q','kg.m-2.s-1',
3591	& zdqssed(:,igcm_dust_mass))
3592	call write_output('zdqsed_rdsq'
3593	& ,'sedimentation q','kg.m-2.s-1',
3594	& zdqsed(:,:,igcm_stormdust_mass))
3595	call write_output('zdqsed_dustN'
3596	& ,'sedimentation N','Nbr.m-2.s-1',
3597	& zdqsed(:,:,igcm_dust_number))
3598	call write_output('rdust','rdust',
3599	& 'm',rdust(:,:))
3600	call write_output('rstormdust','rstormdust',
3601	& 'm',rstormdust(:,:))
3602	call write_output('totaldustq','total dust mass',
3603	& 'kg/kg',qdusttotal(:,:))
3604
3605	select case (trim(dustiropacity))
3606	case ("tes")
3607	call write_output('dsords',
3608	& 'density scaled extinction opacity of stormdust at 9.3um(TES)',
3609	& 'm2.kg-1',dsords(:,:))
3610	case ("mcs")
3611	call write_output('dsords',
3612	& 'density scaled extinction opacity of stormdust at 21.6um(MCS)',
3613	& 'm2.kg-1',dsords(:,:))
3614	end select
3615	endif ! (rdstorm)
3616
3617	if (topflows) then
3618	call write_output('refftopdust','refftopdust',
3619	& 'm',rtopdust(:,:)*ref_r0)
3620	call write_output('topdustq','top Dust mass mr',
3621	& 'kg/kg',pq(:,:,igcm_topdust_mass))
3622	call write_output('topdustN','top Dust number',
3623	& 'part/kg',pq(:,:,igcm_topdust_number))
3624	select case (trim(dustiropacity))
3625	case ("tes")
3626	call write_output('dsotop',
3627	& 'density scaled extinction opacity of topdust at 9.3um(TES)',
3628	& 'm2.kg-1',dsotop(:,:))
3629	case ("mcs")
3630	call write_output('dsotop',
3631	& 'density scaled extinction opacity of topdust at 21.6um(MCS)',
3632	& 'm2.kg-1',dsotop(:,:))
3633	end select
3634	endif ! (topflows)
3635
3636	if (dustscaling_mode==2) then
3637	call write_output("dust_rad_adjust",
3638	& "radiative adjustment coefficient for dust",
3639	& "",dust_rad_adjust(:))
3640	endif
3641
3642	if (scavenging) then
3643	call write_output('ccnq','CCN mass mr',
3644	& 'kg/kg',qccn(:,:))
3645	call write_output('ccnN','CCN number',
3646	& 'part/kg',nccn(:,:))
3647	call write_output('surfccnq','Surf nuclei mass mr',
3648	& 'kg.m-2',qsurf(:,igcm_ccn_mass,iflat))
3649	do islope=1,nslope
3650	write(str2(1:2),'(i2.2)') islope
3651	call write_output('surfccnq_slope'//str2,
3652	& 'Surf nuclei mass mr',
3653	& 'kg.m-2',qsurf(:,igcm_ccn_mass,islope))
3654	ENDDO
3655	call write_output('surfccnN','Surf nuclei number',
3656	& 'kg.m-2',qsurf(:,igcm_ccn_number,iflat))
3657	do islope=1,nslope
3658	write(str2(1:2),'(i2.2)') islope
3659	call write_output('surfccnN_slope'//str2,
3660	& 'Surf nuclei number',
3661	& 'kg.m-2',qsurf(:,igcm_ccn_number,islope))
3662	ENDDO
3663	endif ! (scavenging)
3664
3665	#else
3666	! !!! to initialize mesoscale we need scaled variables
3667	! !!! because this must correspond to starting point for tracers
3668	! call write_output('dustq','Dust mass mr',
3669	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass))
3670	! call write_output('dustN','Dust number',
3671	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number))
3672	! call write_output('ccn','Nuclei mass mr',
3673	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass))
3674	! call write_output('ccnN','Nuclei number',
3675	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number))
3676	if (freedust) then
3677	call write_output('dustq','Dust mass mr',
3678	& 'kg/kg',qdust)
3679	call write_output('dustN','Dust number',
3680	& 'part/kg',ndust)
3681	call write_output('ccn','CCN mass mr',
3682	& 'kg/kg',qccn)
3683	call write_output('ccnN','CCN number',
3684	& 'part/kg',nccn)
3685	else
3686	call write_output('dustq','Dust mass mr',
3687	& 'kg/kg',pq(:,:,igcm_dust_mass))
3688	call write_output('dustN','Dust number',
3689	& 'part/kg',pq(:,:,igcm_dust_number))
3690	call write_output('ccn','Nuclei mass mr',
3691	& 'kg/kg',pq(:,:,igcm_ccn_mass))
3692	call write_output('ccnN','Nuclei number',
3693	& 'part/kg',pq(:,:,igcm_ccn_number))
3694	endif
3695	#endif
3696
3697	end if ! (dustbin.ne.0)
3698
3699	c ----------------------------------------------------------
3700	c GW non-oro outputs
3701	c ----------------------------------------------------------
3702
3703	if(calllott_nonoro) then
3704	call write_output("dugwno","GW non-oro dU","m/s2",
3705	$ d_u_hin(:,:)/ptimestep)
3706	call write_output("dvgwno","GW non-oro dV","m/s2",
3707	$ d_v_hin(:,:)/ptimestep)
3708	endif !(calllott_nonoro)
3709
3710	c ----------------------------------------------------------
3711	c Thermospheric outputs
3712	c ----------------------------------------------------------
3713
3714	if(callthermos) then
3715
3716	call write_output("q15um","15 um cooling","K/s",
3717	$ zdtnlte(:,:))
3718	call write_output("quv","UV heating","K/s",
3719	$ zdteuv(:,:))
3720	call write_output("cond","Thermal conduction","K/s",
3721	$ zdtconduc(:,:))
3722	call write_output("qnir","NIR heating","K/s",
3723	$ zdtnirco2(:,:))
3724
3725	!H, H2 and D escape fluxes
3726
3727	call write_output("PhiH","H escape flux","s-1",
3728	$ PhiEscH)
3729	call write_output("PhiH2","H2 escape flux","s-1",
3730	$ PhiEscH2)
3731	call write_output("PhiD","D escape flux","s-1",
3732	$ PhiEscD)
3733
3734	endif !(callthermos)
3735
3736	call write_output("q15um","15 um cooling","K/s",
3737	$ zdtnlte(:,:))
3738	call write_output("qnir","NIR heating","K/s",
3739	$ zdtnirco2(:,:))
3740
3741	c ----------------------------------------------------------
3742	c ----------------------------------------------------------
3743	c PBL OUTPUS
3744	c ----------------------------------------------------------
3745	c ----------------------------------------------------------
3746
3747	c ----------------------------------------------------------
3748	c Outputs of thermals
3749	c ----------------------------------------------------------
3750	if (calltherm) then
3751	call write_output('lmax_th',
3752	& 'hauteur du thermique','point',
3753	& lmax_th_out(:))
3754	call write_output('zmax_th',
3755	& 'hauteur du thermique','m',
3756	& zmax_th(:))
3757	call write_output('hfmax_th',
3758	& 'maximum TH heat flux','K.m/s',
3759	& hfmax_th(:))
3760	call write_output('wstar',
3761	& 'maximum TH vertical velocity','m/s',
3762	& wstar(:))
3763	end if
3764
3765	c ----------------------------------------------------------
3766	c ----------------------------------------------------------
3767	c END OF PBL OUTPUS
3768	c ----------------------------------------------------------
3769	c ----------------------------------------------------------
3770
3771
3772	c ----------------------------------------------------------
3773	c Output in netcdf file "diagsoil.nc" for subterranean
3774	c variables (output every "ecritphy", as for writediagfi)
3775	c ----------------------------------------------------------
3776
3777	! Write soil temperature
3778	call write_output("soiltemp","Soil temperature","K",
3779	& tsoil(:,:,iflat))
3780	do islope=1,nslope
3781	write(str2(1:2),'(i2.2)') islope
3782	call write_output("soiltemp_slope"//str2,
3783	& "Soil temperature","K",
3784	& tsoil(:,:,islope))
3785	ENDDO
3786
3787	!PREVIOUSLY IN 1D ONLY
3788	call write_output("dtrad","rad. heat. rate",
3789	& "K.s-1",dtrad(:,:))
3790
3791	if (rdstorm) then
3792	call write_output('aerosol_dust','opacity of env. dust',''
3793	& ,aerosol(:,:,iaer_dust_doubleq))
3794	call write_output('aerosol_stormdust',
3795	& 'opacity of storm dust',''
3796	& ,aerosol(:,:,iaer_stormdust_doubleq))
3797	call write_output('dqsdifdustq','diffusion',
3798	& 'kg.m-2.s-1',zdqsdif(:,igcm_dust_mass,iflat))
3799	do islope=1,nslope
3800	write(str2(1:2),'(i2.2)') islope
3801	call write_output('dqsdifdustq_slope'//str2,
3802	& 'diffusion',
3803	& 'kg.m-2.s-1',zdqsdif(:,igcm_dust_mass,islope))
3804	ENDDO
3805	call write_output('dqsdifrdsq','diffusion',
3806	& 'kg.m-2.s-1',zdqsdif(:,igcm_stormdust_mass,iflat))
3807	do islope=1,nslope
3808	write(str2(1:2),'(i2.2)') islope
3809	call write_output('dqsdifrdsq_slope'//str2,
3810	& 'diffusion',
3811	& 'kg.m-2.s-1',zdqsdif(:,igcm_stormdust_mass,islope))
3812	ENDDO
3813	endif !(rdstorm)
3814
3815	if(water) then
3816	if (.not.scavenging) then
3817	call write_output('zdqcloud_ice','cloud ice',
3818	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_ice))
3819	call write_output('zdqcloud_vap','cloud vap',
3820	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_vap))
3821	call write_output('zdqcloud','cloud ice',
3822	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_ice)
3823	& +zdqcloud(:,:,igcm_h2o_vap))
3824	IF (hdo) THEN
3825	call write_output('zdqcloud_iceD','cloud ice hdo',
3826	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_hdo_ice))
3827	call write_output('zdqcloud_vapD','cloud vap hdo',
3828	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_hdo_vap))
3829	ENDIF ! hdo
3830	endif !not.scavenging
3831	ENDIF ! of IF (water)
3832	!PREVIOUSLY IN 1D ONLY
3833
3834	c ==========================================================
3835	c END OF WRITEDIAGFI
3836	c ==========================================================
3837	#endif
3838	! of ifdef MESOSCALE
3839
3840	c ELSE ! if(ngrid.eq.1)
3841
3842	c#ifndef MESOSCALE
3843	c write(*,
3844	c & '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)')
3845	c & zls*180./pi,odpref,tau_pref_scenario
3846	c#endif
3847
3848	c END IF ! if(ngrid.ne.1)
3849
3850
3851	! test for co2 conservation with co2 microphysics
3852	if (igcm_co2_ice.ne.0) then
3853	co2totB = 0. ! added by C.M.
3854	do ig=1,ngrid
3855	do l=1,nlayer
3856	co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
3857	& (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice)
3858	& +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
3859	enddo
3860	co2totB = co2totB + qsurf(ig,igcm_co2,iflat)
3861	enddo
3862	else
3863	co2totB = 0. ! added by C.M.
3864	do ig=1,ngrid
3865	do l=1,nlayer
3866	co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
3867	& (pq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2)*ptimestep)
3868	enddo
3869	co2totB = co2totB + qsurf(ig,igcm_co2,iflat)
3870	enddo
3871	endif ! of if (igcm_co2_ice.ne.0)
3872	co2conservation = (co2totA-co2totB)/co2totA
3873	call write_output( 'co2conservation',
3874	& 'Total CO2 mass conservation in physic',
3875	& ' ', co2conservation)
3876	! XIOS outputs
3877	#ifdef CPP_XIOS
3878	! Send fields to XIOS: (NB these fields must also be defined as
3879	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
3880	CALL send_xios_field("ls",zls*180./pi)
3881
3882	CALL send_xios_field("controle",tab_cntrl_mod,1)
3883
3884	CALL send_xios_field("ap",ap,1)
3885	CALL send_xios_field("bp",bp,1)
3886	CALL send_xios_field("aps",aps,1)
3887	CALL send_xios_field("bps",bps,1)
3888
3889	CALL send_xios_field("phisinit",phisfi)
3890
3891	CALL send_xios_field("ps",ps)
3892	CALL send_xios_field("area",cell_area)
3893
3894	! CALL send_xios_field("ISR",fluxtop_dn_sw_tot)
3895	CALL send_xios_field("OLR",fluxtop_lw)
3896
3897	CALL send_xios_field("tsurf",tsurf(:,iflat))
3898	! CALL send_xios_field("inertiedat",inertiedat)
3899	CALL send_xios_field("tsoil",tsoil(:,:,iflat))
3900	CALL send_xios_field("co2ice",qsurf(:,igcm_co2,iflat))
3901
3902	! CALL send_xios_field("temp",zt)
3903	CALL send_xios_field("u",zu)
3904	CALL send_xios_field("v",zv)
3905
3906	! CALL send_xios_field("rho",rho)
3907	! Orographic Gravity waves tendencies
3908	! if (calllott) then
3909	! CALL send_xios_field("dugw",zdugw/ptimestep)
3910	! CALL send_xios_field("dvgw",zdvgw/ptimestep)
3911	! CALL send_xios_field("dtgw",zdtgw/ptimestep)
3912	! endif
3913	!CREATE IF CO2CYCLE
3914	! if (igcm_co2.ne.0) then
3915	! CALL send_xios_field("co2",zq(:,:,igcm_co2))
3916	! endif
3917	! Water cycle
3918	! if (water) then
3919	! CALL send_xios_field("watercap",watercap)
3920	!CALL send_xios_field("watercaptag",watercaptag)
3921	! CALL send_xios_field("mtot",mtot)
3922	! CALL send_xios_field("icetot",icetot)
3923	! if (igcm_h2o_vap.ne.0 .and. igcm_h2o_ice.ne.0) then
3924	! CALL send_xios_field("h2o_vap",zq(:,:,igcm_h2o_vap))
3925	! CALL send_xios_field("h2o_ice",zq(:,:,igcm_h2o_ice))
3926	! endif
3927	! endif
3928	! if (.not.activice) then
3929	! CALL send_xios_field("tauTESap",tauTES)
3930	! else
3931	! CALL send_xios_field("tauTES",taucloudtes)
3932	! endif
3933
3934	! CALL send_xios_field("h2o_ice_s",qsurf(:,igcm_h2o_ice))
3935
3936
3937	if (lastcall.and.is_omp_master) then
3938	write(,) "physiq lastcall: call xios_context_finalize"
3939	call xios_context_finalize
3940	endif
3941	#endif
3942
3943	if (check_physics_outputs) then
3944	! Check the validity of updated fields at the end of the physics step
3945	call check_physics_fields("end of physiq:",zt,zu,zv,zplev,zq)
3946	endif
3947
3948	icount=icount+1
3949
3950	END SUBROUTINE physiq
3951
3952	END MODULE physiq_mod

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

Download in other formats:

Original Format