Context Navigation

physiq_mod.F @ 3369

Last change on this file since 3369 was 3369, checked in by emillour, 6 months ago
Mars PCM: Change the way the rate of outputs for diagfi.nc files is specified: IMPORTANT: Specifying "ecritphy" no longer possible and will trigger an error. Use "outputs_per_sol" to specify output rate instead. This should makes things (hopefully) clearer for users and also better enforces a cleaner and clearer separation between dynamics and physics. EM
File size: 162.6 KB

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

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

Download in other formats:

Original Format