Context Navigation

physiq_mod.F @ 2997

Last change on this file since 2997 was 2965, checked in by abierjon, 20 months ago

Mars GCM:
Fixed r2963 which was preventing to compile the model without XIOS
+ changed the computation of variables rhowater_* so that they are real densities (factor 1/rvap missing ; this doesn't affect the previous PEM results as these densities were only compared between each other)
+ added comments and units for ice table variables in physiq_mod.F
+ made Clapeyron coefficient names in physiq_mod.F coherent with how they are defined in the PEM
+ fixed a reference in constants_marspem_mod.F90
+ fixed unit attribute of surface/soil water densities in field_def_physics_mars.xml

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

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

Download in other formats:

Original Format