Context Navigation

physiq_mod.F @ 3468

Last change on this file since 3468 was 3468, checked in by emillour, 4 weeks ago
Mars PCM: Remove obsolete/depreciated lwrite flag (which would trigger some very specific extra text outputs), in code and in reference callphys.def files. EM
File size: 161.9 KB

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

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

Download in other formats:

Original Format