Context Navigation

physiq_mod.F @ 2800

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

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

Download in other formats:

Original Format