Context Navigation

physiq_mod.F @ 2880

Last change on this file since 2880 was 2880, checked in by llange, 22 months ago
MARS PCM Following R2879, remove stuff that were displayed on the screen LL
File size: 163.2 KB

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

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

Download in other formats:

Original Format