Context Navigation

physiq_mod.F @ 2827

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

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

Download in other formats:

Original Format