Context Navigation

physiq_mod.F @ 3242

Last change on this file since 3242 was 3230, checked in by llange, 17 months ago

Mars PCM

Move soil_tifeedback into a module waterice_tifeedback_mod.F90
The flag to call tifeedback has been changed from "tifeedback" to surfaceice_tifeedback for clarity
Add the possibility to change the thermal inertia while pore ices is forming in the soil: to do so, use the flag poreice_tifeedback. The computation is done in waterice_tifeedback_mod.F90.

For now, surfaceice_tifeedback and poreice_tifeedback can not be called together, we might think about how to merge later.
LL

File size: 162.4 KB

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

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

Download in other formats:

Original Format