Context Navigation

physiq_mod.F @ 3759

Last change on this file since 3759 was 3740, checked in by emillour, 10 months ago
Mars PCM: Code tidying; removed unused mufract.F (mucorr.F does the job); turned ambiguously named sig.F to sig_h2o.F90 and make it a module. EM
File size: 163.2 KB

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

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

Download in other formats:

Original Format