Context Navigation

physiq_mod.F @ 3016

Last change on this file since 3016 was 3016, checked in by emillour, 17 months ago
Mars PCM: Code cleanup: nltedata.h is only included in nltecool.F so turn this data into module data there. Also convert lteparams.h into module nlteparams_h.F90. And while at it also turn nlthermeq.F and blendrad.F into modules. EM
File size: 158.3 KB

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

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

Download in other formats:

Original Format