Context Navigation

physiq_mod.F @ 2942

Last change on this file since 2942 was 2942, checked in by llange, 20 months ago
Mars PCM Add the possibility to use a different thermal inertia (field 'inertiesoil') than inertiedat in the PCM (for paleoclimate studies). By defaut, if there is not inertiesoil, inertiedat is used. Soil_tifeedback still work with inertiedat Newstart adapted, start2archive will be modified by Romain Vandemeulebrouck. LL
File size: 155.2 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, inertiesoil,! dat: soil thermal inertia for present climate, inertiesoil is the TI read in the start
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_tifeedback(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	integer :: islope
543	real :: zdqsdif_meshavg_tmp(ngrid,nq) ! temporary for dust lifting
544
545	logical :: write_restart
546
547	c=======================================================================
548	pdq(:,:,:) = 0.
549
550	c 1. Initialisation:
551	c -----------------
552	c 1.1 Initialisation only at first call
553	c ---------------------------------------
554
555	IF (firstcall) THEN
556
557	call getin_p("check_physics_inputs",check_physics_inputs)
558	call getin_p("check_physics_outputs",check_physics_outputs)
559
560	c variables set to 0
561	c ~~~~~~~~~~~~~~~~~~
562	aerosol(:,:,:)=0
563	dtrad(:,:)=0
564
565	#ifndef MESOSCALE
566	fluxrad(:,:)=0
567	wstar(:)=0.
568	#endif
569
570	#ifdef CPP_XIOS
571	! Initialize XIOS context
572	write(,) "physiq: call wxios_context_init"
573	CALL wxios_context_init
574	#endif
575
576	c read startfi
577	c ~~~~~~~~~~~~
578	#ifndef MESOSCALE
579
580	! GCM. Read netcdf initial physical parameters.
581	CALL phyetat0 ("startfi.nc",0,0,
582	& nsoilmx,ngrid,nlayer,nq,
583	& day_ini,time_phys,
584	& tsurf,tsoil,albedo,emis,
585	& q2,qsurf,tauscaling,totcloudfrac,wstar,
586	& watercap,def_slope,def_slope_mean,subslope_dist)
587
588	! Sky view:
589	DO islope=1,nslope
590	sky_slope(islope) = (1.+cos(pi*def_slope_mean(islope)/180.))/2.
591	END DO
592	! Determine the 'flatest' slopes
593	iflat = 1
594	DO islope=2,nslope
595	IF(abs(def_slope_mean(islope)).lt.
596	& abs(def_slope_mean(iflat)))THEN
597	iflat = islope
598	ENDIF
599	ENDDO
600	PRINT*,'Flat slope for islope = ',iflat
601	PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
602
603	#else
604	! MESOSCALE. Supposedly everything is already set in modules.
605	! So we just check. And we fill day_ini
606	print*,"check: --- in physiq.F"
607	print*,"check: rad,cpp,g,r,rcp,daysec"
608	print*,rad,cpp,g,r,rcp,daysec
609	PRINT*,'check: tsurf ',tsurf(1,:),tsurf(ngrid,:)
610	PRINT*,'check: tsoil ',tsoil(1,1,:),tsoil(ngrid,nsoilmx,:)
611	PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)
612	PRINT*,'check: midlayer,layer ', mlayer(:),layer(:)
613	PRINT*,'check: tracernames ', noms
614	PRINT*,'check: emis ',emis(1,:),emis(ngrid,:)
615	PRINT*,'check: q2 ',q2(1,1),q2(ngrid,nlayer+1)
616	PRINT*,'check: qsurf ',qsurf(1,1,:),qsurf(ngrid,nq,:)
617	PRINT*,'check: co2ice ',qsurf(1,igcm_co2,:),qsurf(ngrid,igcm_co2,:)
618	!!!
619	day_ini = pday
620	!!! a couple initializations (dummy for mesoscale) done in phyetat0
621	!!! --- maybe this should be done in update_inputs_physiq_mod
622
623	tauscaling(:)=1.0 !! probably important
624	totcloudfrac(:)=1.0
625	DO islope = 1,nslope
626	albedo(:,1,islope)=albedodat(:)
627	albedo(:,2,islope)=albedo(:,1,islope)
628	inertiesoil(:,:,islope) = inertiedat(:,:)
629	watercap(:,:)=0.0
630	ENDDO
631	#endif
632	#ifndef MESOSCALE
633	if (.not.startphy_file) then
634	! starting without startfi.nc and with callsoil
635	! is not yet possible as soildepth default is not defined
636	if (callsoil) then
637	! default mlayer distribution, following a power law:
638	! mlayer(k)=lay1alpha*(k-1/2)
639	lay1=2.e-4
640	alpha=2
641	do iloop=0,nsoilmx-1
642	mlayer(iloop)=lay1(alpha*(iloop-0.5))
643	enddo
644	lay1=sqrt(mlayer(0)*mlayer(1))
645	alpha=mlayer(1)/mlayer(0)
646	do iloop=1,nsoilmx
647	layer(iloop)=lay1(alpha*(iloop-1))
648	enddo
649	endif
650	! additionnal "academic" initialization of physics
651	do islope = 1,nslope
652	tsurf(:,islope)=pt(:,1)
653	enddo
654	write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
655	do isoil=1,nsoilmx
656	tsoil(1:ngrid,isoil,:)=tsurf(1:ngrid,:)
657	enddo
658	write(,) "Physiq: initializing inertiedat !!"
659	inertiedat(:,:)=400.
660	inertiesoil(:,:,:)=400.
661	write(,) "Physiq: initializing day_ini to pdat !"
662	day_ini=pday
663	endif
664	#endif
665	if (pday.ne.day_ini) then
666	write(,) "PHYSIQ: ERROR: bad synchronization between ",
667	& "physics and dynamics"
668	write(,) "dynamics day [pday]: ",pday
669	write(,) "physics day [day_ini]: ",day_ini
670	call abort_physic("physiq","dynamics day /= physics day",1)
671	endif
672
673	write (,) 'In physiq day_ini =', day_ini
674
675	c initialize tracers
676	c ~~~~~~~~~~~~~~~~~~
677
678	CALL initracer(ngrid,nq,qsurf)
679
680	c Initialize albedo and orbital calculation
681	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
682	CALL surfini(ngrid,qsurf)
683
684	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
685
686	c initialize soil
687	c ~~~~~~~~~~~~~~~
688	IF (callsoil) THEN
689	c Thermal inertia feedback:
690	IF (tifeedback) THEN
691	DO islope = 1,nslope
692	CALL soil_tifeedback(ngrid,nsoilmx,
693	s qsurf(:,:,islope),
694	s inertiesoil_tifeedback(:,:,islope))
695	ENDDO
696	CALL soil(ngrid,nsoilmx,firstcall,
697	s inertiesoil_tifeedback,
698	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
699	ELSE
700	CALL soil(ngrid,nsoilmx,firstcall,inertiesoil,
701	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
702	ENDIF ! of IF (tifeedback)
703	ELSE
704	PRINT*,
705	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
706	DO ig=1,ngrid
707	capcal(ig,:)=1.e5
708	fluxgrd(ig,:)=0.
709	ENDDO
710	ENDIF
711	icount=1
712
713	#ifndef MESOSCALE
714	c Initialize thermospheric parameters
715	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
716
717	if (callthermos) then
718	call fill_data_thermos
719	call allocate_param_thermos(nlayer)
720	call allocate_param_iono(nlayer,nreact)
721	call param_read_e107
722	endif
723	#endif
724	c Initialize R and Cp as constant
725
726	if (.not.callthermos .and. .not.photochem) then
727	do l=1,nlayer
728	do ig=1,ngrid
729	rnew(ig,l)=r
730	cpnew(ig,l)=cpp
731	mmean(ig,l)=mugaz
732	enddo
733	enddo
734	endif
735
736	if(callnlte.and.nltemodel.eq.2) call nlte_setup
737	if(callnirco2.and.nircorr.eq.1) call NIR_leedat
738
739
740	IF (water.AND.(ngrid.NE.1)) THEN
741	write(,)"physiq: water_param Surface water frost albedo:",
742	. albedo_h2o_frost
743	write(,)"physiq: water_param Surface watercap albedo:",
744	. albedo_h2o_cap
745	ENDIF
746
747	#ifndef MESOSCALE
748	if(ngrid.ne.1) then
749	! no need to compute slopes when in 1D; it is an input
750	if (callslope) call getslopes(ngrid,phisfi)
751	endif !1D
752	if (ecritstart.GT.0) then
753	call physdem0("restartfi.nc",longitude,latitude,
754	& nsoilmx,ngrid,nlayer,nq,
755	& ptimestep,pday,0.,cell_area,
756	& albedodat,inertiedat,def_slope,
757	& subslope_dist)
758	else
759	call physdem0("restartfi.nc",longitude,latitude,
760	& nsoilmx,ngrid,nlayer,nq,
761	& ptimestep,float(day_end),0.,cell_area,
762	& albedodat,inertiedat,def_slope,
763	& subslope_dist)
764	endif
765
766	c Initialize mountain mesh fraction for the entrainment by top flows param.
767	c ~~~~~~~~~~~~~~~
768	if (topflows) call topmons_setup(ngrid)
769
770	#endif
771
772	#ifdef CPP_XIOS
773	! XIOS outputs
774	write(,) "physiq firstcall: call initialize_xios_output"
775	call initialize_xios_output(pday,ptime,ptimestep,daysec,
776	& presnivs,pseudoalt,mlayer)
777	#endif
778	ENDIF ! (end of "if firstcall")
779
780	if (check_physics_inputs) then
781	! Check the validity of input fields coming from the dynamics
782	call check_physics_fields("begin physiq:",pt,pu,pv,pplev,pq)
783	endif
784
785	c ---------------------------------------------------
786	c 1.2 Initializations done at every physical timestep:
787	c ---------------------------------------------------
788	c
789
790	#ifdef CPP_XIOS
791	! update XIOS time/calendar
792	call update_xios_timestep
793	#endif
794
795
796
797	c Initialize various variables
798	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
799	pdv(:,:)=0
800	pdu(:,:)=0
801	pdt(:,:)=0
802	pdq(:,:,:)=0
803	pdpsrf(:)=0
804	zflubid(:,:)=0
805	zdtsurf(:,:)=0
806	dqsurf(:,:,:)=0
807	dsodust(:,:)=0.
808	dsords(:,:)=0.
809	dsotop(:,:)=0.
810	dwatercap(:,:)=0
811
812	call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf,
813	& albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg)
814
815	! Dust scenario conversion coefficient from IRabs to VISext
816	IRtoVIScoef(1:ngrid)=2.6 ! initialized with former value from Montabone et al 2015
817	! recomputed in aeropacity if reff_driven_IRtoVIS_scenario=.true.
818
819	#ifdef DUSTSTORM
820	pq_tmp(:,:,:)=0
821	#endif
822	igout=ngrid/2+1
823
824
825	zday=pday+ptime ! compute time, in sols (and fraction thereof)
826	! Compute local time at each grid point
827	DO ig=1,ngrid
828	local_time(ig)=modulo(1.+(zday-INT(zday))
829	& +(longitude_deg(ig)/15)/24,1.)
830	ENDDO
831	c Compute Solar Longitude (Ls) :
832	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
833	if (season) then
834	call solarlong(zday,zls)
835	else
836	call solarlong(float(day_ini),zls)
837	end if
838
839	c Initialize pressure levels
840	c ~~~~~~~~~~~~~~~~~~
841	zplev(:,:) = pplev(:,:)
842	zplay(:,:) = pplay(:,:)
843	ps(:) = pplev(:,1)
844
845	c Compute geopotential at interlayers
846	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
847	c ponderation des altitudes au niveau des couches en dp/p
848
849	DO l=1,nlayer
850	DO ig=1,ngrid
851	zzlay(ig,l)=pphi(ig,l)/g
852	ENDDO
853	ENDDO
854	DO ig=1,ngrid
855	zzlev(ig,1)=0.
856	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
857	ENDDO
858	DO l=2,nlayer
859	DO ig=1,ngrid
860	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
861	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
862	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
863	ENDDO
864	ENDDO
865
866
867	! Potential temperature calculation not the same in physiq and dynamic
868
869	c Compute potential temperature
870	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
871	DO l=1,nlayer
872	DO ig=1,ngrid
873	zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp
874	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
875	ENDDO
876	ENDDO
877
878	#ifndef MESOSCALE
879	c-----------------------------------------------------------------------
880	c 1.2.5 Compute mean mass, cp, and R
881	c --------------------------------
882
883	if(photochem.or.callthermos) then
884	call concentrations(ngrid,nlayer,nq,
885	& zplay,pt,pdt,pq,pdq,ptimestep)
886	endif
887	#endif
888
889	! Compute vertical velocity (m/s) from vertical mass flux
890	! w = F / (rhoarea) and rho = P/(rT)
891	! but first linearly interpolate mass flux to mid-layers
892	do l=1,nlayer-1
893	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
894	enddo
895	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
896	do l=1,nlayer
897	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) /
898	& (pplay(1:ngrid,l)*cell_area(1:ngrid))
899	! NB: here we use r and not rnew since this diagnostic comes
900	! from the dynamics
901	enddo
902
903	! test for co2 conservation with co2 microphysics
904	if (igcm_co2_ice.ne.0) then
905	! calculates the amount of co2 at the beginning of physics
906	co2totA = 0.
907	do ig=1,ngrid
908	do l=1,nlayer
909	co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g*
910	& (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice)
911	& +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
912	end do
913	do islope = 1,nslope
914	co2totA = co2totA + qsurf(ig,igcm_co2,islope)*
915	& subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
916	enddo
917	end do
918	else
919	co2totA = 0.
920	do ig=1,ngrid
921	do l=1,nlayer
922	co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g*
923	& (pq(ig,l,igcm_co2)
924	& +pdq(ig,l,igcm_co2)*ptimestep)
925	end do
926	do islope = 1,nslope
927	co2totA = co2totA + qsurf(ig,igcm_co2,islope)*
928	& subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
929	enddo
930	end do
931	endif ! of if (igcm_co2_ice.ne.0)
932	c-----------------------------------------------------------------------
933	c 2. Compute radiative tendencies :
934	c------------------------------------
935
936	IF (callrad) THEN
937
938	c Local Solar zenith angle
939	c ~~~~~~~~~~~~~~~~~~~~~~~~
940
941	CALL orbite(zls,dist_sol,declin)
942
943	IF (diurnal) THEN
944	ztim1=SIN(declin)
945	ztim2=COS(declin)COS(2.pi*(zday-.5))
946	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
947
948	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
949	& ztim1,ztim2,ztim3, mu0,fract)
950
951	ELSE
952	CALL mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad)
953	ENDIF ! of IF (diurnal)
954
955	IF( MOD(icount-1,iradia).EQ.0) THEN
956
957	c NLTE cooling from CO2 emission
958	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
959	IF(callnlte) then
960	if(nltemodel.eq.0.or.nltemodel.eq.1) then
961	CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte)
962	else if(nltemodel.eq.2) then
963	co2vmr_gcm(1:ngrid,1:nlayer)=
964	& pq(1:ngrid,1:nlayer,igcm_co2)*
965	& mmean(1:ngrid,1:nlayer)/mmol(igcm_co2)
966	n2vmr_gcm(1:ngrid,1:nlayer)=
967	& pq(1:ngrid,1:nlayer,igcm_n2)*
968	& mmean(1:ngrid,1:nlayer)/mmol(igcm_n2)
969	covmr_gcm(1:ngrid,1:nlayer)=
970	& pq(1:ngrid,1:nlayer,igcm_co)*
971	& mmean(1:ngrid,1:nlayer)/mmol(igcm_co)
972	ovmr_gcm(1:ngrid,1:nlayer)=
973	& pq(1:ngrid,1:nlayer,igcm_o)*
974	& mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
975
976	CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
977	$ pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
978	$ ovmr_gcm, zdtnlte,ierr_nlte,varerr )
979	if(ierr_nlte.gt.0) then
980	write(,)
981	$ 'WARNING: nlte_tcool output with error message',
982	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
983	write(,)'I will continue anyway'
984	endif
985
986	zdtnlte(1:ngrid,1:nlayer)=
987	& zdtnlte(1:ngrid,1:nlayer)/86400.
988	endif
989	ELSE
990	zdtnlte(:,:)=0.
991	ENDIF !end callnlte
992
993	c Find number of layers for LTE radiation calculations
994	IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
995	& CALL nlthermeq(ngrid,nlayer,zplev,zplay)
996
997	c rocketstorm : compute dust storm mesh fraction
998	IF (rdstorm) THEN
999	CALL calcstormfract(ngrid,nlayer,nq,pq,
1000	& totstormfract)
1001	ENDIF
1002
1003	c Note: Dustopacity.F has been transferred to callradite.F
1004
1005	#ifdef DUSTSTORM
1006	!! specific case: save the quantity of dust before adding perturbation
1007
1008	if (firstcall) then
1009	pq_tmp(1:ngrid,1:nlayer,1)=pq(1:ngrid,1:nlayer,igcm_dust_mass)
1010	pq_tmp(1:ngrid,1:nlayer,2)=pq(1:ngrid,1:nlayer,igcm_dust_number)
1011	endif
1012	#endif
1013
1014	c Call main radiative transfer scheme
1015	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1016	c Transfer through CO2 (except NIR CO2 absorption)
1017	c and aerosols (dust and water ice)
1018	! callradite for background dust (out of the rdstorm fraction)
1019	clearatm=.true.
1020	!! callradite for background dust (out of the topflows fraction)
1021	nohmons=.true.
1022
1023	c Radiative transfer
1024	c ------------------
1025	! callradite for the part with clouds
1026	clearsky=.false. ! part with clouds for both cases CLFvarying true/false
1027	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
1028	& albedo_meshavg,emis_meshavg,
1029	& mu0,zplev,zplay,pt,tsurf_meshavg,fract,dist_sol,igout,
1030	& zdtlw,zdtsw,fluxsurf_lw(:,iflat),fluxsurf_dn_sw(:,:,iflat),
1031	& fluxsurf_up_sw,
1032	& fluxtop_lw,fluxtop_dn_sw,fluxtop_up_sw,
1033	& tau_pref_scenario,tau_pref_gcm,
1034	& tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
1035	& taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,
1036	& qsurf_meshavg(:,igcm_co2),rstormdust,rtopdust,totstormfract,
1037	& clearatm,dsords,dsotop,nohmons,clearsky,totcloudfrac)
1038
1039	DO islope=1,nslope
1040	fluxsurf_lw(:,islope) =fluxsurf_lw(:,iflat)
1041	fluxsurf_dn_sw(:,:,islope) =fluxsurf_dn_sw(:,:,iflat)
1042	ENDDO
1043
1044	! case of sub-grid water ice clouds: callradite for the clear case
1045	IF (CLFvarying) THEN
1046	! ---> PROBLEMS WITH ALLOCATED ARRAYS
1047	! (temporary solution in callcorrk: do not deallocate
1048	! if
1049	! CLFvarying ...) ?? AP ??
1050	clearsky=.true.
1051	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
1052	& albedo_meshavg,emis_meshavg,mu0,zplev,zplay,pt,
1053	& tsurf_meshavg,fract,
1054	& dist_sol,igout,zdtlwclf,zdtswclf,
1055	& fluxsurf_lwclf,fluxsurf_dn_swclf,fluxsurf_up_swclf,
1056	& fluxtop_lwclf,fluxtop_dn_swclf,fluxtop_up_swclf,
1057	& tau_pref_scenario,tau_pref_gcm,tau,aerosol,
1058	& dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
1059	& taucloudtesclf,rdust,
1060	& rice,nuice,riceco2, nuiceco2,
1061	& qsurf_meshavg(:,igcm_co2),
1062	& rstormdust,rtopdust,totstormfract,
1063	& clearatm,dsords,dsotop,
1064	& nohmons,clearsky,totcloudfrac)
1065	clearsky = .false. ! just in case.
1066	! Sum the fluxes and heating rates from cloudy/clear
1067	! cases
1068	DO ig=1,ngrid
1069	tf_clf=totcloudfrac(ig)
1070	ntf_clf=1.-tf_clf
1071	DO islope=1,nslope
1072	fluxsurf_lw(ig,islope) = ntf_clf*fluxsurf_lwclf(ig)
1073	& + tf_clf*fluxsurf_lw(ig,islope)
1074	fluxsurf_dn_sw(ig,1:2,islope) =
1075	& ntf_clf*fluxsurf_dn_swclf(ig,1:2)
1076	& + tf_clf*fluxsurf_dn_sw(ig,1:2,islope)
1077	ENDDO
1078	fluxsurf_up_sw(ig,1:2) =
1079	& ntf_clf*fluxsurf_up_swclf(ig,1:2)
1080	& + tf_clf*fluxsurf_up_sw(ig,1:2)
1081	fluxtop_lw(ig) = ntf_clf*fluxtop_lwclf(ig)
1082	& + tf_clf*fluxtop_lw(ig)
1083	fluxtop_dn_sw(ig,1:2)=ntf_clf*fluxtop_dn_swclf(ig,1:2)
1084	& + tf_clf*fluxtop_dn_sw(ig,1:2)
1085	fluxtop_up_sw(ig,1:2)=ntf_clf*fluxtop_up_swclf(ig,1:2)
1086	& + tf_clf*fluxtop_up_sw(ig,1:2)
1087	taucloudtes(ig) = ntf_clf*taucloudtesclf(ig)
1088	& + tf_clf*taucloudtes(ig)
1089	zdtlw(ig,1:nlayer) = ntf_clf*zdtlwclf(ig,1:nlayer)
1090	& + tf_clf*zdtlw(ig,1:nlayer)
1091	zdtsw(ig,1:nlayer) = ntf_clf*zdtswclf(ig,1:nlayer)
1092	& + tf_clf*zdtsw(ig,1:nlayer)
1093	ENDDO
1094
1095	ENDIF ! (CLFvarying)
1096
1097	!============================================================================
1098
1099	#ifdef DUSTSTORM
1100	!! specific case: compute the added quantity of dust for perturbation
1101	if (firstcall) then
1102	pdq(1:ngrid,1:nlayer,igcm_dust_mass)=
1103	& pdq(1:ngrid,1:nlayer,igcm_dust_mass)
1104	& - pq_tmp(1:ngrid,1:nlayer,1)
1105	& + pq(1:ngrid,1:nlayer,igcm_dust_mass)
1106	pdq(1:ngrid,1:nlayer,igcm_dust_number)=
1107	& pdq(1:ngrid,1:nlayer,igcm_dust_number)
1108	& - pq_tmp(1:ngrid,1:nlayer,2)
1109	& + pq(1:ngrid,1:nlayer,igcm_dust_number)
1110	endif
1111	#endif
1112
1113	c Outputs for basic check (middle of domain)
1114	c ------------------------------------------
1115	write(*,'("Ls =",f11.6," check lat =",f10.6,
1116	& " lon =",f11.6)')
1117	& zls180./pi,latitude(igout)180/pi,
1118	& longitude(igout)*180/pi
1119
1120	write(*,'(" tau_pref_gcm(",f4.0," Pa) =",f9.6,
1121	& " tau(",f4.0," Pa) =",f9.6)')
1122	& odpref,tau_pref_gcm(igout),
1123	& odpref,tau(igout,1)*odpref/zplev(igout,1)
1124
1125
1126	c ---------------------------------------------------------
1127	c Call slope parameterization for direct and scattered flux
1128	c ---------------------------------------------------------
1129	IF(callslope) THEN
1130	! assume that in this case, nslope = 1
1131	if(nslope.ne.1) then
1132	call abort_physic(
1133	& "physiq","callslope=true but nslope.ne.1",1)
1134	endif
1135	print *, 'Slope scheme is on and computing...'
1136	DO ig=1,ngrid
1137	sl_the = theta_sl(ig)
1138	IF (sl_the .ne. 0.) THEN
1139	ztim1=fluxsurf_dn_sw(ig,1,iflat)
1140	& +fluxsurf_dn_sw(ig,2,iflat)
1141	DO l=1,2
1142	sl_lct = ptime24. + 180.longitude(ig)/pi/15.
1143	sl_ra = pi*(1.0-sl_lct/12.)
1144	sl_lat = 180.*latitude(ig)/pi
1145	sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
1146	sl_alb = albedo(ig,l,iflat)
1147	sl_psi = psi_sl(ig)
1148	sl_fl0 = fluxsurf_dn_sw(ig,l,iflat)
1149	sl_di0 = 0.
1150	if ((mu0(ig) .gt. 0.).and.(ztim1.gt.0.)) then
1151	sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
1152	sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
1153	sl_di0 = sl_di0/ztim1
1154	sl_di0 = fluxsurf_dn_sw(ig,l,iflat)*sl_di0
1155	endif
1156	! you never know (roundup concern...)
1157	if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
1158	!!!!!!!!!!!!!!!!!!!!!!!!!!
1159	CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
1160	& sl_tau, sl_alb, sl_the, sl_psi,
1161	& sl_di0, sl_fl0, sl_flu )
1162	!!!!!!!!!!!!!!!!!!!!!!!!!!
1163	fluxsurf_dn_sw(ig,l,1) = sl_flu
1164	ENDDO
1165	!!! compute correction on IR flux as well
1166	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
1167	fluxsurf_lw(ig,:)= fluxsurf_lw(ig,:)*sky
1168	ENDIF
1169	ENDDO
1170	ELSE ! not calling subslope, nslope might be > 1
1171	DO islope = 1,nslope
1172	IF(def_slope_mean(islope).ge.0.) THEN
1173	sl_psi = 0. !Northward slope
1174	ELSE
1175	sl_psi = 180. !Southward slope
1176	ENDIF
1177	sl_the=abs(def_slope_mean(islope))
1178	IF (sl_the .gt. 1e-6) THEN
1179	DO ig=1,ngrid
1180	ztim1=fluxsurf_dn_sw(ig,1,islope)
1181	s +fluxsurf_dn_sw(ig,2,islope)
1182	DO l=1,2
1183	sl_lct = ptime24. + 180.longitude(ig)/pi/15.
1184	sl_ra = pi*(1.0-sl_lct/12.)
1185	sl_lat = 180.*latitude(ig)/pi
1186	sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
1187	sl_alb = albedo(ig,l,islope)
1188	sl_psi = psi_sl(ig)
1189	sl_fl0 = fluxsurf_dn_sw(ig,l,islope)
1190	sl_di0 = 0.
1191	if (mu0(ig) .gt. 0.) then
1192	sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
1193	sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
1194	sl_di0 = sl_di0/ztim1
1195	sl_di0 = fluxsurf_dn_sw(ig,l,islope)*sl_di0
1196	endif
1197	! you never know (roundup concern...)
1198	if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
1199	!!!!!!!!!!!!!!!!!!!!!!!!!!
1200	CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
1201	& sl_tau, sl_alb, sl_the, sl_psi,
1202	& sl_di0, sl_fl0, sl_flu )
1203	!!!!!!!!!!!!!!!!!!!!!!!!!!
1204	fluxsurf_dn_sw(ig,l,islope) = sl_flu
1205	ENDDO
1206	!!! compute correction on IR flux as well
1207
1208	fluxsurf_lw(ig,islope)= fluxsurf_lw(ig,islope)
1209	& *sky_slope(islope)
1210	ENDDO
1211	ENDIF ! sub grid is not flat
1212	ENDDO ! islope = 1,nslope
1213	ENDIF ! callslope
1214
1215	c CO2 near infrared absorption
1216	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1217	zdtnirco2(:,:)=0
1218	if (callnirco2) then
1219	call nirco2abs (ngrid,nlayer,zplay,dist_sol,nq,pq,
1220	. mu0,fract,declin, zdtnirco2)
1221	endif
1222
1223	c Radiative flux from the sky absorbed by the surface (W.m-2)
1224	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1225	DO ig=1,ngrid
1226	DO islope = 1,nslope
1227	fluxrad_sky(ig,islope) =
1228	$ emis(ig,islope)*fluxsurf_lw(ig,islope)
1229	$ +fluxsurf_dn_sw(ig,1,islope)*(1.-albedo(ig,1,islope))
1230	$ +fluxsurf_dn_sw(ig,2,islope)*(1.-albedo(ig,2,islope))
1231	ENDDO
1232	ENDDO
1233
1234
1235	c Net atmospheric radiative heating rate (K.s-1)
1236	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1237	IF(callnlte) THEN
1238	CALL blendrad(ngrid, nlayer, zplay,
1239	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
1240	ELSE
1241	DO l=1,nlayer
1242	DO ig=1,ngrid
1243	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
1244	& +zdtnirco2(ig,l)
1245	ENDDO
1246	ENDDO
1247	ENDIF
1248
1249	ENDIF ! of if(mod(icount-1,iradia).eq.0)
1250
1251	c Transformation of the radiative tendencies:
1252	c -------------------------------------------
1253
1254	c Net radiative surface flux (W.m-2)
1255	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
1256	c
1257	DO ig=1,ngrid
1258	DO islope = 1,nslope
1259	zplanck(ig)=tsurf(ig,islope)*tsurf(ig,islope)
1260	zplanck(ig)=emis(ig,islope)*
1261	$ stephanzplanck(ig)zplanck(ig)
1262	fluxrad(ig,islope)=fluxrad_sky(ig,islope)-zplanck(ig)
1263	IF(callslope) THEN
1264	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
1265	fluxrad(ig,nslope)=fluxrad(ig,nslope)+
1266	$ (1.-sky)*zplanck(ig)
1267	ELSE
1268	fluxrad(ig,islope)=fluxrad(ig,islope) +
1269	$ (1.-sky_slope(iflat))emis(ig,iflat)
1270	$ stephantsurf(ig,iflat)*4
1271	ENDIF
1272	ENDDO
1273	ENDDO
1274
1275	DO l=1,nlayer
1276	DO ig=1,ngrid
1277	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
1278	ENDDO
1279	ENDDO
1280
1281	ENDIF ! of IF (callrad)
1282
1283	c 3.1 Rocket dust storm
1284	c -------------------------------------------
1285	IF (rdstorm) THEN
1286	clearatm=.false.
1287	pdqrds(:,:,:)=0.
1288	qdusttotal0(:)=0.
1289	qdusttotal1(:)=0.
1290	do ig=1,ngrid
1291	do l=1,nlayer
1292	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) ! updated potential
1293	! temperature tendency
1294	! for diagnostics
1295	! qdustrds0(ig,l)=pq(ig,l,igcm_dust_mass)+
1296	! & pdq(ig,l,igcm_dust_mass)*ptimestep
1297	! qstormrds0(ig,l)=pq(ig,l,igcm_stormdust_mass)+
1298	! & pdq(ig,l,igcm_stormdust_mass)*ptimestep
1299	! qdusttotal0(ig)=qdusttotal0(ig)+(qdustrds0(ig,l)+
1300	! & qstormrds0(ig,l))*(zplev(ig,l)-
1301	! & zplev(ig,l+1))/g
1302	enddo
1303	enddo
1304	! call write_output('qdustrds0','qdust before rds',
1305	! & 'kg/kg ',qdustrds0(:,:))
1306	! call write_output('qstormrds0','qstorm before rds',
1307	! & 'kg/kg ',qstormrds0(:,:))
1308
1309	CALL rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep,
1310	& pq,pdq,pt,pdt,zplev,zplay,zzlev,
1311	& zzlay,zdtsw,zdtlw,
1312	c for radiative transfer
1313	& clearatm,icount,zday,zls,
1314	& tsurf,qsurf_meshavg(:,igcm_co2),igout,
1315	& totstormfract,tauscaling,
1316	& dust_rad_adjust,IRtoVIScoef,
1317	c input sub-grid scale cloud
1318	& clearsky,totcloudfrac,
1319	c input sub-grid scale topography
1320	& nohmons,
1321	c output
1322	& pdqrds,wspeed,dsodust,dsords,dsotop,
1323	& tau_pref_scenario,tau_pref_gcm)
1324
1325	c update the tendencies of both dust after vertical transport
1326	DO l=1,nlayer
1327	DO ig=1,ngrid
1328	pdq(ig,l,igcm_stormdust_mass)=
1329	& pdq(ig,l,igcm_stormdust_mass)+
1330	& pdqrds(ig,l,igcm_stormdust_mass)
1331	pdq(ig,l,igcm_stormdust_number)=
1332	& pdq(ig,l,igcm_stormdust_number)+
1333	& pdqrds(ig,l,igcm_stormdust_number)
1334
1335	pdq(ig,l,igcm_dust_mass)=
1336	& pdq(ig,l,igcm_dust_mass)+ pdqrds(ig,l,igcm_dust_mass)
1337	pdq(ig,l,igcm_dust_number)=
1338	& pdq(ig,l,igcm_dust_number)+
1339	& pdqrds(ig,l,igcm_dust_number)
1340
1341	ENDDO
1342	ENDDO
1343	do l=1,nlayer
1344	do ig=1,ngrid
1345	qdustrds1(ig,l)=pq(ig,l,igcm_dust_mass)+
1346	& pdq(ig,l,igcm_dust_mass)*ptimestep
1347	qstormrds1(ig,l)=pq(ig,l,igcm_stormdust_mass)+
1348	& pdq(ig,l,igcm_stormdust_mass)*ptimestep
1349	qdusttotal1(ig)=qdusttotal1(ig)+(qdustrds1(ig,l)+
1350	& qstormrds1(ig,l))*(zplev(ig,l)-
1351	& zplev(ig,l+1))/g
1352	enddo
1353	enddo
1354
1355	c for diagnostics
1356	! call write_output('qdustrds1','qdust after rds',
1357	! & 'kg/kg ',qdustrds1(:,:))
1358	! call write_output('qstormrds1','qstorm after rds',
1359	! & 'kg/kg ',qstormrds1(:,:))
1360	!
1361	! call write_output('qdusttotal0','q sum before rds',
1362	! & 'kg/m2 ',qdusttotal0(:))
1363	! call write_output('qdusttotal1','q sum after rds',
1364	! & 'kg/m2 ',qdusttotal1(:))
1365
1366	ENDIF ! end of if(rdstorm)
1367
1368	c 3.2 Dust entrained from the PBL up to the top of sub-grid scale topography
1369	c -------------------------------------------
1370	IF (topflows) THEN
1371	clearatm=.true. ! stormdust is not accounted in the extra heating on top of the mountains
1372	nohmons=.false.
1373	pdqtop(:,:,:)=0.
1374	CALL topmons(ngrid,nlayer,nq,ptime,ptimestep,
1375	& pq,pdq,pt,pdt,zplev,zplay,zzlev,
1376	& zzlay,zdtsw,zdtlw,
1377	& icount,zday,zls,tsurf(:,iflat),
1378	& qsurf_meshavg(:,igcm_co2),
1379	& igout,aerosol,tauscaling,dust_rad_adjust,
1380	& IRtoVIScoef,totstormfract,clearatm,
1381	& clearsky,totcloudfrac,
1382	& nohmons,
1383	& pdqtop,wtop,dsodust,dsords,dsotop,
1384	& tau_pref_scenario,tau_pref_gcm)
1385
1386	c update the tendencies of both dust after vertical transport
1387	DO l=1,nlayer
1388	DO ig=1,ngrid
1389	pdq(ig,l,igcm_topdust_mass)=
1390	& pdq(ig,l,igcm_topdust_mass)+
1391	& pdqtop(ig,l,igcm_topdust_mass)
1392	pdq(ig,l,igcm_topdust_number)=
1393	& pdq(ig,l,igcm_topdust_number)+
1394	& pdqtop(ig,l,igcm_topdust_number)
1395	pdq(ig,l,igcm_dust_mass)=
1396	& pdq(ig,l,igcm_dust_mass)+ pdqtop(ig,l,igcm_dust_mass)
1397	pdq(ig,l,igcm_dust_number)=
1398	& pdq(ig,l,igcm_dust_number)+pdqtop(ig,l,igcm_dust_number)
1399
1400	ENDDO
1401	ENDDO
1402
1403	ENDIF ! end of if (topflows)
1404
1405	c 3.3 Dust injection from the surface
1406	c -------------------------------------------
1407	if (dustinjection.gt.0) then
1408
1409	CALL compute_dtau(ngrid,nlayer,
1410	& zday,pplev,tau_pref_gcm,
1411	& ptimestep,local_time,IRtoVIScoef,
1412	& dustliftday)
1413	endif ! end of if (dustinjection.gt.0)
1414
1415
1416	c-----------------------------------------------------------------------
1417	c 4. Gravity wave and subgrid scale topography drag :
1418	c -------------------------------------------------
1419
1420
1421	IF(calllott)THEN
1422	CALL calldrag_noro(ngrid,nlayer,ptimestep,
1423	& zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
1424
1425	DO l=1,nlayer
1426	DO ig=1,ngrid
1427	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
1428	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
1429	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
1430	ENDDO
1431	ENDDO
1432	ENDIF
1433
1434	c-----------------------------------------------------------------------
1435	c 5. Vertical diffusion (turbulent mixing):
1436	c -----------------------------------------
1437
1438	IF (calldifv) THEN
1439	DO ig=1,ngrid
1440	DO islope = 1,nslope
1441	zflubid(ig,islope)=fluxrad(ig,islope)
1442	& +fluxgrd(ig,islope)
1443	ENDDO
1444	ENDDO
1445	zdum1(:,:)=0
1446	zdum2(:,:)=0
1447	do l=1,nlayer
1448	do ig=1,ngrid
1449	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
1450	enddo
1451	enddo
1452
1453	c ----------------------
1454	c Treatment of a special case : using new surface layer (Richardson based)
1455	c without using the thermals in gcm and mesoscale can yield problems in
1456	c weakly unstable situations when winds are near to 0. For those cases, we add
1457	c a unit subgrid gustiness. Remember that thermals should be used we using the
1458	c Richardson based surface layer model.
1459	IF ( .not.calltherm
1460	. .and. callrichsl
1461	. .and. .not.turb_resolved) THEN
1462
1463	DO ig=1, ngrid
1464	IF (zh(ig,1) .lt. tsurf_meshavg(ig)) THEN
1465	wstar(ig)=1.
1466	hfmax_th(ig)=0.2
1467	ELSE
1468	wstar(ig)=0.
1469	hfmax_th(ig)=0.
1470	ENDIF
1471	ENDDO
1472	ENDIF
1473	c ----------------------
1474
1475	IF (tke_heat_flux .ne. 0.) THEN
1476
1477	zz1(:)=(pt(:,1)+pdt(:,1)ptimestep)(r/g)*
1478	& (-alog(zplay(:,1)/zplev(:,1)))
1479	pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
1480	ENDIF
1481
1482	c Calling vdif (Martian version WITH CO2 condensation)
1483	dwatercap_dif(:,:) = 0.
1484	zcdh(:) = 0.
1485	zcdv(:) = 0.
1486
1487	CALL vdifc(ngrid,nlayer,nq,zpopsk,
1488	$ ptimestep,capcal,lwrite,
1489	$ zplay,zplev,zzlay,zzlev,z0,
1490	$ pu,pv,zh,pq,tsurf,emis,qsurf,
1491	$ zdum1,zdum2,zdh,pdq,zflubid,
1492	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
1493	& zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,
1494	& zcondicea_co2microp,sensibFlux,
1495	& dustliftday,local_time,watercap,dwatercap_dif)
1496
1497	DO ig=1,ngrid
1498	zdtsurf(ig,:)=zdtsurf(ig,:)+zdtsdif(ig,:)
1499	dwatercap(ig,:)=dwatercap(ig,:)+dwatercap_dif(ig,:)
1500	ENDDO
1501	call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,zdqsdif,
1502	& albedo_meshavg,emis_meshavg,tsurf_meshavg,zdqsdif_meshavg_tmp)
1503	IF (.not.turb_resolved) THEN
1504	DO l=1,nlayer
1505	DO ig=1,ngrid
1506	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
1507	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
1508	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
1509
1510	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
1511	ENDDO
1512	ENDDO
1513
1514	DO iq=1, nq
1515	DO l=1,nlayer
1516	DO ig=1,ngrid
1517	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
1518	ENDDO
1519	ENDDO
1520	ENDDO
1521	DO iq=1, nq
1522	DO ig=1,ngrid
1523	dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:)
1524	ENDDO
1525	ENDDO
1526
1527	ELSE
1528	write (,) '******************************************'
1529	write (,) '** LES mode: the difv part is only used to'
1530	write (,) '** - provide HFX and UST to the dynamics'
1531	write (,) '** - update TSURF'
1532	write (,) '******************************************'
1533	!! Specific treatment for lifting in turbulent-resolving mode (AC)
1534	IF (lifting .and. doubleq) THEN
1535	!! lifted dust is injected in the first layer.
1536	!! Sedimentation must be called after turbulent mixing, i.e. on next step, after WRF.
1537	!! => lifted dust is not incremented before the sedimentation step.
1538	zdqdif(1:ngrid,1,1:nq)=0.
1539	zdqdif(1:ngrid,1,igcm_dust_number) =
1540	. -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_number)
1541	zdqdif(1:ngrid,1,igcm_dust_mass) =
1542	. -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_mass)
1543	zdqdif(1:ngrid,2:nlayer,1:nq) = 0.
1544	DO iq=1, nq
1545	IF ((iq .ne. igcm_dust_mass)
1546	& .and. (iq .ne. igcm_dust_number)) THEN
1547	zdqsdif(:,iq,:)=0.
1548	ENDIF
1549	ENDDO
1550	ELSE
1551	zdqdif(1:ngrid,1:nlayer,1:nq) = 0.
1552	zdqsdif(1:ngrid,1:nq,1:nslope) = 0.
1553	ENDIF
1554	ENDIF
1555	ELSE
1556	DO ig=1,ngrid
1557	DO islope=1,nslope
1558	zdtsurf(ig,islope)=zdtsurf(ig,islope)+
1559	s (fluxrad(ig,islope)+fluxgrd(ig,islope))/capcal(ig,islope)
1560	ENDDO
1561	ENDDO
1562	IF (turb_resolved) THEN
1563	write(,) 'Turbulent-resolving mode !'
1564	write(,) 'Please set calldifv to T in callphys.def'
1565	call abort_physic("physiq","turbulent-resolving mode",1)
1566	ENDIF
1567	ENDIF ! of IF (calldifv)
1568
1569	c-----------------------------------------------------------------------
1570	c 6. Thermals :
1571	c -----------------------------
1572
1573	if(calltherm .and. .not.turb_resolved) then
1574
1575	call calltherm_interface(ngrid,nlayer,nq,igcm_co2,
1576	$ zzlev,zzlay,
1577	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
1578	$ zplay,zplev,pphi,zpopsk,
1579	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
1580	$ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
1581
1582	DO l=1,nlayer
1583	DO ig=1,ngrid
1584	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
1585	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
1586	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
1587	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
1588	ENDDO
1589	ENDDO
1590
1591	DO ig=1,ngrid
1592	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
1593	ENDDO
1594
1595	DO iq=1,nq
1596	DO l=1,nlayer
1597	DO ig=1,ngrid
1598	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
1599	ENDDO
1600	ENDDO
1601	ENDDO
1602
1603	lmax_th_out(:)=real(lmax_th(:))
1604
1605	else !of if calltherm
1606	lmax_th(:)=0
1607	wstar(:)=0.
1608	hfmax_th(:)=0.
1609	lmax_th_out(:)=0.
1610	end if
1611
1612	c-----------------------------------------------------------------------
1613	c 7. Dry convective adjustment:
1614	c -----------------------------
1615
1616	IF(calladj) THEN
1617
1618	DO l=1,nlayer
1619	DO ig=1,ngrid
1620	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
1621	ENDDO
1622	ENDDO
1623	zduadj(:,:)=0
1624	zdvadj(:,:)=0
1625	zdhadj(:,:)=0
1626
1627	CALL convadj(ngrid,nlayer,nq,ptimestep,
1628	$ zplay,zplev,zpopsk,lmax_th,
1629	$ pu,pv,zh,pq,
1630	$ pdu,pdv,zdh,pdq,
1631	$ zduadj,zdvadj,zdhadj,
1632	$ zdqadj)
1633
1634	DO l=1,nlayer
1635	DO ig=1,ngrid
1636	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
1637	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
1638	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
1639
1640	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
1641	ENDDO
1642	ENDDO
1643
1644	DO iq=1, nq
1645	DO l=1,nlayer
1646	DO ig=1,ngrid
1647	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
1648	ENDDO
1649	ENDDO
1650	ENDDO
1651	ENDIF ! of IF(calladj)
1652
1653	c-----------------------------------------------------
1654	c 8. Non orographic Gravity waves :
1655	c -------------------------------------------------
1656
1657	IF (calllott_nonoro) THEN
1658
1659	CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep,
1660	& cpnew,rnew,
1661	& zplay,
1662	& zmax_th, ! max altitude reached by thermals (m)
1663	& pt, pu, pv,
1664	& pdt, pdu, pdv,
1665	& zustrhi,zvstrhi,
1666	& d_t_hin, d_u_hin, d_v_hin)
1667
1668	! Update tendencies
1669	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)
1670	& +d_t_hin(1:ngrid,1:nlayer)
1671	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)
1672	& +d_u_hin(1:ngrid,1:nlayer)
1673	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)
1674	& +d_v_hin(1:ngrid,1:nlayer)
1675
1676	ENDIF ! of IF (calllott_nonoro)
1677
1678	c-----------------------------------------------------------------------
1679	c 9. Specific parameterizations for tracers
1680	c: -----------------------------------------
1681
1682
1683	c 9a. Water and ice
1684	c ---------------
1685
1686	c ---------------------------------------
1687	c Water ice condensation in the atmosphere
1688	c ----------------------------------------
1689	IF (water) THEN
1690
1691	call watercloud(ngrid,nlayer,ptimestep,
1692	& zplev,zplay,pdpsrf,zzlay, pt,pdt,
1693	& pq,pdq,zdqcloud,zdtcloud,
1694	& nq,tau,tauscaling,rdust,rice,nuice,
1695	& rsedcloud,rhocloud,totcloudfrac)
1696	c Temperature variation due to latent heat release
1697	if (activice) then
1698	pdt(1:ngrid,1:nlayer) =
1699	& pdt(1:ngrid,1:nlayer) +
1700	& zdtcloud(1:ngrid,1:nlayer)
1701	endif
1702
1703	! increment water vapour and ice atmospheric tracers tendencies
1704	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) =
1705	& pdq(1:ngrid,1:nlayer,igcm_h2o_vap) +
1706	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_vap)
1707	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
1708	& pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
1709	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice)
1710
1711	if (hdo) then
1712	! increment HDO vapour and ice atmospheric tracers tendencies
1713	pdq(1:ngrid,1:nlayer,igcm_hdo_vap) =
1714	& pdq(1:ngrid,1:nlayer,igcm_hdo_vap) +
1715	& zdqcloud(1:ngrid,1:nlayer,igcm_hdo_vap)
1716	pdq(1:ngrid,1:nlayer,igcm_hdo_ice) =
1717	& pdq(1:ngrid,1:nlayer,igcm_hdo_ice) +
1718	& zdqcloud(1:ngrid,1:nlayer,igcm_hdo_ice)
1719	endif !hdo
1720
1721	! increment dust and ccn masses and numbers
1722	! We need to check that we have Nccn & Ndust > 0
1723	! This is due to single precision rounding problems
1724	if (microphys) then
1725	pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
1726	& pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
1727	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
1728	pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
1729	& pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
1730	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
1731	where (pq(:,:,igcm_ccn_mass) +
1732	& ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
1733	pdq(:,:,igcm_ccn_mass) =
1734	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1735	pdq(:,:,igcm_ccn_number) =
1736	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1737	end where
1738	where (pq(:,:,igcm_ccn_number) +
1739	& ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
1740	pdq(:,:,igcm_ccn_mass) =
1741	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1742	pdq(:,:,igcm_ccn_number) =
1743	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1744	end where
1745	endif
1746
1747	if (scavenging) then
1748	pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
1749	& pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
1750	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_mass)
1751	pdq(1:ngrid,1:nlayer,igcm_dust_number) =
1752	& pdq(1:ngrid,1:nlayer,igcm_dust_number) +
1753	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_number)
1754	where (pq(:,:,igcm_dust_mass) +
1755	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
1756	pdq(:,:,igcm_dust_mass) =
1757	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1758	pdq(:,:,igcm_dust_number) =
1759	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1760	end where
1761	where (pq(:,:,igcm_dust_number) +
1762	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
1763	pdq(:,:,igcm_dust_mass) =
1764	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1765	pdq(:,:,igcm_dust_number) =
1766	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1767	end where
1768	endif ! of if scavenging
1769
1770	END IF ! of IF (water)
1771
1772	c 9a bis. CO2 clouds (CL & JA)
1773	c ---------------------------------------
1774	c CO2 ice cloud condensation in the atmosphere
1775	c ----------------------------------------
1776	c flag needed in callphys.def:
1777	c co2clouds=.true. is mandatory (default is .false.)
1778	c co2useh2o=.true. if you want to allow co2 condensation
1779	c on water ice particles
1780	c meteo_flux=.true. if you want to add a meteoritic
1781	c supply of CCN
1782	c CLFvaryingCO2=.true. if you want to have a sub-grid
1783	c temperature distribution
1784	c spantCO2=integer (i.e. 3) amplitude of the sub-grid T disti
1785	c nuiceco2_sed=0.2 variance of the size distribution for the
1786	c sedimentation
1787	c nuiceco2_ref=0.2 variance of the size distribution for the
1788	c nucleation
1789	c imicroco2=50 micro-timestep is 1/50 of physical timestep
1790	zdqssed_co2(:) = 0.
1791	zdqssed_ccn(:,:) = 0.
1792
1793	IF (co2clouds) THEN
1794	call co2cloud(ngrid,nlayer,ptimestep,
1795	& zplev,zplay,pdpsrf,zzlay,pt,pdt,
1796	& pq,pdq,zdqcloudco2,zdtcloudco2,
1797	& nq,tau,tauscaling,rdust,rice,riceco2,nuice,
1798	& rhocloud, rsedcloudco2,rhocloudco2,zzlev,zdqssed_co2,
1799	& zdqssed_ccn,pdu,pu,zcondicea_co2microp)
1800
1801	DO iq=1, nq
1802	DO ig=1,ngrid
1803	DO islope = 1,nslope
1804	dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope)+
1805	& zdqssed_ccn(ig,iq)cos(pidef_slope_mean(islope)/180.)
1806	ENDDO !(islope)
1807	ENDDO ! (ig)
1808	ENDDO ! (iq)q)
1809	c Temperature variation due to latent heat release
1810	pdt(1:ngrid,1:nlayer) =
1811	& pdt(1:ngrid,1:nlayer) +
1812	& zdtcloudco2(1:ngrid,1:nlayer)
1813
1814	! increment dust and ccn masses and numbers
1815	! We need to check that we have Nccn & Ndust > 0
1816	! This is due to single precision rounding problems
1817	! increment dust tracers tendancies
1818	pdq(:,:,igcm_dust_mass) = pdq(:,:,igcm_dust_mass)
1819	& + zdqcloudco2(:,:,igcm_dust_mass)
1820
1821	pdq(:,:,igcm_dust_number) = pdq(:,:,igcm_dust_number)
1822	& + zdqcloudco2(:,:,igcm_dust_number)
1823
1824	pdq(:,:,igcm_co2) = pdq(:,:,igcm_co2)
1825	& + zdqcloudco2(:,:,igcm_co2)
1826
1827	pdq(:,:,igcm_co2_ice) = pdq(:,:,igcm_co2_ice)
1828	& + zdqcloudco2(:,:,igcm_co2_ice)
1829
1830	pdq(:,:,igcm_ccnco2_mass) = pdq(:,:,igcm_ccnco2_mass)
1831	& + zdqcloudco2(:,:,igcm_ccnco2_mass)
1832
1833	pdq(:,:,igcm_ccnco2_number) = pdq(:,:,igcm_ccnco2_number)
1834	& + zdqcloudco2(:,:,igcm_ccnco2_number)
1835
1836	if (meteo_flux) then
1837	pdq(:,:,igcm_ccnco2_meteor_mass) =
1838	& pdq(:,:,igcm_ccnco2_meteor_mass) +
1839	& zdqcloudco2(:,:,igcm_ccnco2_meteor_mass)
1840
1841	pdq(:,:,igcm_ccnco2_meteor_number) =
1842	& pdq(:,:,igcm_ccnco2_meteor_number)
1843	& + zdqcloudco2(:,:,igcm_ccnco2_meteor_number)
1844	end if
1845	!Update water ice clouds values as well
1846	if (co2useh2o) then
1847	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
1848	& pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
1849	& zdqcloudco2(1:ngrid,1:nlayer,igcm_h2o_ice)
1850	pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
1851	& pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
1852	& zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_mass)
1853	pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
1854	& pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
1855	& zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
1856
1857	pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
1858	& pdq(:,:,igcm_ccnco2_h2o_mass_ice) +
1859	& zdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ice)
1860
1861	pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
1862	& pdq(:,:,igcm_ccnco2_h2o_mass_ccn) +
1863	& zdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ccn)
1864
1865	pdq(:,:,igcm_ccnco2_h2o_number) =
1866	& pdq(:,:,igcm_ccnco2_h2o_number) +
1867	& zdqcloudco2(:,:,igcm_ccnco2_h2o_number)
1868
1869	c Negative values?
1870	where (pq(:,:,igcm_ccn_mass) +
1871	& ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
1872	pdq(:,:,igcm_ccn_mass) =
1873	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1874	pdq(:,:,igcm_ccn_number) =
1875	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1876	end where
1877	c Negative values?
1878	where (pq(:,:,igcm_ccn_number) +
1879	& ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
1880	pdq(:,:,igcm_ccn_mass) =
1881	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1882	pdq(:,:,igcm_ccn_number) =
1883	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1884	end where
1885	where (pq(:,:,igcm_ccnco2_h2o_mass_ice) +
1886	& pq(:,:,igcm_ccnco2_h2o_mass_ccn) +
1887	& (pdq(:,:,igcm_ccnco2_h2o_mass_ice) +
1888	& pdq(:,:,igcm_ccnco2_h2o_mass_ccn)
1889	& )*ptimestep < 0.)
1890	pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
1891	& - pq(:,:,igcm_ccnco2_h2o_mass_ice)
1892	& /ptimestep + 1.e-30
1893	pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
1894	& - pq(:,:,igcm_ccnco2_h2o_mass_ccn)
1895	& /ptimestep + 1.e-30
1896	pdq(:,:,igcm_ccnco2_h2o_number) =
1897	& - pq(:,:,igcm_ccnco2_h2o_number)
1898	& /ptimestep + 1.e-30
1899	end where
1900
1901	where (pq(:,:,igcm_ccnco2_h2o_number) +
1902	& (pdq(:,:,igcm_ccnco2_h2o_number)
1903	& )*ptimestep < 0.)
1904	pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
1905	& - pq(:,:,igcm_ccnco2_h2o_mass_ice)
1906	& /ptimestep + 1.e-30
1907	pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
1908	& - pq(:,:,igcm_ccnco2_h2o_mass_ccn)
1909	& /ptimestep + 1.e-30
1910	pdq(:,:,igcm_ccnco2_h2o_number) =
1911	& - pq(:,:,igcm_ccnco2_h2o_number)
1912	& /ptimestep + 1.e-30
1913	end where
1914	endif ! of if (co2useh2o)
1915	c Negative values?
1916	where (pq(:,:,igcm_ccnco2_mass) +
1917	& ptimestep*pdq(:,:,igcm_ccnco2_mass) < 0.)
1918	pdq(:,:,igcm_ccnco2_mass) =
1919	& - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
1920	pdq(:,:,igcm_ccnco2_number) =
1921	& - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
1922	end where
1923	where (pq(:,:,igcm_ccnco2_number) +
1924	& ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.)
1925	pdq(:,:,igcm_ccnco2_mass) =
1926	& - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
1927	pdq(:,:,igcm_ccnco2_number) =
1928	& - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
1929	end where
1930
1931	c Negative values?
1932	where (pq(:,:,igcm_dust_mass) +
1933	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
1934	pdq(:,:,igcm_dust_mass) =
1935	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1936	pdq(:,:,igcm_dust_number) =
1937	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1938	end where
1939	where (pq(:,:,igcm_dust_number) +
1940	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
1941	pdq(:,:,igcm_dust_mass) =
1942	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1943	pdq(:,:,igcm_dust_number) =
1944	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1945	end where
1946	if (meteo_flux) then
1947	where (pq(:,:,igcm_ccnco2_meteor_mass) +
1948	& ptimestep*pdq(:,:,igcm_ccnco2_meteor_mass) < 0.)
1949	pdq(:,:,igcm_ccnco2_meteor_mass) =
1950	& - pq(:,:,igcm_ccnco2_meteor_mass)/ptimestep + 1.e-30
1951	pdq(:,:,igcm_ccnco2_meteor_number) =
1952	& - pq(:,:,igcm_ccnco2_meteor_number)/ptimestep + 1.e-30
1953	end where
1954	where (pq(:,:,igcm_ccnco2_meteor_number) +
1955	& ptimestep*pdq(:,:,igcm_ccnco2_meteor_number) < 0.)
1956	pdq(:,:,igcm_ccnco2_meteor_mass) =
1957	& - pq(:,:,igcm_ccnco2_meteor_mass)/ptimestep + 1.e-30
1958	pdq(:,:,igcm_ccnco2_meteor_number) =
1959	& - pq(:,:,igcm_ccnco2_meteor_number)/ptimestep + 1.e-30
1960	end where
1961	end if
1962	END IF ! of IF (co2clouds)
1963
1964	c 9b. Aerosol particles
1965	c -------------------
1966	c ----------
1967	c Dust devil :
1968	c ----------
1969	IF(callddevil) then
1970	call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
1971	& zdqdev,zdqsdev)
1972
1973	if (dustbin.ge.1) then
1974	do iq=1,nq
1975	DO l=1,nlayer
1976	DO ig=1,ngrid
1977	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
1978	ENDDO
1979	ENDDO
1980	enddo
1981	do iq=1,nq
1982	DO ig=1,ngrid
1983	DO islope = 1,nslope
1984	dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) +
1985	& zdqsdev(ig,iq)cos(pidef_slope_mean(islope)/180.)
1986	ENDDO
1987	ENDDO
1988	enddo
1989	endif ! of if (dustbin.ge.1)
1990
1991	END IF ! of IF (callddevil)
1992
1993	c -------------
1994	c Sedimentation : acts also on water ice
1995	c -------------
1996	IF (sedimentation) THEN
1997	zdqsed(1:ngrid,1:nlayer,1:nq)=0
1998	zdqssed(1:ngrid,1:nq)=0
1999
2000	c Sedimentation for co2 clouds tracers are inside co2cloud microtimestep
2001	c Zdqssed isn't
2002	call callsedim(ngrid,nlayer,ptimestep,
2003	& zplev,zzlev,zzlay,pt,pdt,
2004	& rdust,rstormdust,rtopdust,
2005	& rice,rsedcloud,rhocloud,
2006	& pq,pdq,zdqsed,zdqssed,nq,
2007	& tau,tauscaling)
2008
2009	c Flux at the surface of co2 ice computed in co2cloud microtimestep
2010	IF (rdstorm) THEN
2011	c Storm dust cannot sediment to the surface
2012	DO ig=1,ngrid
2013	zdqsed(ig,1,igcm_stormdust_mass)=
2014	& zdqsed(ig,1,igcm_stormdust_mass)+
2015	& zdqssed(ig,igcm_stormdust_mass) /
2016	& ((pplev(ig,1)-pplev(ig,2))/g)
2017	zdqsed(ig,1,igcm_stormdust_number)=
2018	& zdqsed(ig,1,igcm_stormdust_number)+
2019	& zdqssed(ig,igcm_stormdust_number) /
2020	& ((pplev(ig,1)-pplev(ig,2))/g)
2021	zdqssed(ig,igcm_stormdust_mass)=0.
2022	zdqssed(ig,igcm_stormdust_number)=0.
2023	ENDDO
2024	ENDIF !rdstorm
2025
2026	DO iq=1, nq
2027	DO l=1,nlayer
2028	DO ig=1,ngrid
2029	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
2030	ENDDO
2031	ENDDO
2032	ENDDO
2033	DO iq=1, nq
2034	DO ig=1,ngrid
2035	DO islope = 1,nslope
2036	dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) +
2037	& zdqssed(ig,iq)cos(pidef_slope_mean(islope)/180.)
2038	ENDDO
2039	ENDDO
2040	ENDDO
2041
2042	END IF ! of IF (sedimentation)
2043
2044	c Add lifted dust to tendancies after sedimentation in the LES (AC)
2045	IF (turb_resolved) THEN
2046	DO iq=1, nq
2047	DO l=1,nlayer
2048	DO ig=1,ngrid
2049	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
2050	ENDDO
2051	ENDDO
2052	ENDDO
2053	DO iq=1, nq
2054	DO ig=1,ngrid
2055	dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:)
2056	ENDDO
2057	ENDDO
2058	ENDIF
2059	c
2060	c 9c. Chemical species
2061	c ------------------
2062
2063	#ifndef MESOSCALE
2064	c --------------
2065	c photochemistry :
2066	c --------------
2067	IF (photochem) then
2068
2069	if (modulo(icount-1,ichemistry).eq.0) then
2070	! compute chemistry every ichemistry physics step
2071
2072	! dust and ice surface area
2073	call surfacearea(ngrid, nlayer, naerkind,
2074	$ ptimestep, zplay, zzlay,
2075	$ pt, pq, pdq, nq,
2076	$ rdust, rice, tau, tauscaling,
2077	$ surfdust, surfice)
2078	! call photochemistry
2079	DO ig = 1,ngrid
2080	qsurf_tmp(ig,:) = qsurf(ig,:,major_slope(ig))
2081	ENDDO
2082	call calchim(ngrid,nlayer,nq,
2083	& ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
2084	$ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
2085	$ zdqcloud,zdqscloud,tau(:,1),
2086	$ qsurf_tmp(:,igcm_co2),
2087	$ pu,pdu,pv,pdv,surfdust,surfice)
2088
2089	endif ! of if (modulo(icount-1,ichemistry).eq.0)
2090
2091	! increment values of tracers:
2092	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
2093	! tracers is zero anyways
2094	DO l=1,nlayer
2095	DO ig=1,ngrid
2096	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
2097	ENDDO
2098	ENDDO
2099	ENDDO ! of DO iq=1,nq
2100
2101	! add condensation tendency for H2O2
2102	if (igcm_h2o2.ne.0) then
2103	DO l=1,nlayer
2104	DO ig=1,ngrid
2105	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
2106	& +zdqcloud(ig,l,igcm_h2o2)
2107	ENDDO
2108	ENDDO
2109	endif
2110
2111	! increment surface values of tracers:
2112	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
2113	! tracers is zero anyways
2114	DO ig=1,ngrid
2115	DO islope = 1,nslope
2116	dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope) +
2117	& zdqschim(ig,iq)cos(pidef_slope_mean(islope)/180.)
2118	ENDDO
2119	ENDDO
2120	ENDDO ! of DO iq=1,nq
2121
2122	! add condensation tendency for H2O2
2123	if (igcm_h2o2.ne.0) then
2124	DO ig=1,ngrid
2125	DO islope = 1,nslope
2126	dqsurf(ig,igcm_h2o2,islope)=dqsurf(ig,igcm_h2o2,islope)+
2127	& zdqscloud(ig,igcm_h2o2)cos(pidef_slope_mean(islope)/180.)
2128	ENDDO
2129	ENDDO
2130	endif
2131
2132	END IF ! of IF (photochem)
2133	#endif
2134
2135
2136	#ifndef MESOSCALE
2137	c-----------------------------------------------------------------------
2138	c 10. THERMOSPHERE CALCULATION
2139	c-----------------------------------------------------------------------
2140
2141	if (callthermos) then
2142	call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol,
2143	$ mu0,ptimestep,ptime,zday,tsurf_meshavg,zzlev,zzlay,
2144	& pt,pq,pu,pv,pdt,pdq,
2145	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff,
2146	$ PhiEscH,PhiEscH2,PhiEscD)
2147
2148	DO l=1,nlayer
2149	DO ig=1,ngrid
2150	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
2151	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)+zdteuv(ig,l)
2152	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
2153	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
2154	DO iq=1, nq
2155	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
2156	ENDDO
2157	ENDDO
2158	ENDDO
2159
2160	endif ! of if (callthermos)
2161	#endif
2162	c-----------------------------------------------------------------------
2163	c 11. Carbon dioxide condensation-sublimation:
2164	c (should be the last atmospherical physical process to be computed)
2165	c -------------------------------------------
2166	IF (tituscap) THEN
2167	!!! get the actual co2 seasonal cap from Titus observations
2168	CALL geticecover(ngrid, 180.*zls/pi,
2169	. 180.longitude/pi, 180.latitude/pi,
2170	. qsurf_tmp(:,igcm_co2) )
2171	qsurf_tmp(:,igcm_co2) = qsurf_tmp(:,igcm_co2) * 10000.
2172	ENDIF
2173
2174
2175	IF (callcond) THEN
2176	zdtc(:,:) = 0.
2177	zdtsurfc(:,:) = 0.
2178	zduc(:,:) = 0.
2179	zdvc(:,:) = 0.
2180	zdqc(:,:,:) = 0.
2181	zdqsc(:,:,:) = 0.
2182	CALL co2condens(ngrid,nlayer,nq,ptimestep,
2183	$ capcal,zplay,zplev,tsurf,pt,
2184	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
2185	$ qsurf(:,igcm_co2,:),albedo,emis,rdust,
2186	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
2187	$ fluxsurf_dn_sw,zls,
2188	$ zdqssed_co2,zcondicea_co2microp,
2189	& zdqsc)
2190
2191	DO iq=1, nq
2192	DO ig=1,ngrid
2193	dqsurf(ig,iq,:)=dqsurf(ig,iq,:)+zdqsc(ig,iq,:)
2194	ENDDO ! (ig)
2195	ENDDO ! (iq)
2196	DO l=1,nlayer
2197	DO ig=1,ngrid
2198	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
2199	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
2200	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
2201	ENDDO
2202	ENDDO
2203	DO ig=1,ngrid
2204	zdtsurf(ig,:) = zdtsurf(ig,:) + zdtsurfc(ig,:)
2205	ENDDO
2206
2207	DO iq=1, nq
2208	DO l=1,nlayer
2209	DO ig=1,ngrid
2210	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
2211	ENDDO
2212	ENDDO
2213	ENDDO
2214
2215	#ifndef MESOSCALE
2216	! update surface pressure
2217	DO ig=1,ngrid
2218	ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
2219	ENDDO
2220	! update pressure levels
2221	DO l=1,nlayer
2222	DO ig=1,ngrid
2223	zplay(ig,l) = aps(l) + bps(l)*ps(ig)
2224	zplev(ig,l) = ap(l) + bp(l)*ps(ig)
2225	ENDDO
2226	ENDDO
2227	zplev(:,nlayer+1) = 0.
2228	! update layers altitude
2229	DO l=2,nlayer
2230	DO ig=1,ngrid
2231	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
2232	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
2233	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
2234	ENDDO
2235	ENDDO
2236	#endif
2237	ENDIF ! of IF (callcond)
2238
2239	c-----------------------------------------------------------------------
2240	c Updating tracer budget on surface
2241	c-----------------------------------------------------------------------
2242	DO iq=1, nq
2243	DO ig=1,ngrid
2244	DO islope = 1,nslope
2245	qsurf(ig,iq,islope)=qsurf(ig,iq,islope)+
2246	& ptimestep*dqsurf(ig,iq,islope)
2247	ENDDO
2248	ENDDO ! (ig)
2249	ENDDO ! (iq)
2250	c-----------------------------------------------------------------------
2251	c 12. Surface and sub-surface soil temperature
2252	c-----------------------------------------------------------------------
2253	c
2254	c
2255	c 12.1 Increment Surface temperature:
2256	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2257
2258	DO ig=1,ngrid
2259	DO islope = 1,nslope
2260	tsurf(ig,islope)=tsurf(ig,islope)+
2261	& ptimestep*zdtsurf(ig,islope)
2262	ENDDO
2263	ENDDO
2264
2265	c Prescribe a cold trap at south pole (except at high obliquity !!)
2266	c Temperature at the surface is set there to be the temperature
2267	c corresponding to equilibrium temperature between phases of CO2
2268
2269
2270	IF (water) THEN
2271	!#ifndef MESOSCALE
2272	! if (caps.and.(obliquit.lt.27.)) then => now done in co2condens
2273	! NB: Updated surface pressure, at grid point 'ngrid', is
2274	! ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
2275	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
2276	! & (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
2277	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095ps(ngrid)))
2278	! endif
2279	!#endif
2280	c -------------------------------------------------------------
2281	c Change of surface albedo in case of ground frost
2282	c everywhere except on the north permanent cap and in regions
2283	c covered by dry ice.
2284	c ALWAYS PLACE these lines after co2condens !!!
2285	c -------------------------------------------------------------
2286	do ig=1,ngrid
2287	do islope = 1,nslope
2288	if ((qsurf(ig,igcm_co2,islope).eq.0).and.
2289	& (qsurf(ig,igcm_h2o_ice,islope)
2290	& .gt.frost_albedo_threshold)) then
2291	if ((watercaptag(ig)).and.(cst_cap_albedo)) then
2292	albedo(ig,1,islope) = albedo_h2o_cap
2293	albedo(ig,2,islope) = albedo_h2o_cap
2294	else
2295	albedo(ig,1,islope) = albedo_h2o_frost
2296	albedo(ig,2,islope) = albedo_h2o_frost
2297	endif !((watercaptag(ig)).and.(cst_cap_albedo)) then
2298	c write(,) "frost thickness", qsurf(ig,igcm_h2o_ice)
2299	c write(,) "physiq.F frost :"
2300	c & ,latitude(ig)180./pi, longitude(ig)180./pi
2301	endif
2302	enddo ! islope
2303	enddo ! of do ig=1,ngrid
2304	ENDIF ! of IF (water)
2305
2306	c
2307	c 12.2 Compute soil temperatures and subsurface heat flux:
2308	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2309	IF (callsoil) THEN
2310	c Thermal inertia feedback
2311	IF (tifeedback) THEN
2312	DO islope = 1,nslope
2313	CALL soil_tifeedback(ngrid,nsoilmx,
2314	s qsurf(:,:,islope),
2315	s inertiesoil_tifeedback(:,:,islope))
2316	ENDDO
2317	CALL soil(ngrid,nsoilmx,.false.,inertiesoil_tifeedback,
2318	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
2319	ELSE
2320	CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
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,inertiesoil,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","",
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","",emis(:,islope))
3141	ENDDO
3142	call write_output("zzlay","Midlayer altitude",
3143	& "m",zzlay(:,:))
3144	call write_output("zzlev","Interlayer altitude",
3145	& "m",zzlev(:,1:nlayer))
3146	call write_output("pphi","Geopotential","m2s-2",
3147	& pphi(:,:))
3148	call write_output("phisfi","Surface geopotential",
3149	& "m2s-2",phisfi(:))
3150	call write_output("tsurf","Surface temperature","K",
3151	& tsurf(:,iflat))
3152	do islope=1,nslope
3153	write(str2(1:2),'(i2.2)') islope
3154	call write_output("tsurf_slope"//str2,
3155	& "Surface temperature","K",
3156	& tsurf(:,islope))
3157	ENDDO
3158	call write_output("ps","surface pressure","Pa",ps(:))
3159	call write_output("co2ice","co2 ice thickness"
3160	& ,"kg.m-2",qsurf(:,igcm_co2,iflat))
3161	do islope=1,nslope
3162	write(str2(1:2),'(i2.2)') islope
3163	call write_output("co2ice_slope"//str2,"co2 ice thickness"
3164	& ,"kg.m-2",qsurf(:,igcm_co2,islope))
3165	ENDDO
3166	call write_output("watercap","Perennial water ice thickness"
3167	& ,"kg.m-2",watercap(:,iflat))
3168	do islope=1,nslope
3169	write(str2(1:2),'(i2.2)') islope
3170	call write_output("watercap_slope"//str2,
3171	& "Perennial water ice thickness"
3172	& ,"kg.m-2",watercap(:,islope))
3173	ENDDO
3174	call write_output("temp_layer1","temperature in layer 1",
3175	& "K",zt(:,1))
3176	call write_output("temp7","temperature in layer 7",
3177	& "K",zt(:,7))
3178	call write_output("fluxsurf_lw","fluxsurf_lw","W.m-2",
3179	& fluxsurf_lw(:,iflat))
3180	do islope=1,nslope
3181	write(str2(1:2),'(i2.2)') islope
3182	call write_output("fluxsurf_lw_slope"//str2,
3183	& "fluxsurf_lw","W.m-2",
3184	& fluxsurf_lw(:,islope))
3185	ENDDO
3186	call write_output("fluxsurf_dn_sw","fluxsurf_dn_sw",
3187	& "W.m-2",fluxsurf_dn_sw_tot(:,iflat))
3188	do islope=1,nslope
3189	write(str2(1:2),'(i2.2)') islope
3190	call write_output("fluxsurf_dn_sw_slope"//str2,
3191	& "fluxsurf_dn_sw",
3192	& "W.m-2",fluxsurf_dn_sw_tot(:,islope))
3193	ENDDO
3194	call write_output("fluxtop_lw","fluxtop_lw","W.m-2",
3195	& fluxtop_lw(:))
3196	call write_output("fluxtop_up_sw","fluxtop_up_sw","W.m-2",
3197	& fluxtop_up_sw_tot(:))
3198	call write_output("temp","temperature","K",zt(:,:))
3199	call write_output("Sols","Time","sols",zday)
3200	call write_output("Ls","Solar longitude","deg",
3201	& zls*180./pi)
3202	call write_output("u","Zonal wind","m.s-1",zu(:,:))
3203	call write_output("v","Meridional wind","m.s-1",zv(:,:))
3204	call write_output("w","Vertical wind","m.s-1",pw(:,:))
3205	call write_output("rho","density","kg.m-3",rho(:,:))
3206	call write_output("pressure","Pressure","Pa",zplay(:,:))
3207	call write_output("zplev","Interlayer pressure","Pa",
3208	& zplev(:,1:nlayer))
3209	call write_output('sw_htrt','sw heat. rate',
3210	& 'K/s',zdtsw(:,:))
3211	call write_output('lw_htrt','lw heat. rate',
3212	& 'K/s',zdtlw(:,:))
3213	call write_output("local_time","Local time",
3214	& 'sol',local_time(:))
3215	if (.not.activice) then
3216	CALL write_output('tauTESap',
3217	& 'tau abs 825 cm-1',
3218	& '',tauTES(:))
3219	else
3220	CALL write_output('tauTES',
3221	& 'tau abs 825 cm-1',
3222	& '',taucloudtes(:))
3223	endif
3224	#else
3225	!!! this is to ensure correct initialisation of mesoscale model
3226	call write_output("tsurf","Surface temperature","K",
3227	& tsurf(:,iflat))
3228	call write_output("ps","surface pressure","Pa",ps(:))
3229	call write_output("co2ice","co2 ice thickness","kg.m-2",
3230	& qsurf(:,igcm_co2,iflat))
3231	call write_output("temp","temperature","K",zt(:,:))
3232	call write_output("u","Zonal wind","m.s-1",zu(:,:))
3233	call write_output("v","Meridional wind","m.s-1",zv(:,:))
3234	call write_output("emis","Surface emissivity","",
3235	& emis(:,iflat))
3236	call write_output("tsoil","Soil temperature",
3237	& "K",tsoil(:,:,iflat))
3238	call write_output("inertiedat","Soil inertia",
3239	& "K",inertiedat(:,:))
3240	#endif
3241
3242	c ----------------------------------------------------------
3243	c Outputs of the CO2 cycle
3244	c ----------------------------------------------------------
3245
3246	if (igcm_co2.ne.0) then
3247	call write_output("co2","co2 mass mixing ratio",
3248	& "kg.kg-1",zq(:,:,igcm_co2))
3249
3250	if (co2clouds) then
3251	call write_output('ccnqco2','CCNco2 mmr',
3252	& 'kg.kg-1',zq(:,:,igcm_ccnco2_mass))
3253
3254	call write_output('ccnNco2','CCNco2 number',
3255	& 'part.kg-1',zq(:,:,igcm_ccnco2_number))
3256
3257	call write_output('co2_ice','co2_ice mmr in atm',
3258	& 'kg.kg-1',zq(:,:,igcm_co2_ice))
3259
3260	call write_output("mtotco2","total mass atm of co2",
3261	& "kg.m-2",mtotco2(:))
3262	call write_output("icetotco2","total mass atm of co2 ice",
3263	& "kg.m-2", icetotco2(:))
3264	call write_output("vaptotco2","total mass atm of co2 "//
3265	& "vapor","kg.m-2", vaptotco2(:))
3266	if (co2useh2o) then
3267	call write_output('ccnqco2_h2o_m_ice',
3268	& 'CCNco2_h2o_mass_ice mmr',
3269	& 'kg.kg-1',zq(:,:,igcm_ccnco2_h2o_mass_ice))
3270
3271	call write_output('ccnqco2_h2o_m_ccn',
3272	& 'CCNco2_h2o_mass_ccn mmr',
3273	& 'kg.kg-1',zq(:,:,igcm_ccnco2_h2o_mass_ccn))
3274
3275	call write_output('ccnNco2_h2o','CCNco2_h2o number',
3276	& 'part.kg-1',zq(:,:,igcm_ccnco2_h2o_number))
3277	end if
3278
3279	if (meteo_flux) then
3280	call write_output('ccnqco2_meteor','CCNco2_meteor mmr',
3281	& 'kg.kg-1',zq(:,:,igcm_ccnco2_meteor_mass))
3282
3283	call write_output('ccnNco2_meteor','CCNco2_meteor number',
3284	& 'part.kg-1',zq(:,:,igcm_ccnco2_meteor_number))
3285	end if
3286
3287	end if ! of if (co2clouds)
3288	end if ! of if (igcm_co2.ne.0)
3289
3290	! Output He tracer, if there is one
3291	if (igcm_he.ne.0) then
3292	call write_output("he","helium mass mixing ratio",
3293	& "kg/kg",zq(:,:,igcm_he))
3294	vmr = zq(1:ngrid,1:nlayer,igcm_he)
3295	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_he)
3296	call write_output('vmr_he','helium vol. mixing ratio',
3297	& 'mol/mol',vmr(:,:))
3298	end if
3299
3300	c ----------------------------------------------------------
3301	c Outputs of the water cycle
3302	c ----------------------------------------------------------
3303	if (water) then
3304	#ifdef MESOINI
3305	!!!! waterice = q01, voir readmeteo.F90
3306	call write_output('q01',noms(igcm_h2o_ice),
3307	& 'kg/kg',
3308	& zq(:,:,igcm_h2o_ice))
3309	!!!! watervapor = q02, voir readmeteo.F90
3310	call write_output('q02',noms(igcm_h2o_vap),
3311	& 'kg/kg',
3312	& zq(:,:,igcm_h2o_vap))
3313	!!!! surface waterice qsurf02 (voir readmeteo)
3314	call write_output('qsurf02','surface tracer',
3315	& 'kg.m-2',
3316	& qsurf(:,igcm_h2o_ice,iflat))
3317	do islope=1,nslope
3318	write(str2(1:2),'(i2.2)') islope
3319	call write_output('qsurf02_slope'//str2,
3320	& 'surface tracer','kg.m-2',
3321	& qsurf(:,igcm_h2o_ice,islope))
3322	ENDDO
3323	#endif
3324	call write_output('mtot',
3325	& 'total mass of water vapor',
3326	& 'kg/m2',mtot(:))
3327	call write_output('icetot',
3328	& 'total mass of water ice',
3329	& 'kg/m2',icetot(:))
3330	vmr = zq(1:ngrid,1:nlayer,igcm_h2o_ice)
3331	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
3332	call write_output('vmr_h2oice','h2o ice vmr',
3333	& 'mol/mol',vmr(:,:))
3334	vmr = zq(1:ngrid,1:nlayer,igcm_h2o_vap)
3335	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
3336	call write_output('vmr_h2ovap','h2o vap vmr',
3337	& 'mol/mol',vmr(:,:))
3338	call write_output('reffice',
3339	& 'Mean reff',
3340	& 'm',rave(:))
3341	call write_output('h2o_ice','h2o_ice','kg/kg',
3342	& zq(:,:,igcm_h2o_ice))
3343	call write_output('h2o_vap','h2o_vap','kg/kg',
3344	& zq(:,:,igcm_h2o_vap))
3345
3346	if (hdo) then
3347	vmr=zq(1:ngrid,1:nlayer,igcm_hdo_ice)
3348	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_ice)
3349	CALL write_output('vmr_hdoice','hdo ice vmr',
3350	& 'mol/mol',vmr(:,:))
3351	vmr=zq(1:ngrid,1:nlayer,igcm_hdo_vap)
3352	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_vap)
3353	CALL write_output('vmr_hdovap','hdo vap vmr',
3354	& 'mol/mol',vmr(:,:))
3355	call write_output('hdo_ice','hdo_ice','kg/kg',
3356	& zq(:,:,igcm_hdo_ice))
3357	call write_output('hdo_vap','hdo_vap','kg/kg',
3358	& zq(:,:,igcm_hdo_vap))
3359
3360	CALL write_output('mtotD',
3361	& 'total mass of HDO vapor',
3362	& 'kg/m2',mtotD(:))
3363	CALL write_output('icetotD',
3364	& 'total mass of HDO ice',
3365	& 'kg/m2',icetotD(:))
3366
3367	C Calculation of the D/H ratio
3368	do l=1,nlayer
3369	do ig=1,ngrid
3370	if (zq(ig,l,igcm_h2o_vap).gt.qperemin) then
3371	DoH_vap(ig,l) = ( zq(ig,l,igcm_hdo_vap)/
3372	& zq(ig,l,igcm_h2o_vap) )1./(2.155.76e-6)
3373	else
3374	DoH_vap(ig,l) = 0.
3375	endif
3376	enddo
3377	enddo
3378
3379	do l=1,nlayer
3380	do ig=1,ngrid
3381	if (zq(ig,l,igcm_h2o_ice).gt.qperemin) then
3382	DoH_ice(ig,l) = ( zq(ig,l,igcm_hdo_ice)/
3383	& zq(ig,l,igcm_h2o_ice) )/(2.*155.76e-6)
3384	else
3385	DoH_ice(ig,l) = 0.
3386	endif
3387	enddo
3388	enddo
3389
3390	CALL write_output('DoH_vap',
3391	& 'D/H ratio in vapor',
3392	& ' ',DoH_vap(:,:))
3393	CALL write_output('DoH_ice',
3394	& 'D/H ratio in ice',
3395	& '',DoH_ice(:,:))
3396
3397	endif !hdo
3398
3399	!A. Pottier
3400	! CALL write_output('rmoym',
3401	! & 'alternative reffice',
3402	! & 'm',rave2(:))
3403	call write_output('h2o_saturation',
3404	& 'h2o vap saturation ratio','',satu(:,:))
3405	if (scavenging) then
3406	CALL write_output("Nccntot",
3407	& "condensation nuclei","Nbr/m2",
3408	& Nccntot(:))
3409	CALL write_output("Mccntot",
3410	& "mass condensation nuclei","kg/m2",
3411	& Mccntot(:))
3412	endif
3413	call write_output('rice','Water ice particle size',
3414	& 'm',rice(:,:))
3415	call write_output('h2o_ice_s',
3416	& 'surface h2o_ice',
3417	& 'kg.m-2',qsurf(:,igcm_h2o_ice,iflat))
3418	do islope=1,nslope
3419	write(str2(1:2),'(i2.2)') islope
3420	call write_output('h2o_ice_s_slope'//str2,
3421	& 'surface h2o_ice',
3422	& 'kg.m-2',qsurf(:,igcm_h2o_ice,islope))
3423	ENDDO
3424	if (hdo) then
3425	call write_output('hdo_ice_s',
3426	& 'surface hdo_ice',
3427	& 'kg.m-2',qsurf(:,igcm_hdo_ice,iflat))
3428	do islope=1,nslope
3429	write(str2(1:2),'(i2.2)') islope
3430	call write_output('hdo_ice_s_slope'//str2,
3431	& 'surface hdo_ice',
3432	& 'kg.m-2',qsurf(:,igcm_hdo_ice,islope))
3433	ENDDO
3434
3435	do ig=1,ngrid
3436	if (qsurf_meshavg(ig,igcm_h2o_ice).gt.qperemin) then
3437	DoH_surf(ig) = 0.5*( qsurf_meshavg(ig,igcm_hdo_ice)/
3438	& qsurf_meshavg(ig,igcm_h2o_ice) )/155.76e-6
3439	else
3440	DoH_surf(ig) = 0.
3441	endif
3442	enddo
3443
3444	call write_output('DoH_surf',
3445	& 'surface D/H',
3446	& '',DoH_surf(:))
3447	endif ! hdo
3448
3449	CALL write_output('albedo',
3450	& 'albedo',
3451	& '',albedo(:,1,iflat))
3452	do islope=1,nslope
3453	write(str2(1:2),'(i2.2)') islope
3454	CALL write_output('albedo_slope'//str2,
3455	& 'albedo',
3456	& '',albedo(:,1,islope))
3457	ENDDO
3458	if (tifeedback) then
3459	call write_output("soiltemp",
3460	& "Soil temperature","K",
3461	& tsoil(:,:,iflat))
3462	call write_output('soilti',
3463	& 'Soil Thermal Inertia',
3464	& 'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,iflat))
3465	do islope=1,nslope
3466	write(str2(1:2),'(i2.2)') islope
3467	call write_output('soilti_slope'//str2,
3468	& 'Soil Thermal Inertia',
3469	& 'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,islope))
3470	ENDDO
3471	endif
3472	!A. Pottier
3473	if (CLFvarying) then !AP14 nebulosity
3474	call write_output('totcloudfrac',
3475	& 'Total cloud fraction',
3476	& ' ',totcloudfrac(:))
3477	end if !clf varying
3478	end if !(water)
3479
3480	c ----------------------------------------------------------
3481	c Outputs of the dust cycle
3482	c ----------------------------------------------------------
3483
3484	call write_output('tau_pref_scenario',
3485	& 'Prescribed visible dust optical depth at 610Pa',
3486	& 'NU',tau_pref_scenario(:))
3487
3488	call write_output('tau_pref_gcm',
3489	& 'Visible dust optical depth at 610Pa in the GCM',
3490	& 'NU',tau_pref_gcm(:))
3491
3492	if (reff_driven_IRtoVIS_scenario) then
3493	call write_output('IRtoVIScoef',
3494	& 'Conversion coeff for dust tau from abs9.3um to ext0.67um',
3495	& '/',IRtoVIScoef(:))
3496	endif
3497
3498	if (dustbin.ne.0) then
3499
3500	#ifndef MESOINI
3501	if (doubleq) then
3502	call write_output('dqsdust',
3503	& 'deposited surface dust mass',
3504	& 'kg.m-2.s-1',dqdustsurf(:))
3505	call write_output('dqndust',
3506	& 'deposited surface dust number',
3507	& 'number.m-2.s-1',dndustsurf(:))
3508	call write_output('reffdust','reffdust',
3509	& 'm',rdust(:,:)*ref_r0)
3510	call write_output('dustq','Dust mass mr',
3511	& 'kg/kg',qdust(:,:))
3512	call write_output('dustN','Dust number',
3513	& 'part/kg',ndust(:,:))
3514
3515	select case (trim(dustiropacity))
3516	case ("tes")
3517	call write_output('dsodust_TES',
3518	& 'density scaled extinction opacity of std dust at 9.3um(TES)',
3519	& 'm2.kg-1',dsodust(:,:))
3520	call write_output('dso_TES',
3521	& 'density scaled extinction opacity of all dust at 9.3um(TES)',
3522	& 'm2.kg-1',dsodust(:,:)+dsords(:,:)+dsotop(:,:))
3523	case ("mcs")
3524	call write_output('dsodust',
3525	& 'density scaled extinction opacity of std dust at 21.6um(MCS)',
3526	& 'm2.kg-1',dsodust(:,:))
3527	call write_output('dso',
3528	& 'density scaled extinction opacity of all dust at 21.6um(MCS)',
3529	& 'm2.kg-1',dsodust(:,:)+dsords(:,:)+dsotop(:,:))
3530	end select
3531	else ! (doubleq=.false.)
3532	do iq=1,dustbin
3533	write(str2(1:2),'(i2.2)') iq
3534	call write_output('q'//str2,'mix. ratio',
3535	& 'kg/kg',zq(:,:,iq))
3536	call write_output('qsurf'//str2,'qsurf',
3537	& 'kg.m-2',qsurf(:,iq,iflat))
3538	do islope=1,nslope
3539	write(str2(1:2),'(i2.2)') islope
3540	call write_output('qsurf_slope'//str2,'qsurf',
3541	& 'kg.m-2',qsurf(:,iq,islope))
3542	ENDDO
3543	end do
3544	endif ! (doubleq)
3545
3546	if (rdstorm) then ! writediagfi tendencies stormdust tracers
3547	call write_output('reffstormdust','reffstormdust',
3548	& 'm',rstormdust(:,:)*ref_r0)
3549	call write_output('mstormdtot',
3550	& 'total mass of stormdust only',
3551	& 'kg.m-2',mstormdtot(:))
3552	call write_output('mdusttot',
3553	& 'total mass of dust only',
3554	& 'kg.m-2',mdusttot(:))
3555	call write_output('rdsdqsdust',
3556	& 'deposited surface stormdust mass',
3557	& 'kg.m-2.s-1',rdsdqdustsurf(:))
3558	call write_output('rdsdustq','storm Dust mass mr',
3559	& 'kg/kg',rdsqdust(:,:))
3560	call write_output('rdsdustqmodel','storm Dust massmr',
3561	& 'kg/kg',pq(:,:,igcm_stormdust_mass))
3562	call write_output('rdsdustN','storm Dust number',
3563	& 'part/kg',rdsndust(:,:))
3564	call write_output("stormfract",
3565	& "fraction of the mesh, with stormdust","none",
3566	& totstormfract(:))
3567	! call write_output('qsurf',
3568	! & 'stormdust injection',
3569	! & 'kg.m-2',qsurf(:,igcm_stormdust_mass,iflat))
3570	! do islope=1,nslope
3571	! write(str2(1:2),'(i2.2)') islope
3572	! call write_output('qsurf_slope'//str2,
3573	! & 'stormdust injection',
3574	! & 'kg.m-2',qsurf(:,igcm_stormdust_mass,islope))
3575	! ENDDO
3576	! call write_output('pdqsurf',
3577	! & 'tendancy stormdust mass at surface',
3578	! & 'kg.m-2',dqsurf(:,igcm_stormdust_mass,iflat))
3579	! do islope=1,nslope
3580	! write(str2(1:2),'(i2.2)') islope
3581	! call write_output('pdqsurf_slope'//str2,
3582	! & 'tendancy stormdust mass at surface',
3583	! & 'kg.m-2',dqsurf(:,igcm_stormdust_mass,islope))
3584	! ENDDO
3585	call write_output('wspeed_stormdust',
3586	& 'vertical velocity of stormdust',
3587	& 'm/s',wspeed(:,:))
3588	call write_output('zdqsed_dust_mass'
3589	& ,'sedimentation tendency of background dust mmr'
3590	& ,'kg/kg.s-1',
3591	& zdqsed(:,:,igcm_dust_mass))
3592	call write_output('zdqssed_dust_mass'
3593	& ,'sedimentation tendency of background dust on surface'
3594	& ,'kg.m-2.s-1',
3595	& zdqssed(:,igcm_dust_mass))
3596	call write_output('zdqsed_stormdust_mass'
3597	& ,'sedimentation tendency of stormdust dust mmr'
3598	& ,'kg/kg.s-1',
3599	& zdqsed(:,:,igcm_stormdust_mass))
3600	call write_output('zdqsed_dust_number'
3601	& ,'sedimentation tendency of background dust number'
3602	& ,'nbr/kg.s-1',
3603	& zdqsed(:,:,igcm_dust_number))
3604	call write_output('rdust','rdust',
3605	& 'm',rdust(:,:))
3606	call write_output('rstormdust','rstormdust',
3607	& 'm',rstormdust(:,:))
3608
3609	select case (trim(dustiropacity))
3610	case ("tes")
3611	call write_output('dsords_TES',
3612	& 'density scaled extinction opacity of stormdust at 9.3um(TES)',
3613	& 'm2.kg-1',dsords(:,:))
3614	case ("mcs")
3615	call write_output('dsords',
3616	& 'density scaled extinction opacity of stormdust at 21.6um(MCS)',
3617	& 'm2.kg-1',dsords(:,:))
3618	end select
3619	endif ! (rdstorm)
3620
3621	if (topflows) then
3622	call write_output('refftopdust',
3623	& 'Topdust dust effective radius',
3624	& 'm',rtopdust(:,:)*ref_r0)
3625	call write_output('topdustq','top Dust mass mr',
3626	& 'kg/kg',pq(:,:,igcm_topdust_mass))
3627	call write_output('topdustN','top Dust number',
3628	& 'part/kg',pq(:,:,igcm_topdust_number))
3629	select case (trim(dustiropacity))
3630	case ("tes")
3631	call write_output('dsotop_TES',
3632	& 'density scaled extinction opacity of topdust at 9.3um(TES)',
3633	& 'm2.kg-1',dsotop(:,:))
3634	case ("mcs")
3635	call write_output('dsotop',
3636	& 'density scaled extinction opacity of topdust at 21.6um(MCS)',
3637	& 'm2.kg-1',dsotop(:,:))
3638	end select
3639	endif ! (topflows)
3640
3641	if (dustscaling_mode==2) then
3642	call write_output("dust_rad_adjust",
3643	& "radiative adjustment coefficient for dust",
3644	& "",dust_rad_adjust(:))
3645	endif
3646
3647	! if (scavenging) then ! these outputs should be in the scavenging routine
3648	! call write_output('ccnq','CCN mass mr',
3649	! & 'kg/kg',qccn(:,:))
3650	! call write_output('ccnN','CCN number',
3651	! & 'part/kg',nccn(:,:))
3652	! call write_output('surfccnq','Surf nuclei mass mr',
3653	! & 'kg.m-2',qsurf(:,igcm_ccn_mass,iflat))
3654	! do islope=1,nslope
3655	! write(str2(1:2),'(i2.2)') islope
3656	! call write_output('surfccnq_slope'//str2,
3657	! & 'Surf nuclei mass mr',
3658	! & 'kg.m-2',qsurf(:,igcm_ccn_mass,islope))
3659	! ENDDO
3660	! call write_output('surfccnN','Surf nuclei number',
3661	! & 'kg.m-2',qsurf(:,igcm_ccn_number,iflat))
3662	! do islope=1,nslope
3663	! write(str2(1:2),'(i2.2)') islope
3664	! call write_output('surfccnN_slope'//str2,
3665	! & 'Surf nuclei number',
3666	! & 'kg.m-2',qsurf(:,igcm_ccn_number,islope))
3667	! ENDDO
3668	! endif ! (scavenging)
3669
3670	#else
3671	! !!! to initialize mesoscale we need scaled variables
3672	! !!! because this must correspond to starting point for tracers
3673	! call write_output('dustq','Dust mass mr',
3674	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass))
3675	! call write_output('dustN','Dust number',
3676	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number))
3677	! call write_output('ccn','Nuclei mass mr',
3678	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass))
3679	! call write_output('ccnN','Nuclei number',
3680	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number))
3681	if (freedust) then
3682	call write_output('dustq','Dust mass mr',
3683	& 'kg/kg',qdust)
3684	call write_output('dustN','Dust number',
3685	& 'part/kg',ndust)
3686	call write_output('ccn','CCN mass mr',
3687	& 'kg/kg',qccn)
3688	call write_output('ccnN','CCN number',
3689	& 'part/kg',nccn)
3690	else
3691	call write_output('dustq','Dust mass mr',
3692	& 'kg/kg',pq(:,:,igcm_dust_mass))
3693	call write_output('dustN','Dust number',
3694	& 'part/kg',pq(:,:,igcm_dust_number))
3695	call write_output('ccn','Nuclei mass mr',
3696	& 'kg/kg',pq(:,:,igcm_ccn_mass))
3697	call write_output('ccnN','Nuclei number',
3698	& 'part/kg',pq(:,:,igcm_ccn_number))
3699	endif
3700	#endif
3701
3702	end if ! (dustbin.ne.0)
3703
3704	c ----------------------------------------------------------
3705	c Thermospheric outputs
3706	c ----------------------------------------------------------
3707
3708	if(callthermos) then
3709
3710	call write_output("quv","UV heating","K/s",
3711	$ zdteuv(:,:))
3712	call write_output("cond","Thermal conduction","K/s",
3713	$ zdtconduc(:,:))
3714
3715	!H, H2 and D escape fluxes
3716
3717	call write_output("PhiH","H escape flux","s-1",
3718	$ PhiEscH)
3719	call write_output("PhiH2","H2 escape flux","s-1",
3720	$ PhiEscH2)
3721	call write_output("PhiD","D escape flux","s-1",
3722	$ PhiEscD)
3723
3724	endif !(callthermos)
3725
3726	call write_output("q15um","15 um cooling","K/s",
3727	$ zdtnlte(:,:))
3728	call write_output("qnir","NIR heating","K/s",
3729	$ zdtnirco2(:,:))
3730
3731	c ----------------------------------------------------------
3732	c ----------------------------------------------------------
3733	c PBL OUTPUS
3734	c ----------------------------------------------------------
3735	c ----------------------------------------------------------
3736
3737	c ----------------------------------------------------------
3738	c Outputs of thermals
3739	c ----------------------------------------------------------
3740	if (calltherm) then
3741	call write_output('lmax_th',
3742	& 'index of vertical extension of thermals',
3743	& 'grid level',lmax_th_out(:))
3744	call write_output('zmax_th',
3745	& 'vertical extension of thermals','m',
3746	& zmax_th(:))
3747	call write_output('hfmax_th',
3748	& 'maximum heat flux in thermals','K.m/s',
3749	& hfmax_th(:))
3750	call write_output('wstar',
3751	& 'maximum thermals vertical velocity','m/s',
3752	& wstar(:))
3753	end if
3754
3755	c ----------------------------------------------------------
3756	c ----------------------------------------------------------
3757	c END OF PBL OUTPUS
3758	c ----------------------------------------------------------
3759	c ----------------------------------------------------------
3760
3761
3762	c ----------------------------------------------------------
3763	c Output in netcdf file "diagsoil.nc" for subterranean
3764	c variables (output every "ecritphy", as for writediagfi)
3765	c ----------------------------------------------------------
3766
3767	! Write soil temperature
3768	call write_output("soiltemp","Soil temperature","K",
3769	& tsoil(:,:,iflat))
3770	do islope=1,nslope
3771	write(str2(1:2),'(i2.2)') islope
3772	call write_output("soiltemp_slope"//str2,
3773	& "Soil temperature","K",
3774	& tsoil(:,:,islope))
3775	ENDDO
3776
3777	!PREVIOUSLY IN 1D ONLY
3778	call write_output("dtrad","rad. heat. rate",
3779	& "K.s-1",dtrad(:,:))
3780
3781	if (rdstorm) then
3782	call write_output('aerosol_dust','opacity of env. dust',''
3783	& ,aerosol(:,:,iaer_dust_doubleq))
3784	call write_output('aerosol_stormdust',
3785	& 'opacity of storm dust',''
3786	& ,aerosol(:,:,iaer_stormdust_doubleq))
3787	call write_output('dqsdifdustq',
3788	&'tendency due to vertical diffusion of background dust on surface'
3789	& ,'kg.m-2.s-1',zdqsdif(:,igcm_dust_mass,iflat))
3790	do islope=1,nslope
3791	write(str2(1:2),'(i2.2)') islope
3792	call write_output('dqsdifdustq_slope'//str2,
3793	&'tendency due to vertical diffusion of background dust on surface'
3794	& ,'kg.m-2.s-1',zdqsdif(:,igcm_dust_mass,islope))
3795	ENDDO
3796	call write_output('dqsdifrdsq',
3797	& 'tendency due to vertical diffusion of stormdust on surface',
3798	& 'kg.m-2.s-1',zdqsdif(:,igcm_stormdust_mass,iflat))
3799	do islope=1,nslope
3800	write(str2(1:2),'(i2.2)') islope
3801	call write_output('dqsdifrdsq_slope'//str2,
3802	& 'tendency due to vertical diffusion of stormdust on surface',
3803	& 'kg.m-2.s-1',zdqsdif(:,igcm_stormdust_mass,islope))
3804	ENDDO
3805	endif !(rdstorm)
3806
3807	if(water) then
3808	if (.not.scavenging) then
3809	call write_output('zdqcloud_ice','cloud ice',
3810	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_ice))
3811	call write_output('zdqcloud_vap','cloud vap',
3812	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_vap))
3813	call write_output('zdqcloud','cloud ice',
3814	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_ice)
3815	& +zdqcloud(:,:,igcm_h2o_vap))
3816	IF (hdo) THEN
3817	call write_output('zdqcloud_iceD','cloud ice hdo',
3818	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_hdo_ice))
3819	call write_output('zdqcloud_vapD','cloud vap hdo',
3820	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_hdo_vap))
3821	ENDIF ! hdo
3822	endif !not.scavenging
3823	ENDIF ! of IF (water)
3824	!PREVIOUSLY IN 1D ONLY
3825
3826	c ==========================================================
3827	c END OF WRITEDIAGFI
3828	c ==========================================================
3829	#endif
3830	! of ifdef MESOSCALE
3831
3832	c ELSE ! if(ngrid.eq.1)
3833
3834	c#ifndef MESOSCALE
3835	c write(*,
3836	c & '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)')
3837	c & zls*180./pi,odpref,tau_pref_scenario
3838	c#endif
3839
3840	c END IF ! if(ngrid.ne.1)
3841
3842
3843	! test for co2 conservation with co2 microphysics
3844	if (igcm_co2_ice.ne.0) then
3845	co2totB = 0. ! added by C.M.
3846	do ig=1,ngrid
3847	do l=1,nlayer
3848	co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
3849	& (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice)
3850	& +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
3851	enddo
3852	co2totB = co2totB + qsurf(ig,igcm_co2,iflat)
3853	enddo
3854	else
3855	co2totB = 0. ! added by C.M.
3856	do ig=1,ngrid
3857	do l=1,nlayer
3858	co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
3859	& (pq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2)*ptimestep)
3860	enddo
3861	co2totB = co2totB + qsurf(ig,igcm_co2,iflat)
3862	enddo
3863	endif ! of if (igcm_co2_ice.ne.0)
3864	co2conservation = (co2totA-co2totB)/co2totA
3865	call write_output( 'co2conservation',
3866	& 'Total CO2 mass conservation in physic',
3867	& 'kg', co2conservation)
3868	! XIOS outputs
3869	#ifdef CPP_XIOS
3870	! Send fields to XIOS: (NB these fields must also be defined as
3871	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
3872
3873	CALL send_xios_field("controle",tab_cntrl_mod,1)
3874
3875	CALL send_xios_field("ap",ap,1)
3876	CALL send_xios_field("bp",bp,1)
3877	CALL send_xios_field("aps",aps,1)
3878	CALL send_xios_field("bps",bps,1)
3879
3880	if (lastcall.and.is_omp_master) then
3881	write(,) "physiq lastcall: call xios_context_finalize"
3882	call xios_context_finalize
3883	endif
3884	#endif
3885
3886	if (check_physics_outputs) then
3887	! Check the validity of updated fields at the end of the physics step
3888	call check_physics_fields("end of physiq:",zt,zu,zv,zplev,zq)
3889	endif
3890
3891	icount=icount+1
3892
3893	END SUBROUTINE physiq
3894
3895	END MODULE physiq_mod

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

Download in other formats:

Original Format