Context Navigation

physiq_mod.F @ 4022

Last change on this file since 4022 was 4022, checked in by aslmd, 9 days ago
MESOSCALE: sensible flux computation is broken, has to fix it by putting the right zcdh back in physiq_mod again
File size: 163.4 KB

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

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

Download in other formats:

Original Format