Context Navigation

physiq_mod.F @ 4058

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

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

Download in other formats:

Original Format