Context Navigation

physiq_mod.F @ 2814

Last change on this file since 2814 was 2808, checked in by emillour, 3 years ago
Mars GCM: Bug fix in jthermcalc_e107 (wrong indexes used in some computations involving NO2 and H2). AM
File size: 163.4 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,"pphi","Geopotential","m2s-2",3,
3002	& pphi)
3003	call writediagfi(ngrid,"phisfi","Surface geopotential",
3004	& "m2s-2",2,phisfi)
3005	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
3006	& tsurf)
3007	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
3008	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
3009	& ,"kg.m-2",2,co2ice)
3010	call WRITEDIAGFI(ngrid,"watercap","Water ice thickness"
3011	& ,"kg.m-2",2,watercap)
3012
3013	call WRITEDIAGFI(ngrid,"temp_layer1","temperature in layer 1",
3014	& "K",2,zt(1,1))
3015	call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
3016	& "K",2,zt(1,7))
3017	call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
3018	& fluxsurf_lw)
3019	call WRITEDIAGFI(ngrid,"fluxsurf_dn_sw","fluxsurf_dn_sw",
3020	& "W.m-2",2,fluxsurf_dn_sw_tot)
3021	call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
3022	& fluxtop_lw)
3023	call WRITEDIAGFI(ngrid,"fluxtop_up_sw","fluxtop_up_sw","W.m-2",
3024	& 2,fluxtop_up_sw_tot)
3025	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
3026	call WRITEDIAGFI(ngrid,"Sols","Time","sols",0,[zday])
3027	call WRITEDIAGFI(ngrid,"Ls","Solar longitude","deg",
3028	& 0,[zls*180./pi])
3029	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
3030	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
3031	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
3032	call WRITEDIAGFI(ngrid,"rho","density","kg.m-3",3,rho)
3033	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
3034	c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
3035	call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,zplay)
3036	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
3037	c & zstress)
3038	call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
3039	& 'w.m-2',3,zdtsw)
3040	call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
3041	& 'w.m-2',3,zdtlw)
3042	call writediagfi(ngrid,"local_time","Local time",
3043	& 'sol',2,local_time)
3044	if (.not.activice) then
3045	CALL WRITEDIAGFI(ngrid,'tauTESap',
3046	& 'tau abs 825 cm-1',
3047	& '',2,tauTES)
3048	else
3049	CALL WRITEDIAGFI(ngrid,'tauTES',
3050	& 'tau abs 825 cm-1',
3051	& '',2,taucloudtes)
3052	endif
3053	#else
3054	!!! this is to ensure correct initialisation of mesoscale model
3055	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
3056	& tsurf)
3057	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
3058	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
3059	& co2ice)
3060	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
3061	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
3062	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
3063	call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
3064	& emis)
3065	call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
3066	& "K",3,tsoil)
3067	call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
3068	& "K",3,inertiedat)
3069	#endif
3070
3071	c ----------------------------------------------------------
3072	c Outputs of the CO2 cycle
3073	c ----------------------------------------------------------
3074
3075	if (tracer.and.(igcm_co2.ne.0)) then
3076	call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
3077	& "kg.kg-1",3,zq(:,:,igcm_co2))
3078
3079	if (co2clouds) then
3080	call WRITEDIAGFI(ngrid,'ccnqco2','CCNco2 mmr',
3081	& 'kg.kg-1',3,zq(:,:,igcm_ccnco2_mass))
3082
3083	call WRITEDIAGFI(ngrid,'ccnNco2','CCNco2 number',
3084	& 'part.kg-1',3,zq(:,:,igcm_ccnco2_number))
3085
3086	call WRITEDIAGFI(ngrid,'co2_ice','co2_ice mmr in atm',
3087	& 'kg.kg-1', 3, zq(:,:,igcm_co2_ice))
3088
3089	call WRITEDIAGFI(ngrid,"mtotco2","total mass atm of co2",
3090	& "kg.m-2",2, mtotco2)
3091	call WRITEDIAGFI(ngrid,"icetotco2","total mass atm of co2 ice",
3092	& "kg.m-2", 2, icetotco2)
3093	call WRITEDIAGFI(ngrid,"vaptotco2","total mass atm of co2
3094	& vapor","kg.m-2", 2, vaptotco2)
3095	call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
3096	& emis)
3097	if (co2useh2o) then
3098	call WRITEDIAGFI(ngrid,'ccnqco2_h2o_m_ice',
3099	& 'CCNco2_h2o_mass_ice mmr',
3100	& 'kg.kg-1',3,zq(:,:,igcm_ccnco2_h2o_mass_ice))
3101
3102	call WRITEDIAGFI(ngrid,'ccnqco2_h2o_m_ccn',
3103	& 'CCNco2_h2o_mass_ccn mmr',
3104	& 'kg.kg-1',3,zq(:,:,igcm_ccnco2_h2o_mass_ccn))
3105
3106	call WRITEDIAGFI(ngrid,'ccnNco2_h2o','CCNco2_h2o number',
3107	& 'part.kg-1',3,zq(:,:,igcm_ccnco2_h2o_number))
3108	end if
3109
3110	if (meteo_flux) then
3111	call WRITEDIAGFI(ngrid,'ccnqco2_meteor','CCNco2_meteor mmr',
3112	& 'kg.kg-1',3,zq(:,:,igcm_ccnco2_meteor_mass))
3113
3114	call WRITEDIAGFI(ngrid,'ccnNco2_meteor','CCNco2_meteor number',
3115	& 'part.kg-1',3,zq(:,:,igcm_ccnco2_meteor_number))
3116	end if
3117
3118	end if ! of if (co2clouds)
3119	end if ! of if (tracer.and.(igcm_co2.ne.0))
3120
3121	! Output He tracer, if there is one
3122	if (tracer.and.(igcm_he.ne.0)) then
3123	call WRITEDIAGFI(ngrid,"he","helium mass mixing ratio",
3124	& "kg/kg",3,zq(1,1,igcm_he))
3125	vmr = zq(1:ngrid,1:nlayer,igcm_he)
3126	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_he)
3127	call WRITEDIAGFI(ngrid,'vmr_he','helium vol. mixing ratio',
3128	& 'mol/mol',3,vmr)
3129	end if
3130
3131	c ----------------------------------------------------------
3132	c Outputs of the water cycle
3133	c ----------------------------------------------------------
3134	if (tracer) then
3135	if (water) then
3136	#ifdef MESOINI
3137	!!!! waterice = q01, voir readmeteo.F90
3138	call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice),
3139	& 'kg/kg',3,
3140	& zq(1:ngrid,1:nlayer,igcm_h2o_ice))
3141	!!!! watervapor = q02, voir readmeteo.F90
3142	call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap),
3143	& 'kg/kg',3,
3144	& zq(1:ngrid,1:nlayer,igcm_h2o_vap))
3145	!!!! surface waterice qsurf02 (voir readmeteo)
3146	call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
3147	& 'kg.m-2',2,
3148	& qsurf(1:ngrid,igcm_h2o_ice))
3149	#endif
3150	call WRITEDIAGFI(ngrid,'mtot',
3151	& 'total mass of water vapor',
3152	& 'kg/m2',2,mtot)
3153	call WRITEDIAGFI(ngrid,'icetot',
3154	& 'total mass of water ice',
3155	& 'kg/m2',2,icetot)
3156	vmr = zq(1:ngrid,1:nlayer,igcm_h2o_ice)
3157	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
3158	call WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
3159	& 'mol/mol',3,vmr)
3160	vmr = zq(1:ngrid,1:nlayer,igcm_h2o_vap)
3161	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
3162	call WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
3163	& 'mol/mol',3,vmr)
3164	call WRITEDIAGFI(ngrid,'reffice',
3165	& 'Mean reff',
3166	& 'm',2,rave)
3167	call WRITEDIAGFI(ngrid,'h2o_ice','h2o_ice','kg/kg',
3168	& 3,zq(:,:,igcm_h2o_ice))
3169	call WRITEDIAGFI(ngrid,'h2o_vap','h2o_vap','kg/kg',
3170	& 3,zq(:,:,igcm_h2o_vap))
3171
3172	if (hdo) then
3173	vmr=zq(1:ngrid,1:nlayer,igcm_hdo_ice)
3174	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_ice)
3175	CALL WRITEDIAGFI(ngrid,'vmr_hdoice','hdo ice vmr',
3176	& 'mol/mol',3,vmr)
3177	vmr=zq(1:ngrid,1:nlayer,igcm_hdo_vap)
3178	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_vap)
3179	CALL WRITEDIAGFI(ngrid,'vmr_hdovap','hdo vap vmr',
3180	& 'mol/mol',3,vmr)
3181	call WRITEDIAGFI(ngrid,'hdo_ice','hdo_ice','kg/kg',
3182	& 3,zq(:,:,igcm_hdo_ice))
3183	call WRITEDIAGFI(ngrid,'hdo_vap','hdo_vap','kg/kg',
3184	& 3,zq(:,:,igcm_hdo_vap))
3185
3186	CALL WRITEDIAGFI(ngrid,'mtotD',
3187	& 'total mass of HDO vapor',
3188	& 'kg/m2',2,mtotD)
3189	CALL WRITEDIAGFI(ngrid,'icetotD',
3190	& 'total mass of HDO ice',
3191	& 'kg/m2',2,icetotD)
3192
3193	C Calculation of the D/H ratio
3194	do l=1,nlayer
3195	do ig=1,ngrid
3196	if (zq(ig,l,igcm_h2o_vap).gt.qperemin) then
3197	DoH_vap(ig,l) = ( zq(ig,l,igcm_hdo_vap)/
3198	& zq(ig,l,igcm_h2o_vap) )1./(2.155.76e-6)
3199	else
3200	DoH_vap(ig,l) = 0.
3201	endif
3202	enddo
3203	enddo
3204
3205	do l=1,nlayer
3206	do ig=1,ngrid
3207	if (zq(ig,l,igcm_h2o_ice).gt.qperemin) then
3208	DoH_ice(ig,l) = ( zq(ig,l,igcm_hdo_ice)/
3209	& zq(ig,l,igcm_h2o_ice) )/(2.*155.76e-6)
3210	else
3211	DoH_ice(ig,l) = 0.
3212	endif
3213	enddo
3214	enddo
3215
3216	CALL WRITEDIAGFI(ngrid,'DoH_vap',
3217	& 'D/H ratio in vapor',
3218	& ' ',3,DoH_vap)
3219	CALL WRITEDIAGFI(ngrid,'DoH_ice',
3220	& 'D/H ratio in ice',
3221	& '',3,DoH_ice)
3222
3223	endif !hdo
3224
3225	!A. Pottier
3226	CALL WRITEDIAGFI(ngrid,'rmoym',
3227	& 'alternative reffice',
3228	& 'm',2,rave2)
3229	call WRITEDIAGFI(ngrid,'saturation',
3230	& 'h2o vap saturation ratio','dimless',3,satu)
3231	if (scavenging) then
3232	CALL WRITEDIAGFI(ngrid,"Nccntot",
3233	& "condensation nuclei","Nbr/m2",
3234	& 2,Nccntot)
3235	CALL WRITEDIAGFI(ngrid,"Mccntot",
3236	& "mass condensation nuclei","kg/m2",
3237	& 2,Mccntot)
3238	endif
3239	call WRITEDIAGFI(ngrid,'rice','Ice particle size',
3240	& 'm',3,rice)
3241	call WRITEDIAGFI(ngrid,'h2o_ice_s',
3242	& 'surface h2o_ice',
3243	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
3244	if (hdo) then
3245	call WRITEDIAGFI(ngrid,'hdo_ice_s',
3246	& 'surface hdo_ice',
3247	& 'kg.m-2',2,qsurf(1,igcm_hdo_ice))
3248
3249	do ig=1,ngrid
3250	if (qsurf(ig,igcm_h2o_ice).gt.qperemin) then
3251	DoH_surf(ig) = 0.5*( qsurf(ig,igcm_hdo_ice)/
3252	& qsurf(ig,igcm_h2o_ice) )/155.76e-6
3253	else
3254	DoH_surf(ig) = 0.
3255	endif
3256	enddo
3257
3258	call WRITEDIAGFI(ngrid,'DoH_surf',
3259	& 'surface D/H',
3260	& '',2,DoH_surf)
3261	endif ! hdo
3262
3263	CALL WRITEDIAGFI(ngrid,'albedo',
3264	& 'albedo',
3265	& '',2,albedo(1,1))
3266	if (tifeedback) then
3267	call WRITEDIAGSOIL(ngrid,"soiltemp",
3268	& "Soil temperature","K",
3269	& 3,tsoil)
3270	call WRITEDIAGSOIL(ngrid,'soilti',
3271	& 'Soil Thermal Inertia',
3272	& 'J.s-1/2.m-2.K-1',3,inertiesoil)
3273	endif
3274	!A. Pottier
3275	if (CLFvarying) then !AP14 nebulosity
3276	call WRITEDIAGFI(ngrid,'totcloudfrac',
3277	& 'Total cloud fraction',
3278	& ' ',2,totcloudfrac)
3279	end if !clf varying
3280	end if !(water)
3281
3282	if (water.and..not.photochem) then
3283	iq = nq
3284	c write(str2(1:2),'(i2.2)') iq
3285	c call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
3286	c & 'kg.m-2',2,zdqscloud(1,iq))
3287	c call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
3288	c & 'kg/kg',3,zdqchim(1,1,iq))
3289	c call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
3290	c & 'kg/kg',3,zdqdif(1,1,iq))
3291	c call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
3292	c & 'kg/kg',3,zdqadj(1,1,iq))
3293	c call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
3294	c & 'kg/kg',3,zdqc(1,1,iq))
3295	end if !(water.and..not.photochem)
3296	end if !tracer
3297
3298	c ----------------------------------------------------------
3299	c Outputs of the dust cycle
3300	c ----------------------------------------------------------
3301
3302	call WRITEDIAGFI(ngrid,'tau_pref_scenario',
3303	& 'Prescribed visible dust optical depth at 610Pa',
3304	& 'NU',2,tau_pref_scenario)
3305
3306	call WRITEDIAGFI(ngrid,'tau_pref_gcm',
3307	& 'Visible dust optical depth at 610Pa in the GCM',
3308	& 'NU',2,tau_pref_gcm)
3309
3310	if (reff_driven_IRtoVIS_scenario) then
3311	call WRITEDIAGFI(ngrid,'IRtoVIScoef',
3312	& 'conversion coeff for dust tau from abs9.3um to ext0.67um',
3313	& '/',2,IRtoVIScoef)
3314	endif
3315
3316	if (tracer.and.(dustbin.ne.0)) then
3317
3318	call WRITEDIAGFI(ngrid,'tau','taudust','SI',2,tau(1,1))
3319
3320	#ifndef MESOINI
3321	if (doubleq) then
3322	c call WRITEDIAGFI(ngrid,'qsurf','qsurf',
3323	c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))
3324	c call WRITEDIAGFI(ngrid,'Nsurf','N particles',
3325	c & 'N.m-2',2,qsurf(1,igcm_dust_number))
3326	c call WRITEDIAGFI(ngrid,'dqsdev','ddevil lift',
3327	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
3328	c call WRITEDIAGFI(ngrid,'dqssed','sedimentation',
3329	c & 'kg.m-2.s-1',2,zdqssed(1,1))
3330	c call WRITEDIAGFI(ngrid,'dqsdif','diffusion',
3331	c & 'kg.m-2.s-1',2,zdqsdif(1,1))
3332	c call WRITEDIAGFI(ngrid,'sedice','sedimented ice',
3333	c & 'kg.m-2.s-1',2,zdqssed(:,igcm_h2o_ice))
3334	c call WRITEDIAGFI(ngrid,'subice','sublimated ice',
3335	c & 'kg.m-2.s-1',2,zdqsdif(:,igcm_h2o_ice))
3336	call WRITEDIAGFI(ngrid,'dqsdust',
3337	& 'deposited surface dust mass',
3338	& 'kg.m-2.s-1',2,dqdustsurf)
3339	call WRITEDIAGFI(ngrid,'dqndust',
3340	& 'deposited surface dust number',
3341	& 'number.m-2.s-1',2,dndustsurf)
3342	call WRITEDIAGFI(ngrid,'reffdust','reffdust',
3343	& 'm',3,rdust*ref_r0)
3344	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
3345	& 'kg/kg',3,qdust)
3346	call WRITEDIAGFI(ngrid,'dustN','Dust number',
3347	& 'part/kg',3,ndust)
3348
3349	select case (trim(dustiropacity))
3350	case ("tes")
3351	call WRITEDIAGFI(ngrid,'dsodust',
3352	& 'density scaled extinction opacity of std dust at 9.3um(TES)',
3353	& 'm2.kg-1',3,dsodust)
3354	call WRITEDIAGFI(ngrid,'dso',
3355	& 'density scaled extinction opacity of all dust at 9.3um(TES)',
3356	& 'm2.kg-1',3,dsodust+dsords+dsotop)
3357	case ("mcs")
3358	call WRITEDIAGFI(ngrid,'dsodust',
3359	& 'density scaled extinction opacity of std dust at 21.6um(MCS)',
3360	& 'm2.kg-1',3,dsodust)
3361	call WRITEDIAGFI(ngrid,'dso',
3362	& 'density scaled extinction opacity of all dust at 21.6um(MCS)',
3363	& 'm2.kg-1',3,dsodust+dsords+dsotop)
3364	end select
3365	else ! (doubleq=.false.)
3366	do iq=1,dustbin
3367	write(str2(1:2),'(i2.2)') iq
3368	call WRITEDIAGFI(ngrid,'q'//str2,'mix. ratio',
3369	& 'kg/kg',3,zq(1,1,iq))
3370	call WRITEDIAGFI(ngrid,'qsurf'//str2,'qsurf',
3371	& 'kg.m-2',2,qsurf(1,iq))
3372	end do
3373	endif ! (doubleq)
3374
3375	if (rdstorm) then ! writediagfi tendencies stormdust tracers
3376	call WRITEDIAGFI(ngrid,'reffstormdust','reffstormdust',
3377	& 'm',3,rstormdust*ref_r0)
3378	call WRITEDIAGFI(ngrid,'mstormdtot',
3379	& 'total mass of stormdust only',
3380	& 'kg.m-2',2,mstormdtot)
3381	call WRITEDIAGFI(ngrid,'mdusttot',
3382	& 'total mass of dust only',
3383	& 'kg.m-2',2,mdusttot)
3384	call WRITEDIAGFI(ngrid,'rdsdqsdust',
3385	& 'deposited surface stormdust mass',
3386	& 'kg.m-2.s-1',2,rdsdqdustsurf)
3387	call WRITEDIAGFI(ngrid,'rdsdustq','storm Dust mass mr',
3388	& 'kg/kg',3,rdsqdust)
3389	call WRITEDIAGFI(ngrid,'rdsdustqmodel','storm Dust massmr',
3390	& 'kg/kg',3,pq(:,:,igcm_stormdust_mass))
3391	call WRITEDIAGFI(ngrid,'rdsdustN','storm Dust number',
3392	& 'part/kg',3,rdsndust)
3393	call WRITEDIAGFI(ngrid,"stormfract",
3394	& "fraction of the mesh, with stormdust","none",
3395	& 2,totstormfract)
3396	call WRITEDIAGFI(ngrid,'qsurf',
3397	& 'stormdust injection',
3398	& 'kg.m-2',2,qsurf(:,igcm_stormdust_mass))
3399	call WRITEDIAGFI(ngrid,'pdqsurf',
3400	& 'tendancy stormdust mass at surface',
3401	& 'kg.m-2',2,dqsurf(:,igcm_stormdust_mass))
3402	call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust',
3403	& 'm/s',3,wspeed(:,1:nlayer))
3404	call WRITEDIAGFI(ngrid,'zdqsed_dustq'
3405	& ,'sedimentation q','kg.m-2.s-1',3,
3406	& zdqsed(:,:,igcm_dust_mass))
3407	call WRITEDIAGFI(ngrid,'zdqssed_dustq'
3408	& ,'sedimentation q','kg.m-2.s-1',2,
3409	& zdqssed(:,igcm_dust_mass))
3410	call WRITEDIAGFI(ngrid,'zdqsed_rdsq'
3411	& ,'sedimentation q','kg.m-2.s-1',3,
3412	& zdqsed(:,:,igcm_stormdust_mass))
3413	call WRITEDIAGFI(ngrid,'rdust','rdust',
3414	& 'm',3,rdust)
3415	call WRITEDIAGFI(ngrid,'rstormdust','rstormdust',
3416	& 'm',3,rstormdust)
3417	call WRITEDIAGFI(ngrid,'totaldustq','total dust mass',
3418	& 'kg/kg',3,qdusttotal)
3419
3420	select case (trim(dustiropacity))
3421	case ("tes")
3422	call WRITEDIAGFI(ngrid,'dsords',
3423	& 'density scaled extinction opacity of stormdust at 9.3um(TES)',
3424	& 'm2.kg-1',3,dsords)
3425	case ("mcs")
3426	call WRITEDIAGFI(ngrid,'dsords',
3427	& 'density scaled extinction opacity of stormdust at 21.6um(MCS)',
3428	& 'm2.kg-1',3,dsords)
3429	end select
3430	endif ! (rdstorm)
3431
3432	if (topflows) then
3433	call WRITEDIAGFI(ngrid,'refftopdust','refftopdust',
3434	& 'm',3,rtopdust*ref_r0)
3435	call WRITEDIAGFI(ngrid,'topdustq','top Dust mass mr',
3436	& 'kg/kg',3,pq(:,:,igcm_topdust_mass))
3437	call WRITEDIAGFI(ngrid,'topdustN','top Dust number',
3438	& 'part/kg',3,pq(:,:,igcm_topdust_number))
3439	select case (trim(dustiropacity))
3440	case ("tes")
3441	call WRITEDIAGFI(ngrid,'dsotop',
3442	& 'density scaled extinction opacity of topdust at 9.3um(TES)',
3443	& 'm2.kg-1',3,dsotop)
3444	case ("mcs")
3445	call WRITEDIAGFI(ngrid,'dsotop',
3446	& 'density scaled extinction opacity of topdust at 21.6um(MCS)',
3447	& 'm2.kg-1',3,dsotop)
3448	end select
3449	endif ! (topflows)
3450
3451	if (dustscaling_mode==2) then
3452	call writediagfi(ngrid,"dust_rad_adjust",
3453	& "radiative adjustment coefficient for dust",
3454	& "",2,dust_rad_adjust)
3455	endif
3456
3457	if (scavenging) then
3458	call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr',
3459	& 'kg/kg',3,qccn)
3460	call WRITEDIAGFI(ngrid,'ccnN','CCN number',
3461	& 'part/kg',3,nccn)
3462	call WRITEDIAGFI(ngrid,'surfccnq','Surf nuclei mass mr',
3463	& 'kg.m-2',2,qsurf(1,igcm_ccn_mass))
3464	call WRITEDIAGFI(ngrid,'surfccnN','Surf nuclei number',
3465	& 'kg.m-2',2,qsurf(1,igcm_ccn_number))
3466	endif ! (scavenging)
3467
3468	c if (submicron) then
3469	c call WRITEDIAGFI(ngrid,'dustsubm','subm mass mr',
3470	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
3471	c endif ! (submicron)
3472
3473	#else
3474	! !!! to initialize mesoscale we need scaled variables
3475	! !!! because this must correspond to starting point for tracers
3476	! call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
3477	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass))
3478	! call WRITEDIAGFI(ngrid,'dustN','Dust number',
3479	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number))
3480	! call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
3481	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass))
3482	! call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
3483	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number))
3484	if (freedust) then
3485	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
3486	& 'kg/kg',3,qdust)
3487	call WRITEDIAGFI(ngrid,'dustN','Dust number',
3488	& 'part/kg',3,ndust)
3489	call WRITEDIAGFI(ngrid,'ccn','CCN mass mr',
3490	& 'kg/kg',3,qccn)
3491	call WRITEDIAGFI(ngrid,'ccnN','CCN number',
3492	& 'part/kg',3,nccn)
3493	else
3494	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
3495	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
3496	call WRITEDIAGFI(ngrid,'dustN','Dust number',
3497	& 'part/kg',3,pq(1,1,igcm_dust_number))
3498	call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
3499	& 'kg/kg',3,pq(1,1,igcm_ccn_mass))
3500	call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
3501	& 'part/kg',3,pq(1,1,igcm_ccn_number))
3502	endif
3503	#endif
3504
3505	end if ! (tracer.and.(dustbin.ne.0))
3506
3507	c ----------------------------------------------------------
3508	c GW non-oro outputs
3509	c ----------------------------------------------------------
3510
3511	if(calllott_nonoro) then
3512	call WRITEDIAGFI(ngrid,"dugwno","GW non-oro dU","m/s2",
3513	$ 3,d_u_hin/ptimestep)
3514	call WRITEDIAGFI(ngrid,"dvgwno","GW non-oro dV","m/s2",
3515	$ 3,d_v_hin/ptimestep)
3516	endif !(calllott_nonoro)
3517
3518	c ----------------------------------------------------------
3519	c Thermospheric outputs
3520	c ----------------------------------------------------------
3521
3522	if(callthermos) then
3523
3524	call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s",
3525	$ 3,zdtnlte)
3526	call WRITEDIAGFI(ngrid,"quv","UV heating","K/s",
3527	$ 3,zdteuv)
3528	call WRITEDIAGFI(ngrid,"cond","Thermal conduction","K/s",
3529	$ 3,zdtconduc)
3530	call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s",
3531	$ 3,zdtnirco2)
3532
3533	!H, H2 and D escape fluxes
3534
3535	call WRITEDIAGFI(ngrid,"PhiH","H escape flux","s-1",
3536	$ 0,[PhiEscH])
3537	call WRITEDIAGFI(ngrid,"PhiH2","H2 escape flux","s-1",
3538	$ 0,[PhiEscH2])
3539	call WRITEDIAGFI(ngrid,"PhiD","D escape flux","s-1",
3540	$ 0,[PhiEscD])
3541
3542	! call wstats(ngrid,"PhiH","H escape flux","s-1",
3543	! $ 0,[PhiEscH])
3544	! call wstats(ngrid,"PhiH2","H2 escape flux","s-1",
3545	! $ 0,[PhiEscH2])
3546	! call wstats(ngrid,"PhiD","D escape flux","s-1",
3547	! $ 0,[PhiEscD])
3548
3549	! call wstats(ngrid,"q15um","15 um cooling","K/s",
3550	! $ 3,zdtnlte)
3551	! call wstats(ngrid,"quv","UV heating","K/s",
3552	! $ 3,zdteuv)
3553	! call wstats(ngrid,"cond","Thermal conduction","K/s",
3554	! $ 3,zdtconduc)
3555	! call wstats(ngrid,"qnir","NIR heating","K/s",
3556	! $ 3,zdtnirco2)
3557
3558	endif !(callthermos)
3559
3560	call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s",
3561	$ 3,zdtnlte)
3562	call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s",
3563	$ 3,zdtnirco2)
3564
3565	c ----------------------------------------------------------
3566	c ----------------------------------------------------------
3567	c PBL OUTPUS
3568	c ----------------------------------------------------------
3569	c ----------------------------------------------------------
3570
3571	c ----------------------------------------------------------
3572	c Outputs of thermals
3573	c ----------------------------------------------------------
3574	if (calltherm) then
3575	! call WRITEDIAGFI(ngrid,'dtke',
3576	! & 'tendance tke thermiques','m2/s2',
3577	! & 3,dtke_th)
3578	! call WRITEDIAGFI(ngrid,'d_u_ajs',
3579	! & 'tendance u thermiques','m/s',
3580	! & 3,pdu_th*ptimestep)
3581	! call WRITEDIAGFI(ngrid,'d_v_ajs',
3582	! & 'tendance v thermiques','m/s',
3583	! & 3,pdv_th*ptimestep)
3584	! if (tracer) then
3585	! if (nq .eq. 2) then
3586	! call WRITEDIAGFI(ngrid,'deltaq_th',
3587	! & 'delta q thermiques','kg/kg',
3588	! & 3,ptimestep*pdq_th(:,:,2))
3589	! end if
3590	! end if
3591
3592	call WRITEDIAGFI(ngrid,'zmax_th',
3593	& 'hauteur du thermique','m',
3594	& 2,zmax_th)
3595	call WRITEDIAGFI(ngrid,'hfmax_th',
3596	& 'maximum TH heat flux','K.m/s',
3597	& 2,hfmax_th)
3598	call WRITEDIAGFI(ngrid,'wstar',
3599	& 'maximum TH vertical velocity','m/s',
3600	& 2,wstar)
3601	end if
3602
3603	c ----------------------------------------------------------
3604	c ----------------------------------------------------------
3605	c END OF PBL OUTPUS
3606	c ----------------------------------------------------------
3607	c ----------------------------------------------------------
3608
3609
3610	c ----------------------------------------------------------
3611	c Output in netcdf file "diagsoil.nc" for subterranean
3612	c variables (output every "ecritphy", as for writediagfi)
3613	c ----------------------------------------------------------
3614
3615	! Write soil temperature
3616	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
3617	! & 3,tsoil)
3618	! Write surface temperature
3619	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
3620	! & 2,tsurf)
3621
3622	c ==========================================================
3623	c END OF WRITEDIAGFI
3624	c ==========================================================
3625	#endif
3626	! of ifdef MESOSCALE
3627
3628	ELSE ! if(ngrid.eq.1)
3629
3630	#ifndef MESOSCALE
3631	write(*,
3632	& '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)')
3633	& zls*180./pi,odpref,tau_pref_scenario
3634	c ----------------------------------------------------------------------
3635	c Output in grads file "g1d" (ONLY when using testphys1d)
3636	c (output at every X physical timestep)
3637	c ----------------------------------------------------------------------
3638	c
3639	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
3640	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
3641	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
3642	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
3643	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
3644	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
3645	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
3646
3647	! THERMALS STUFF 1D
3648	if(calltherm) then
3649
3650	call WRITEDIAGFI(ngrid,'lmax_th',
3651	& 'hauteur du thermique','point',
3652	& 0,lmax_th_out)
3653	call WRITEDIAGFI(ngrid,'zmax_th',
3654	& 'hauteur du thermique','m',
3655	& 0,zmax_th)
3656	call WRITEDIAGFI(ngrid,'hfmax_th',
3657	& 'maximum TH heat flux','K.m/s',
3658	& 0,hfmax_th)
3659	call WRITEDIAGFI(ngrid,'wstar',
3660	& 'maximum TH vertical velocity','m/s',
3661	& 0,wstar)
3662
3663	end if ! of if (calltherm)
3664
3665	call WRITEDIAGFI(ngrid,'w','vertical velocity'
3666	& ,'m/s',1,pw)
3667	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
3668	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
3669	& tsurf)
3670	call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
3671	call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
3672
3673	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
3674	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
3675	call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
3676	call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate",
3677	& "K.s-1",1,dtrad)
3678	call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
3679	& 'w.m-2',1,zdtsw)
3680	call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
3681	& 'w.m-2',1,zdtlw)
3682	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
3683	& ,"kg.m-2",0,co2ice)
3684
3685	if (igcm_co2.ne.0) then
3686	call co2sat(ngrid*nlayer,zt,zqsatco2)
3687	do ig=1,ngrid
3688	do l=1,nlayer
3689	satuco2(ig,l) = zq(ig,l,igcm_co2)*
3690	& (mmean(ig,l)/44.01)*zplay(ig,l)/zqsatco2(ig,l)
3691
3692	c write(,) "In PHYSIQMOD, pt,zt,time ",pt(ig,l)
3693	c & ,zt(ig,l),ptime
3694	enddo
3695	enddo
3696	endif
3697
3698	call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
3699	call WRITEDIAGFI(ngrid,'temp','Temperature ',
3700	& 'K JA',1,zt)
3701	c call WRITEDIAGFI(ngrid,'temp2','Temperature ',
3702	c & 'K JA2',1,pt)
3703
3704	if(tracer) then
3705	c CALL writeg1d(ngrid,1,tau,'tau','SI')
3706	do iq=1,nq
3707	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
3708	call WRITEDIAGFI(ngrid,trim(noms(iq)),
3709	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
3710	end do
3711	if (doubleq) then
3712	call WRITEDIAGFI(ngrid,'rdust','rdust',
3713	& 'm',1,rdust)
3714	endif ! doubleq 1D
3715	if (rdstorm) then
3716	call writediagfi(1,'aerosol_dust','opacity of env. dust',''
3717	& ,1,aerosol(:,:,iaer_dust_doubleq))
3718	call writediagfi(1,'aerosol_stormdust',
3719	& 'opacity of storm dust',''
3720	& ,1,aerosol(:,:,iaer_stormdust_doubleq))
3721	call WRITEDIAGFI(ngrid,'dqsdifdustq','diffusion',
3722	& 'kg.m-2.s-1',0,zdqsdif(1,igcm_dust_mass))
3723	call WRITEDIAGFI(ngrid,'dqsdifrdsq','diffusion',
3724	& 'kg.m-2.s-1',0,zdqsdif(1,igcm_stormdust_mass))
3725	call WRITEDIAGFI(ngrid,'mstormdtot',
3726	& 'total mass of stormdust only',
3727	& 'kg.m-2',0,mstormdtot)
3728	call WRITEDIAGFI(ngrid,'mdusttot',
3729	& 'total mass of dust only',
3730	& 'kg.m-2',0,mdusttot)
3731	call WRITEDIAGFI(ngrid,'tau_pref_scenario',
3732	& 'Prescribed dust ref opt depth at 610 Pa',
3733	& 'NU',0,tau_pref_scenario)
3734	call WRITEDIAGFI(ngrid,'tau_pref_gcm',
3735	& 'Dust ref opt depth at 610 Pa in the GCM',
3736	& 'NU',0,tau_pref_gcm)
3737	call WRITEDIAGFI(ngrid,'rdsdqsdust',
3738	& 'deposited surface stormdust mass',
3739	& 'kg.m-2.s-1',0,rdsdqdustsurf)
3740	call WRITEDIAGFI(ngrid,'rdsdustq','storm Dust mass mr',
3741	& 'kg/kg',1,rdsqdust)
3742	call WRITEDIAGFI(ngrid,"stormfract",
3743	& "fraction of the mesh,with stormdust",
3744	& "none",0,totstormfract)
3745	call WRITEDIAGFI(ngrid,'rdsqsurf',
3746	& 'stormdust at surface',
3747	& 'kg.m-2',0,qsurf(:,igcm_stormdust_mass))
3748	call WRITEDIAGFI(ngrid,'qsurf',
3749	& 'dust mass at surface',
3750	& 'kg.m-2',0,qsurf(:,igcm_dust_mass))
3751	call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust',
3752	& 'm/s',1,wspeed)
3753	call WRITEDIAGFI(ngrid,'totaldustq','total dust mass',
3754	& 'kg/kg',1,qdusttotal)
3755	call WRITEDIAGFI(ngrid,'dsords',
3756	& 'density scaled opacity of stormdust',
3757	& 'm2.kg-1',1,dsords)
3758	call WRITEDIAGFI(ngrid,'zdqsed_dustq'
3759	& ,'sedimentation q','kg.m-2.s-1',1,
3760	& zdqsed(1,:,igcm_dust_mass))
3761	call WRITEDIAGFI(ngrid,'zdqsed_rdsq'
3762	& ,'sedimentation q','kg.m-2.s-1',1,
3763	& zdqsed(1,:,igcm_stormdust_mass))
3764	endif !(rdstorm 1D)
3765
3766	if (water.AND.tifeedback) then
3767	call WRITEDIAGFI(ngrid,"soiltemp",
3768	& "Soil temperature","K",
3769	& 1,tsoil)
3770	call WRITEDIAGFI(ngrid,'soilti',
3771	& 'Soil Thermal Inertia',
3772	& 'J.s-1/2.m-2.K-1',1,inertiesoil)
3773	endif
3774	end if
3775
3776	cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc
3777	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3778	IF (water) THEN
3779
3780	if (.not.activice) then
3781
3782	tauTES=0
3783	do l=1,nlayer
3784	Qabsice = min(
3785	& max(0.4e6rice(1,l)(1.+nuice_ref)-0.05 ,0.),1.2
3786	& )
3787	opTES(1,l)= 0.75 * Qabsice *
3788	& zq(1,l,igcm_h2o_ice) *
3789	& (zplev(1,l) - zplev(1,l+1)) / g
3790	& / (rho_ice * rice(1,l) * (1.+nuice_ref))
3791	tauTES=tauTES+ opTES(1,l)
3792	enddo
3793	CALL WRITEDIAGFI(ngrid,'tauTESap',
3794	& 'tau abs 825 cm-1',
3795	& '',0,tauTES)
3796	else
3797
3798	CALL WRITEDIAGFI(ngrid,'tauTES',
3799	& 'tau abs 825 cm-1',
3800	& '',0,taucloudtes)
3801	endif
3802
3803	mtot = 0
3804	icetot = 0
3805	h2otot = qsurf(1,igcm_h2o_ice)
3806	if (hdo) THEN
3807	mtotD = 0
3808	icetotD = 0
3809	hdotot = qsurf(1,igcm_hdo_ice)
3810	ENDIF !hdo
3811
3812	do l=1,nlayer
3813	mtot = mtot + zq(1,l,igcm_h2o_vap)
3814	& * (zplev(1,l) - zplev(1,l+1)) / g
3815	icetot = icetot + zq(1,l,igcm_h2o_ice)
3816	& * (zplev(1,l) - zplev(1,l+1)) / g
3817	if (hdo) THEN
3818	mtotD = mtotD + zq(1,l,igcm_hdo_vap)
3819	& * (zplev(1,l) - zplev(1,l+1)) / g
3820	icetotD = icetotD + zq(1,l,igcm_hdo_ice)
3821	& * (zplev(1,l) - zplev(1,l+1)) / g
3822	ENDIF !hdo
3823	end do
3824	h2otot = h2otot+mtot+icetot
3825	IF (hdo) then
3826	hdotot = hdotot+mtotD+icetotD
3827	ENDIF ! hdo
3828
3829
3830	CALL WRITEDIAGFI(ngrid,'h2otot',
3831	& 'h2otot',
3832	& 'kg/m2',0,h2otot)
3833	CALL WRITEDIAGFI(ngrid,'mtot',
3834	& 'mtot',
3835	& 'kg/m2',0,mtot)
3836	CALL WRITEDIAGFI(ngrid,'icetot',
3837	& 'icetot',
3838	& 'kg/m2',0,icetot)
3839
3840	IF (hdo) THEN
3841	CALL WRITEDIAGFI(ngrid,'mtotD',
3842	& 'mtotD',
3843	& 'kg/m2',0,mtotD)
3844	CALL WRITEDIAGFI(ngrid,'icetotD',
3845	& 'icetotD',
3846	& 'kg/m2',0,icetotD)
3847	CALL WRITEDIAGFI(ngrid,'hdotot',
3848	& 'hdotot',
3849	& 'kg/m2',0,hdotot)
3850
3851	C Calculation of the D/H ratio
3852	do l=1,nlayer
3853	if (zq(1,l,igcm_h2o_vap).gt.qperemin) then
3854	DoH_vap(1,l) = 0.5*( zq(1,l,igcm_hdo_vap)/
3855	& zq(1,l,igcm_h2o_vap) )/155.76e-6
3856	else
3857	DoH_vap(1,l) = 0.
3858	endif
3859	enddo
3860
3861	do l=1,nlayer
3862	if (zq(1,l,igcm_h2o_ice).gt.qperemin) then
3863	DoH_ice(1,l) = 0.5*( zq(1,l,igcm_hdo_ice)/
3864	& zq(1,l,igcm_h2o_ice) )/155.76e-6
3865	else
3866	DoH_ice(1,l) = 0.
3867	endif
3868	enddo
3869
3870	CALL WRITEDIAGFI(ngrid,'DoH_vap',
3871	& 'D/H ratio in vapor',
3872	& ' ',1,DoH_vap)
3873	CALL WRITEDIAGFI(ngrid,'DoH_ice',
3874	& 'D/H ratio in ice',
3875	& '',1,DoH_ice)
3876
3877	ENDIF !Hdo
3878
3879
3880	if (scavenging) then
3881
3882	rave = 0
3883	do l=1,nlayer
3884	cccc Column integrated effective ice radius
3885	cccc is weighted by total ice surface area (BETTER)
3886	rave = rave + tauscaling(1) *
3887	& zq(1,l,igcm_ccn_number) *
3888	& (zplev(1,l) - zplev(1,l+1)) / g *
3889	& rice(1,l) * rice(1,l)* (1.+nuice_ref)
3890	enddo
3891	rave=icetot0.75/max(ravepi*rho_ice,1.e-30) ! surface weight
3892
3893	Nccntot= 0
3894	call watersat(ngrid*nlayer,zt,zplay,zqsat)
3895	do l=1,nlayer
3896	Nccntot = Nccntot +
3897	& zq(1,l,igcm_ccn_number)*tauscaling(1)
3898	& *(zplev(1,l) - zplev(1,l+1)) / g
3899	satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l)
3900	satu(1,l) = (max(satu(1,l),float(1))-1)
3901	! & * zq(1,l,igcm_h2o_vap) *
3902	! & (zplev(1,l) - zplev(1,l+1)) / g
3903	enddo
3904	call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1,
3905	& satu)
3906	CALL WRITEDIAGFI(ngrid,'Nccntot',
3907	& 'Nccntot',
3908	& 'nbr/m2',0,Nccntot)
3909
3910	call WRITEDIAGFI(ngrid,'zdqsed_dustq'
3911	& ,'sedimentation q','kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass))
3912	call WRITEDIAGFI(ngrid,'zdqsed_dustN'
3913	&,'sedimentation N','Nbr.m-2.s-1',1,
3914	& zdqsed(1,:,igcm_dust_number))
3915
3916	else ! of if (scavenging)
3917
3918	cccc Column integrated effective ice radius
3919	cccc is weighted by total ice mass (LESS GOOD)
3920	rave = 0
3921	do l=1,nlayer
3922	rave = rave + zq(1,l,igcm_h2o_ice)
3923	& * (zplev(1,l) - zplev(1,l+1)) / g
3924	& * rice(1,l) * (1.+nuice_ref)
3925	enddo
3926	rave=max(rave/max(icetot,1.e-30),1.e-30) ! mass weight
3927	endif ! of if (scavenging)
3928
3929
3930	CALL WRITEDIAGFI(ngrid,'reffice',
3931	& 'reffice',
3932	& 'm',0,rave)
3933
3934	!Alternative A. Pottier weighting
3935	rave2 = 0.
3936	totrave2 = 0.
3937	do l=1,nlayer
3938	rave2 =rave2+ zq(1,l,igcm_h2o_ice)*rice(1,l)
3939	totrave2 = totrave2 + zq(1,l,igcm_h2o_ice)
3940	end do
3941	rave2=max(rave2/max(totrave2,1.e-30),1.e-30)
3942	CALL WRITEDIAGFI(ngrid,'rmoym',
3943	& 'reffice',
3944	& 'm',0,rave2)
3945
3946	do iq=1,nq
3947	call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s',
3948	& trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq))
3949	end do
3950
3951	call WRITEDIAGFI(ngrid,"watercap","Water ice thickness"
3952	& ,"kg.m-2",0,watercap)
3953	call WRITEDIAGFI(ngrid,'zdqcloud_ice','cloud ice',
3954	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))
3955	call WRITEDIAGFI(ngrid,'zdqcloud_vap','cloud vap',
3956	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap))
3957	call WRITEDIAGFI(ngrid,'zdqcloud','cloud ice',
3958	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)
3959	& +zdqcloud(1,:,igcm_h2o_vap))
3960	IF (hdo) THEN
3961	call WRITEDIAGFI(ngrid,'zdqcloud_iceD','cloud ice hdo',
3962	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_hdo_ice))
3963	call WRITEDIAGFI(ngrid,'zdqcloud_vapD','cloud vap hdo',
3964	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_hdo_vap))
3965
3966	ENDIF ! hdo
3967
3968	call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
3969	& rice)
3970
3971	if (CLFvarying) then
3972	call WRITEDIAGFI(ngrid,'totcloudfrac',
3973	& 'Total cloud fraction',
3974	& ' ',0,totcloudfrac)
3975	endif !clfvarying
3976
3977	ENDIF ! of IF (water)
3978
3979	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3980	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3981
3982	zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g
3983
3984	do l=2,nlayer-1
3985	tmean=zt(1,l)
3986	if(zt(1,l).ne.zt(1,l-1))
3987	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
3988	zlocal(l)= zlocal(l-1)
3989	& -log(zplay(1,l)/zplay(1,l-1))rnew(1,l)tmean/g
3990	enddo
3991	zlocal(nlayer)= zlocal(nlayer-1)-
3992	& log(zplay(1,nlayer)/zplay(1,nlayer-1))*
3993	& rnew(1,nlayer)*tmean/g
3994	#endif
3995
3996	END IF ! if(ngrid.ne.1)
3997	! test for co2 conservation with co2 microphysics
3998	if (igcm_co2_ice.ne.0) then
3999	co2totB = 0. ! added by C.M.
4000	do ig=1,ngrid
4001	do l=1,nlayer
4002	co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
4003	& (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice)
4004	& +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
4005	enddo
4006	co2totB = co2totB + co2ice(ig)
4007	enddo
4008	else
4009	co2totB = 0. ! added by C.M.
4010	do ig=1,ngrid
4011	do l=1,nlayer
4012	co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
4013	& (pq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2)*ptimestep)
4014	enddo
4015	co2totB = co2totB + co2ice(ig)
4016	enddo
4017	endif ! of if (igcm_co2_ice.ne.0)
4018	co2conservation = (co2totA-co2totB)/co2totA
4019	call WRITEDIAGFI(ngrid, 'co2conservation',
4020	& 'Total CO2 mass conservation in physic',
4021	& ' ', 0, co2conservation)
4022	! XIOS outputs
4023	#ifdef CPP_XIOS
4024	! Send fields to XIOS: (NB these fields must also be defined as
4025	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
4026	CALL send_xios_field("ls",zls*180./pi)
4027
4028	CALL send_xios_field("controle",tab_cntrl_mod,1)
4029
4030	CALL send_xios_field("ap",ap,1)
4031	CALL send_xios_field("bp",bp,1)
4032	CALL send_xios_field("aps",aps,1)
4033	CALL send_xios_field("bps",bps,1)
4034
4035	CALL send_xios_field("phisinit",phisfi)
4036
4037	CALL send_xios_field("ps",ps)
4038	CALL send_xios_field("area",cell_area)
4039
4040	! CALL send_xios_field("ISR",fluxtop_dn_sw_tot)
4041	CALL send_xios_field("OLR",fluxtop_lw)
4042
4043	CALL send_xios_field("tsurf",tsurf)
4044	! CALL send_xios_field("inertiedat",inertiedat)
4045	CALL send_xios_field("tsoil",tsoil)
4046	CALL send_xios_field("co2ice",co2ice)
4047
4048	! CALL send_xios_field("temp",zt)
4049	CALL send_xios_field("u",zu)
4050	CALL send_xios_field("v",zv)
4051
4052	! CALL send_xios_field("rho",rho)
4053	! Orographic Gravity waves tendencies
4054	! if (calllott) then
4055	! CALL send_xios_field("dugw",zdugw/ptimestep)
4056	! CALL send_xios_field("dvgw",zdvgw/ptimestep)
4057	! CALL send_xios_field("dtgw",zdtgw/ptimestep)
4058	! endif
4059	!CREATE IF CO2CYCLE
4060	! if (tracer.and.(igcm_co2.ne.0)) then
4061	! CALL send_xios_field("co2",zq(:,:,igcm_co2))
4062	! endif
4063	! Water cycle
4064	! if (water) then
4065	! CALL send_xios_field("watercap",watercap)
4066	!CALL send_xios_field("watercaptag",watercaptag)
4067	! CALL send_xios_field("mtot",mtot)
4068	! CALL send_xios_field("icetot",icetot)
4069	! if (igcm_h2o_vap.ne.0 .and. igcm_h2o_ice.ne.0) then
4070	! CALL send_xios_field("h2o_vap",zq(:,:,igcm_h2o_vap))
4071	! CALL send_xios_field("h2o_ice",zq(:,:,igcm_h2o_ice))
4072	! endif
4073	! endif
4074	! if (.not.activice) then
4075	! CALL send_xios_field("tauTESap",tauTES)
4076	! else
4077	! CALL send_xios_field("tauTES",taucloudtes)
4078	! endif
4079
4080	! CALL send_xios_field("h2o_ice_s",qsurf(:,igcm_h2o_ice))
4081
4082
4083	if (lastcall.and.is_omp_master) then
4084	write(,) "physiq lastcall: call xios_context_finalize"
4085	call xios_context_finalize
4086	endif
4087	#endif
4088
4089	if (check_physics_outputs) then
4090	! Check the validity of updated fields at the end of the physics step
4091	call check_physics_fields("end of physiq:",zt,zu,zv,zplev,zq)
4092	endif
4093
4094	icount=icount+1
4095
4096	END SUBROUTINE physiq
4097
4098	END MODULE physiq_mod

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

Download in other formats:

Original Format