Context Navigation

physiq_mod.F @ 3466

Last change on this file since 3466 was 3466, checked in by emillour, 4 weeks ago

Mars PCM:
More tidying in aeronomars:

remove unused "inv.F" and remove "dtridgl.F" which is not used here and is a duplicate of the "dtridgl" routine in phymars/swr_toon.F
turn aeronomars routines to modules, for those which aren't in modules yet.

File size: 162.7 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,lwrite,
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	c *******************************************************************
2549
2550	c ---------------------
2551	c Outputs to the screen
2552	c ---------------------
2553
2554	IF (lwrite) THEN
2555	write(,)'Global diagnostics for the physics'
2556	write(,)'Variables and their increments x and dx/dt * dt'
2557	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
2558	WRITE(*,'(2f10.5,2f15.5)')
2559	s tsurf(igout,:),zdtsurf(igout,:)*ptimestep,
2560	s zplev(igout,1),pdpsrf(igout)*ptimestep
2561	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
2562	WRITE(*,'(i4,6f10.5)') (l,
2563	s pu(igout,l),pdu(igout,l)*ptimestep,
2564	s pv(igout,l),pdv(igout,l)*ptimestep,
2565	s pt(igout,l),pdt(igout,l)*ptimestep,
2566	s l=1,nlayer)
2567	ENDIF ! of IF (lwrite)
2568
2569	c ----------------------------------------------------------
2570	c ----------------------------------------------------------
2571	c INTERPOLATIONS IN THE SURFACE-LAYER
2572	c ----------------------------------------------------------
2573	c ----------------------------------------------------------
2574
2575	n_out=0 ! number of elements in the z_out array.
2576	! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set
2577	! to 5
2578	IF (n_out .ne. 0) THEN
2579
2580	IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out))
2581	IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out))
2582	IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out))
2583
2584	z_out(:)=[3.,2.,1.,0.5,0.1]
2585	u_out(:,:)=0.
2586	T_out(:,:)=0.
2587
2588	call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
2589	& g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,q2,tsurf(:,iflat),
2590	& zh,zq(:,1,igcm_h2o_vap),qsurf(:,igcm_h2o_ice,iflat),mmean(:,1),
2591	& z_out,n_out,T_out,u_out,ustar,tstar,vhf,vvv)
2592	! pourquoi ustar recalcule ici? fait dans vdifc.
2593
2594	#ifndef MESOSCALE
2595	DO n=1,n_out
2596	write(zstring, '(F8.6)') z_out(n)
2597	call write_output('T_out_'//trim(zstring),
2598	& 'potential temperature at z_out','K',T_out(:,n))
2599	call write_output('u_out_'//trim(zstring),
2600	& 'horizontal velocity norm at z_out','m/s',u_out(:,n))
2601	ENDDO
2602	call write_output('u_star',
2603	& 'friction velocity','m/s',ustar)
2604	call write_output('teta_star',
2605	& 'friction potential temperature','K',tstar)
2606	call write_output('vvv',
2607	& 'Vertical velocity variance at zout','m',vvv)
2608	call write_output('vhf',
2609	& 'Vertical heat flux at zout','m',vhf)
2610	#else
2611	T_out1(:)=T_out(:,1)
2612	u_out1(:)=u_out(:,1)
2613	#endif
2614
2615	ENDIF
2616
2617	c ----------------------------------------------------------
2618	c ----------------------------------------------------------
2619	c END OF SURFACE LAYER INTERPOLATIONS
2620	c ----------------------------------------------------------
2621	c ----------------------------------------------------------
2622
2623	#ifndef MESOSCALE
2624	c -------------------------------------------------------------------
2625	c Writing NetCDF file "RESTARTFI" at the end of the run
2626	c -------------------------------------------------------------------
2627	c Note: 'restartfi' is stored just before dynamics are stored
2628	c in 'restart'. Between now and the writting of 'restart',
2629	c there will have been the itau=itau+1 instruction and
2630	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
2631	c thus we store for time=time+dtvr
2632
2633	! default: not writing a restart file at this time step
2634	write_restart=.false.
2635	IF (ecritstart.GT.0) THEN
2636	! For when we store multiple time steps in the restart file
2637	IF (MODULO(icount*iphysiq,ecritstart).EQ.0) THEN
2638	write_restart=.true.
2639	ENDIF
2640	ENDIF
2641	IF (lastcall) THEN
2642	! Always write a restart at the end of the simulation
2643	write_restart=.true.
2644	ENDIF
2645
2646	IF (write_restart) THEN
2647	IF (grid_type==unstructured) THEN !IF DYNAMICO
2648
2649	! When running Dynamico, no need to add a dynamics time step to ztime_fin
2650	IF (ptime.LE. 1.E-10) THEN
2651	! Residual ptime occurs with Dynamico
2652	ztime_fin = pday !+ ptime + ptimestep/(float(iphysiq)*daysec)
2653	. - day_ini - time_phys
2654	ELSE
2655	ztime_fin = pday + ptime !+ ptimestep/(float(iphysiq)*daysec)
2656	. - day_ini - time_phys
2657	ENDIF
2658	if (ecritstart==0) then
2659	ztime_fin = ztime_fin-(day_end-day_ini)
2660	endif
2661
2662	ELSE ! IF LMDZ
2663
2664	if (ecritstart.GT.0) then !IF MULTIPLE RESTARTS nothing change
2665	ztime_fin = pday - day_ini + ptime
2666	& + ptimestep/(float(iphysiq)*daysec)
2667	else !IF ONE RESTART final time in top of day_end
2668	ztime_fin = pday - day_ini-(day_end-day_ini)
2669	& + ptime + ptimestep/(float(iphysiq)*daysec)
2670	endif
2671
2672	ENDIF ! of IF (grid_type==unstructured)
2673	write(*,'(A,I7,A,F12.5)')
2674	. 'PHYSIQ: writing in restartfi ; icount=',
2675	. icount,' date=',ztime_fin
2676
2677	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,
2678	. ptimestep,ztime_fin,
2679	. tsurf,tsoil,inertiesoil,albedo,
2680	. emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,
2681	. wstar,watercap,perennial_co2ice)
2682	ENDIF ! of IF (write_restart)
2683
2684	#endif
2685
2686	c IF (ngrid.NE.1) then
2687
2688	c -------------------------------------------------------------------
2689	c Calculation of diagnostic variables written in both stats and
2690	c diagfi files
2691	c -------------------------------------------------------------------
2692	do ig=1,ngrid
2693	if(mu0(ig).le.0.01) then
2694	fluxsurf_dir_dn_sw(ig) = 0.
2695	else
2696	if (water) then
2697	! both water and dust contribute
2698	fluxsurf_dir_dn_sw(ig) = flux_1AU/dist_sol/dist_solmu0(ig)
2699	& exp(-(tau(ig,iaer_dust_doubleq)+
2700	& tau(ig,iaer_h2o_ice))/mu0(ig))
2701	else
2702	! only dust contributes
2703	fluxsurf_dir_dn_sw(ig) = flux_1AU/dist_sol/dist_solmu0(ig)
2704	& exp(-(tau(ig,iaer_dust_doubleq))/mu0(ig))
2705	endif ! of if (water)
2706	endif ! of if(mu0(ig).le.0.01)
2707	enddo
2708
2709	! Density-scaled opacities
2710	do ig=1,ngrid
2711	dsodust(ig,:) =
2712	& dsodust(ig,:)*tauscaling(ig)
2713	dsords(ig,:) =
2714	& dsords(ig,:)*tauscaling(ig)
2715	dsotop(ig,:) =
2716	& dsotop(ig,:)*tauscaling(ig)
2717	enddo
2718
2719	if(doubleq) then
2720	do ig=1,ngrid
2721	IF (sedimentation) THEN
2722	dqdustsurf(ig) =
2723	& zdqssed(ig,igcm_dust_mass)*tauscaling(ig)
2724	dndustsurf(ig) =
2725	& zdqssed(ig,igcm_dust_number)*tauscaling(ig)
2726	ENDIF
2727	ndust(ig,:) =
2728	& zq(ig,:,igcm_dust_number)*tauscaling(ig)
2729	qdust(ig,:) =
2730	& zq(ig,:,igcm_dust_mass)*tauscaling(ig)
2731	enddo
2732	if (scavenging) then
2733	do ig=1,ngrid
2734	IF (sedimentation) THEN
2735	dqdustsurf(ig) = dqdustsurf(ig) +
2736	& zdqssed(ig,igcm_ccn_mass)*tauscaling(ig)
2737	dndustsurf(ig) = dndustsurf(ig) +
2738	& zdqssed(ig,igcm_ccn_number)*tauscaling(ig)
2739	ENDIF
2740	nccn(ig,:) =
2741	& zq(ig,:,igcm_ccn_number)*tauscaling(ig)
2742	qccn(ig,:) =
2743	& zq(ig,:,igcm_ccn_mass)*tauscaling(ig)
2744	enddo
2745	endif
2746	endif ! of (doubleq)
2747
2748	if (rdstorm) then ! diagnostics of stormdust tendancies for 1D and 3D
2749	mstormdtot(:)=0
2750	mdusttot(:)=0
2751	qdusttotal(:,:)=0
2752	do ig=1,ngrid
2753	rdsdqdustsurf(ig) =
2754	& zdqssed(ig,igcm_stormdust_mass)*tauscaling(ig)
2755	rdsdndustsurf(ig) =
2756	& zdqssed(ig,igcm_stormdust_number)*tauscaling(ig)
2757	rdsndust(ig,:) =
2758	& pq(ig,:,igcm_stormdust_number)*tauscaling(ig)
2759	rdsqdust(ig,:) =
2760	& pq(ig,:,igcm_stormdust_mass)*tauscaling(ig)
2761	do l=1,nlayer
2762	mstormdtot(ig) = mstormdtot(ig) +
2763	& zq(ig,l,igcm_stormdust_mass) *
2764	& (zplev(ig,l) - zplev(ig,l+1)) / g
2765	mdusttot(ig) = mdusttot(ig) +
2766	& zq(ig,l,igcm_dust_mass) *
2767	& (zplev(ig,l) - zplev(ig,l+1)) / g
2768	qdusttotal(ig,l) = qdust(ig,l)+rdsqdust(ig,l) !calculate total dust
2769	enddo
2770	enddo
2771	endif !(rdstorm)
2772
2773	if (water) then
2774	mtot(:)=0
2775	icetot(:)=0
2776	rave(:)=0
2777	tauTES(:)=0
2778
2779	IF (hdo) then
2780	mtotD(:)=0
2781	icetotD(:)=0
2782	ENDIF !hdo
2783
2784	do ig=1,ngrid
2785	do l=1,nlayer
2786	mtot(ig) = mtot(ig) +
2787	& zq(ig,l,igcm_h2o_vap) *
2788	& (zplev(ig,l) - zplev(ig,l+1)) / g
2789	icetot(ig) = icetot(ig) +
2790	& zq(ig,l,igcm_h2o_ice) *
2791	& (zplev(ig,l) - zplev(ig,l+1)) / g
2792	IF (hdo) then
2793	mtotD(ig) = mtotD(ig) +
2794	& zq(ig,l,igcm_hdo_vap) *
2795	& (zplev(ig,l) - zplev(ig,l+1)) / g
2796	icetotD(ig) = icetotD(ig) +
2797	& zq(ig,l,igcm_hdo_ice) *
2798	& (zplev(ig,l) - zplev(ig,l+1)) / g
2799	ENDIF !hdo
2800
2801	c Computing abs optical depth at 825 cm-1 in each
2802	c layer to simulate NEW TES retrieval
2803	Qabsice = min(
2804	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
2805	& )
2806	opTES(ig,l)= 0.75 * Qabsice *
2807	& zq(ig,l,igcm_h2o_ice) *
2808	& (zplev(ig,l) - zplev(ig,l+1)) / g
2809	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
2810	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
2811	enddo
2812	c rave(ig)=rave(ig)/max(icetot(ig),1.e-30) ! mass weight
2813	c if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
2814	enddo
2815	call watersat(ngrid*nlayer,zt,zplay,zqsat)
2816	satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
2817
2818	if (scavenging) then
2819	Nccntot(:)= 0
2820	Mccntot(:)= 0
2821	rave(:)=0
2822	do ig=1,ngrid
2823	do l=1,nlayer
2824	Nccntot(ig) = Nccntot(ig) +
2825	& zq(ig,l,igcm_ccn_number)*tauscaling(ig)
2826	& *(zplev(ig,l) - zplev(ig,l+1)) / g
2827	Mccntot(ig) = Mccntot(ig) +
2828	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
2829	& *(zplev(ig,l) - zplev(ig,l+1)) / g
2830	cccc Column integrated effective ice radius
2831	cccc is weighted by total ice surface area (BETTER than total ice mass)
2832	rave(ig) = rave(ig) +
2833	& tauscaling(ig) *
2834	& zq(ig,l,igcm_ccn_number) *
2835	& (zplev(ig,l) - zplev(ig,l+1)) / g *
2836	& rice(ig,l) * rice(ig,l)* (1.+nuice_ref)
2837	enddo
2838	rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75
2839	& /max(pi*rave(ig),1.e-30) ! surface weight
2840	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
2841	enddo
2842	else ! of if (scavenging)
2843	rave(:)=0
2844	do ig=1,ngrid
2845	do l=1,nlayer
2846	rave(ig) = rave(ig) +
2847	& zq(ig,l,igcm_h2o_ice) *
2848	& (zplev(ig,l) - zplev(ig,l+1)) / g *
2849	& rice(ig,l) * (1.+nuice_ref)
2850	enddo
2851	rave(ig) = max(rave(ig) /
2852	& max(icetot(ig),1.e-30),1.e-30) ! mass weight
2853	enddo
2854	endif ! of if (scavenging)
2855
2856	!Alternative A. Pottier weighting
2857	rave2(:) = 0.
2858	totrave2(:) = 0.
2859	do ig=1,ngrid
2860	do l=1,nlayer
2861	rave2(ig) =rave2(ig)+ zq(ig,l,igcm_h2o_ice)*rice(ig,l)
2862	totrave2(ig) = totrave2(ig) + zq(ig,l,igcm_h2o_ice)
2863	end do
2864	rave2(ig)=max(rave2(ig)/max(totrave2(ig),1.e-30),1.e-30)
2865	end do
2866
2867	endif ! of if (water)
2868
2869	if (co2clouds) then
2870	mtotco2(1:ngrid) = 0.
2871	icetotco2(1:ngrid) = 0.
2872	vaptotco2(1:ngrid) = 0.
2873	do ig=1,ngrid
2874	do l=1,nlayer
2875	vaptotco2(ig) = vaptotco2(ig) +
2876	& zq(ig,l,igcm_co2) *
2877	& (zplev(ig,l) - zplev(ig,l+1)) / g
2878	icetotco2(ig) = icetot(ig) +
2879	& zq(ig,l,igcm_co2_ice) *
2880	& (zplev(ig,l) - zplev(ig,l+1)) / g
2881	end do
2882	mtotco2(ig) = icetotco2(ig) + vaptotco2(ig)
2883	end do
2884	end if
2885
2886	#ifndef MESOSCALE
2887	c -----------------------------------------------------------------
2888	c WSTATS: Saving statistics
2889	c -----------------------------------------------------------------
2890	c ("stats" stores and accumulates key variables in file "stats.nc"
2891	c which can later be used to make the statistic files of the run:
2892	c if flag "callstats" from callphys.def is .true.)
2893
2894	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
2895	call wstats(ngrid,"tsurf","Surface temperature","K",2
2896	& ,tsurf(:,iflat))
2897	call wstats(ngrid,"co2ice","CO2 ice cover",
2898	& "kg.m-2",2,qsurf(:,igcm_co2,iflat))
2899	call wstats(ngrid,"watercap","H2O ice cover",
2900	& "kg.m-2",2,watercap(:,iflat))
2901	call wstats(ngrid,"tau_pref_scenario",
2902	& "prescribed visible dod at 610 Pa","NU",
2903	& 2,tau_pref_scenario)
2904	call wstats(ngrid,"tau_pref_gcm",
2905	& "visible dod at 610 Pa in the GCM","NU",
2906	& 2,tau_pref_gcm)
2907	call wstats(ngrid,"fluxsurf_lw",
2908	& "Thermal IR radiative flux to surface","W.m-2",2,
2909	& fluxsurf_lw(:,iflat))
2910	call wstats(ngrid,"fluxsurf_dn_sw",
2911	& "Incoming Solar radiative flux to surface","W.m-2",2,
2912	& fluxsurf_dn_sw_tot(:,iflat))
2913	call wstats(ngrid,"fluxsurf_up_sw",
2914	& "Reflected Solar radiative flux from surface","W.m-2",2,
2915	& fluxsurf_up_sw_tot)
2916	call wstats(ngrid,"fluxtop_lw",
2917	& "Thermal IR radiative flux to space","W.m-2",2,
2918	& fluxtop_lw)
2919	call wstats(ngrid,"fluxtop_dn_sw",
2920	& "Incoming Solar radiative flux from space","W.m-2",2,
2921	& fluxtop_dn_sw_tot)
2922	call wstats(ngrid,"fluxtop_up_sw",
2923	& "Outgoing Solar radiative flux to space","W.m-2",2,
2924	& fluxtop_up_sw_tot)
2925	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
2926	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
2927	call wstats(ngrid,"v","Meridional (North-South) wind",
2928	& "m.s-1",3,zv)
2929	call wstats(ngrid,"w","Vertical (down-up) wind",
2930	& "m.s-1",3,pw)
2931	call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
2932	call wstats(ngrid,"pressure","Pressure","Pa",3,zplay)
2933	call wstats(ngrid,"q2",
2934	& "Boundary layer eddy kinetic energy",
2935	& "m2.s-2",3,q2)
2936	call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
2937	& emis(:,iflat))
2938	call wstats(ngrid,"fluxsurf_dir_dn_sw",
2939	& "Direct incoming SW flux at surface",
2940	& "W.m-2",2,fluxsurf_dir_dn_sw)
2941
2942	if (calltherm) then
2943	call wstats(ngrid,"zmax_th","Height of thermals",
2944	& "m",2,zmax_th)
2945	call wstats(ngrid,"hfmax_th","Max thermals heat flux",
2946	& "K.m/s",2,hfmax_th)
2947	call wstats(ngrid,"wstar",
2948	& "Max vertical velocity in thermals",
2949	& "m/s",2,wstar)
2950	endif
2951
2952	if (water) then
2953	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
2954	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2955	call wstats(ngrid,"vmr_h2ovap",
2956	& "H2O vapor volume mixing ratio","mol/mol",
2957	& 3,vmr)
2958	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
2959	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
2960	call wstats(ngrid,"vmr_h2oice",
2961	& "H2O ice volume mixing ratio","mol/mol",
2962	& 3,vmr)
2963	! also store vmr_ice*rice for better diagnostics of rice
2964	vmr(1:ngrid,1:nlayer)=vmr(1:ngrid,1:nlayer)*
2965	& rice(1:ngrid,1:nlayer)
2966	call wstats(ngrid,"vmr_h2oice_rice",
2967	& "H2O ice mixing ratio times ice particule size",
2968	& "(mol/mol)*m",
2969	& 3,vmr)
2970	vmr=zqsat(1:ngrid,1:nlayer)
2971	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2972	call wstats(ngrid,"vmr_h2osat",
2973	& "saturation volume mixing ratio","mol/mol",
2974	& 3,vmr)
2975	call wstats(ngrid,"h2o_ice_s",
2976	& "surface h2o_ice","kg/m2",
2977	& 2,qsurf(1,igcm_h2o_ice,iflat))
2978	call wstats(ngrid,'albedo',
2979	& 'albedo',
2980	& '',2,albedo(1,1,iflat))
2981	call wstats(ngrid,"mtot",
2982	& "total mass of water vapor","kg/m2",
2983	& 2,mtot)
2984	call wstats(ngrid,"icetot",
2985	& "total mass of water ice","kg/m2",
2986	& 2,icetot)
2987	call wstats(ngrid,"reffice",
2988	& "Mean reff","m",
2989	& 2,rave)
2990	call wstats(ngrid,"Nccntot",
2991	& "condensation nuclei","Nbr/m2",
2992	& 2,Nccntot)
2993	call wstats(ngrid,"Mccntot",
2994	& "condensation nuclei mass","kg/m2",
2995	& 2,Mccntot)
2996	call wstats(ngrid,"rice",
2997	& "Ice particle size","m",
2998	& 3,rice)
2999	if (.not.activice) then
3000	call wstats(ngrid,"tauTESap",
3001	& "tau abs 825 cm-1","",
3002	& 2,tauTES)
3003	else
3004	call wstats(ngrid,'tauTES',
3005	& 'tau abs 825 cm-1',
3006	& '',2,taucloudtes)
3007	endif
3008
3009	endif ! of if (water)
3010
3011	if (co2clouds) then
3012	call wstats(ngrid,"mtotco2",
3013	& "total mass atm of co2","kg/m2",
3014	& 2,mtotco2)
3015	call wstats(ngrid,"icetotco2",
3016	& "total mass atm of co2 ice","kg/m2",
3017	& 2,icetotco2)
3018	call wstats(ngrid,"vaptotco2",
3019	& "total mass atm of co2 vapor","kg/m2",
3020	& 2,vaptotco2)
3021	end if
3022
3023
3024	if (dustbin.ne.0) then
3025
3026	call wstats(ngrid,'tau','taudust','SI',2,tau(1,1))
3027
3028	if (doubleq) then
3029	call wstats(ngrid,'dqsdust',
3030	& 'deposited surface dust mass',
3031	& 'kg.m-2.s-1',2,dqdustsurf)
3032	call wstats(ngrid,'dqndust',
3033	& 'deposited surface dust number',
3034	& 'number.m-2.s-1',2,dndustsurf)
3035	call wstats(ngrid,'reffdust','reffdust',
3036	& 'm',3,rdust*ref_r0)
3037	call wstats(ngrid,'dustq','Dust mass mr',
3038	& 'kg/kg',3,qdust)
3039	call wstats(ngrid,'dustN','Dust number',
3040	& 'part/kg',3,ndust)
3041	if (rdstorm) then
3042	call wstats(ngrid,'reffstormdust','reffdust',
3043	& 'm',3,rstormdust*ref_r0)
3044	call wstats(ngrid,'rdsdustq','Dust mass mr',
3045	& 'kg/kg',3,rdsqdust)
3046	call wstats(ngrid,'rdsdustN','Dust number',
3047	& 'part/kg',3,rdsndust)
3048	end if
3049	else
3050	do iq=1,dustbin
3051	write(str2(1:2),'(i2.2)') iq
3052	call wstats(ngrid,'q'//str2,'mix. ratio',
3053	& 'kg/kg',3,zq(1,1,iq))
3054	call wstats(ngrid,'qsurf'//str2,'qsurf',
3055	& 'kg.m-2',2,qsurf(1,iq,iflat))
3056	end do
3057	endif ! (doubleq)
3058
3059	if (scavenging) then
3060	call wstats(ngrid,'ccnq','CCN mass mr',
3061	& 'kg/kg',3,qccn)
3062	call wstats(ngrid,'ccnN','CCN number',
3063	& 'part/kg',3,nccn)
3064	endif ! (scavenging)
3065
3066	endif ! (dustbin.ne.0)
3067
3068	if (photochem) then
3069	do iq=1,nq
3070	if (noms(iq) .ne. "dust_mass" .and.
3071	$ noms(iq) .ne. "dust_number" .and.
3072	$ noms(iq) .ne. "ccn_mass" .and.
3073	$ noms(iq) .ne. "ccn_number" .and.
3074	$ noms(iq) .ne. "ccnco2_mass" .and.
3075	$ noms(iq) .ne. "ccnco2_number" .and.
3076	$ noms(iq) .ne. "stormdust_mass" .and.
3077	$ noms(iq) .ne. "stormdust_number" .and.
3078	$ noms(iq) .ne. "topdust_mass" .and.
3079	$ noms(iq) .ne. "topdust_number") then
3080	! volume mixing ratio
3081
3082	vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
3083	& *mmean(1:ngrid,1:nlayer)/mmol(iq)
3084
3085	call wstats(ngrid,"vmr_"//trim(noms(iq)),
3086	$ "Volume mixing ratio","mol/mol",3,vmr)
3087	if ((noms(iq).eq."o")
3088	$ .or. (noms(iq).eq."co2")
3089	$ .or. (noms(iq).eq."o3")
3090	$ .or. (noms(iq).eq."ar")
3091	$ .or. (noms(iq).eq."o2")
3092	$ .or. (noms(iq).eq."h2o_vap") ) then
3093	call write_output("vmr_"//trim(noms(iq)),
3094	$ "Volume mixing ratio","mol/mol",vmr(:,:))
3095	end if
3096
3097	! number density (molecule.cm-3)
3098
3099	rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
3100	& rho(1:ngrid,1:nlayer)n_avog/
3101	& (1000*mmol(iq))
3102
3103	call wstats(ngrid,"num_"//trim(noms(iq)),
3104	$ "Number density","cm-3",3,rhopart)
3105	call write_output("num_"//trim(noms(iq)),
3106	$ "Number density","cm-3",rhopart(:,:))
3107
3108	! vertical column (molecule.cm-2)
3109
3110	do ig = 1,ngrid
3111	colden(ig,iq) = 0.
3112	end do
3113	do l=1,nlayer
3114	do ig=1,ngrid
3115	colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
3116	$ *(zplev(ig,l)-zplev(ig,l+1))
3117	$ 6.022e22/(mmol(iq)g)
3118	end do
3119	end do
3120
3121	call wstats(ngrid,"c_"//trim(noms(iq)),
3122	$ "column","mol cm-2",2,colden(1,iq))
3123	call write_output("c_"//trim(noms(iq)),
3124	$ "column","mol cm-2",colden(:,iq))
3125
3126	! global mass (g)
3127
3128	call planetwide_sumval(colden(:,iq)/6.022e23
3129	$ mmol(iq)1.e4*cell_area(:),mass(iq))
3130
3131	call write_output("mass_"//trim(noms(iq)),
3132	$ "global mass","g",mass(iq))
3133
3134	end if ! of if (noms(iq) .ne. "dust_mass" ...)
3135	end do ! of do iq=1,nq
3136	end if ! of if (photochem)
3137
3138	IF(lastcall.and.callstats) THEN
3139	write (,) "Writing stats..."
3140	call mkstats(ierr)
3141	ENDIF
3142
3143	c (Store EOF for Mars Climate database software)
3144	IF (calleofdump) THEN
3145	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
3146	ENDIF
3147	#endif
3148	!endif of ifndef MESOSCALE
3149
3150	#ifdef MESOSCALE
3151
3152	!! see comm_wrf.
3153	!! not needed when an array is already in a shared module.
3154	!! --> example : hfmax_th, zmax_th
3155
3156	CALL allocate_comm_wrf(ngrid,nlayer)
3157
3158	!state real HR_SW ikj misc 1 - h "HR_SW" "HEATING RATE SW" "K/s"
3159	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
3160	!state real HR_LW ikj misc 1 - h "HR_LW" "HEATING RATE LW" "K/s"
3161	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
3162	!state real SWDOWNZ ij misc 1 - h "SWDOWNZ" "DOWNWARD SW FLUX AT SURFACE" "W m-2"
3163	comm_SWDOWNZ(1:ngrid) = fluxsurf_dn_sw_tot(1:ngrid)
3164	!state real TAU_DUST ij misc 1 - h "TAU_DUST" "REFERENCE VISIBLE DUST OPACITY" ""
3165	comm_TAU_DUST(1:ngrid) = tau_pref_gcm(1:ngrid)
3166	!state real RDUST ikj misc 1 - h "RDUST" "DUST RADIUS" "m"
3167	comm_RDUST(1:ngrid,1:nlayer) = rdust(1:ngrid,1:nlayer)
3168	!state real QSURFDUST ij misc 1 - h "QSURFDUST" "DUST MASS AT SURFACE" "kg m-2"
3169	IF (igcm_dust_mass .ne. 0) THEN
3170	comm_QSURFDUST(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
3171	ELSE
3172	comm_QSURFDUST(1:ngrid) = 0.
3173	ENDIF
3174	!state real MTOT ij misc 1 - h "MTOT" "TOTAL MASS WATER VAPOR in pmic" "pmic"
3175	comm_MTOT(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
3176	!state real ICETOT ij misc 1 - h "ICETOT" "TOTAL MASS WATER ICE" "kg m-2"
3177	comm_ICETOT(1:ngrid) = icetot(1:ngrid) * 1.e6 / rho_ice
3178	!state real VMR_ICE ikj misc 1 - h "VMR_ICE" "VOL. MIXING RATIO ICE" "ppm"
3179	IF (igcm_h2o_ice .ne. 0) THEN
3180	comm_VMR_ICE(1:ngrid,1:nlayer) = 1.e6
3181	. * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
3182	. * mmean(1:ngrid,1:nlayer) / mmol(igcm_h2o_ice)
3183	ELSE
3184	comm_VMR_ICE(1:ngrid,1:nlayer) = 0.
3185	ENDIF
3186	!state real TAU_ICE ij misc 1 - h "TAU_ICE" "CLOUD OD at 825 cm-1 TES" ""
3187	if (activice) then
3188	comm_TAU_ICE(1:ngrid) = taucloudtes(1:ngrid)
3189	else
3190	comm_TAU_ICE(1:ngrid) = tauTES(1:ngrid)
3191	endif
3192	!state real RICE ikj misc 1 - h "RICE" "ICE RADIUS" "m"
3193	comm_RICE(1:ngrid,1:nlayer) = rice(1:ngrid,1:nlayer)
3194
3195	!! calculate sensible heat flux in W/m2 for outputs
3196	!! -- the one computed in vdifc is not the real one
3197	!! -- vdifc must have been called
3198	if (.not.callrichsl) then
3199	sensibFlux(1:ngrid) = zflubid(1:ngrid)
3200	. - capcal(1:ngrid)*zdtsdif(1:ngrid)
3201	else
3202	sensibFlux(1:ngrid) =
3203	& (pplay(1:ngrid,1)/(rpt(1:ngrid,1)))cpp
3204	& sqrt(pu(1:ngrid,1)pu(1:ngrid,1)+pv(1:ngrid,1)*pv(1:ngrid,1)
3205	& +(log(1.+0.7wstar(1:ngrid) + 2.3wstar(1:ngrid)2))2)
3206	& zcdh(1:ngrid)(tsurf(1:ngrid)-zh(1:ngrid,1))
3207	endif
3208
3209	#else
3210	#ifndef MESOINI
3211
3212	c ==========================================================
3213	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
3214	c any variable for diagnostic (output with period
3215	c "ecritphy", set in "run.def")
3216	c ==========================================================
3217	c WRITEDIAGFI can ALSO be called from any other subroutines
3218	c for any variables !!
3219	call write_output("emis","Surface emissivity","",
3220	& emis(:,iflat))
3221	do islope=1,nslope
3222	write(str2(1:2),'(i2.2)') islope
3223	call write_output("emis_slope"//str2,
3224	& "Surface emissivity","",emis(:,islope))
3225	ENDDO
3226	call write_output("zzlay","Midlayer altitude",
3227	& "m",zzlay(:,:))
3228	call write_output("zzlev","Interlayer altitude",
3229	& "m",zzlev(:,1:nlayer))
3230	call write_output("pphi","Geopotential","m2s-2",
3231	& pphi(:,:))
3232	call write_output("phisfi","Surface geopotential",
3233	& "m2s-2",phisfi(:))
3234	if (grid_type == regular_lonlat) then
3235	call write_output("area","Mesh area","m2",
3236	& cell_area_for_lonlat_outputs)
3237	else ! unstructured grid (e.g. dynamico)
3238	call write_output("area","Mesh area","m2",cell_area)
3239	endif
3240	call write_output("tsurf","Surface temperature","K",
3241	& tsurf(:,iflat))
3242	do islope=1,nslope
3243	write(str2(1:2),'(i2.2)') islope
3244	call write_output("tsurf_slope"//str2,
3245	& "Surface temperature","K",
3246	& tsurf(:,islope))
3247	ENDDO
3248	call write_output("ps","surface pressure","Pa",ps(:))
3249	call write_output("co2ice","co2 ice thickness"
3250	& ,"kg.m-2",qsurf(:,igcm_co2,iflat))
3251	do islope=1,nslope
3252	write(str2(1:2),'(i2.2)') islope
3253	call write_output("co2ice_slope"//str2,"co2 ice thickness"
3254	& ,"kg.m-2",qsurf(:,igcm_co2,islope))
3255	ENDDO
3256	call write_output("watercap","Perennial water ice thickness"
3257	& ,"kg.m-2",watercap(:,iflat))
3258	do islope=1,nslope
3259	write(str2(1:2),'(i2.2)') islope
3260	call write_output("watercap_slope"//str2,
3261	& "Perennial water ice thickness"
3262	& ,"kg.m-2",watercap(:,islope))
3263	ENDDO
3264	call write_output("perennial_co2ice",
3265	& "Perennial co2 ice thickness","kg.m-2",
3266	& perennial_co2ice(:,iflat))
3267	do islope=1,nslope
3268	write(str2(1:2),'(i2.2)') islope
3269	call write_output("perennial_co2ice_slope"//str2,
3270	& "Perennial co2 ice thickness"
3271	& ,"kg.m-2",perennial_co2ice(:,islope))
3272	ENDDO
3273	call write_output("temp_layer1","temperature in layer 1",
3274	& "K",zt(:,1))
3275	call write_output("temp7","temperature in layer 7",
3276	& "K",zt(:,7))
3277	call write_output("fluxsurf_lw","fluxsurf_lw","W.m-2",
3278	& fluxsurf_lw(:,iflat))
3279	do islope=1,nslope
3280	write(str2(1:2),'(i2.2)') islope
3281	call write_output("fluxsurf_lw_slope"//str2,
3282	& "fluxsurf_lw","W.m-2",
3283	& fluxsurf_lw(:,islope))
3284	ENDDO
3285	call write_output("fluxsurf_dn_sw","fluxsurf_dn_sw",
3286	& "W.m-2",fluxsurf_dn_sw_tot(:,iflat))
3287	do islope=1,nslope
3288	write(str2(1:2),'(i2.2)') islope
3289	call write_output("fluxsurf_dn_sw_slope"//str2,
3290	& "fluxsurf_dn_sw",
3291	& "W.m-2",fluxsurf_dn_sw_tot(:,islope))
3292	ENDDO
3293	call write_output("fluxtop_dn_sw","fluxtop_dn_sw",
3294	& "W.m-2",fluxtop_dn_sw(:,1) + fluxtop_dn_sw(:,2))
3295	call write_output("fluxtop_lw","fluxtop_lw","W.m-2",
3296	& fluxtop_lw(:))
3297	call write_output("fluxtop_up_sw","fluxtop_up_sw","W.m-2",
3298	& fluxtop_up_sw_tot(:))
3299	call write_output("temp","temperature","K",zt(:,:))
3300	call write_output("Sols","Time","sols",zday)
3301	call write_output("Ls","Solar longitude","deg",
3302	& zls*180./pi)
3303	call write_output("u","Zonal wind","m.s-1",zu(:,:))
3304	call write_output("v","Meridional wind","m.s-1",zv(:,:))
3305	call write_output("w","Vertical wind","m.s-1",pw(:,:))
3306	call write_output("rho","density","kg.m-3",rho(:,:))
3307	call write_output("pressure","Pressure","Pa",zplay(:,:))
3308	call write_output("zplev","Interlayer pressure","Pa",
3309	& zplev(:,1:nlayer))
3310	call write_output('sw_htrt','sw heat. rate',
3311	& 'K/s',zdtsw(:,:))
3312	call write_output('lw_htrt','lw heat. rate',
3313	& 'K/s',zdtlw(:,:))
3314	call write_output("local_time","Local time",
3315	& 'sol',local_time(:))
3316	if (water) then
3317	if (.not.activice) then
3318	CALL write_output('tauTESap',
3319	& 'tau abs 825 cm-1',
3320	& '',tauTES(:))
3321	else
3322	CALL write_output('tauTES',
3323	& 'tau abs 825 cm-1',
3324	& '',taucloudtes(:))
3325	endif
3326	endif ! of if (water)
3327	#else
3328	!!! this is to ensure correct initialisation of mesoscale model
3329	call write_output("tsurf","Surface temperature","K",
3330	& tsurf(:,iflat))
3331	call write_output("ps","surface pressure","Pa",ps(:))
3332	call write_output("co2ice","co2 ice thickness","kg.m-2",
3333	& qsurf(:,igcm_co2,iflat))
3334	call write_output("temp","temperature","K",zt(:,:))
3335	call write_output("u","Zonal wind","m.s-1",zu(:,:))
3336	call write_output("v","Meridional wind","m.s-1",zv(:,:))
3337	call write_output("emis","Surface emissivity","",
3338	& emis(:,iflat))
3339	call write_output("tsoil","Soil temperature",
3340	& "K",tsoil(:,:,iflat))
3341	call write_output("inertiedat","Soil inertia",
3342	& "K",inertiedat(:,:))
3343	#endif
3344
3345	c ----------------------------------------------------------
3346	c Outputs of the CO2 cycle
3347	c ----------------------------------------------------------
3348	if (igcm_co2.ne.0) then
3349	call write_output("co2","co2 mass mixing ratio",
3350	& "kg.kg-1",zq(:,:,igcm_co2))
3351
3352	if (co2clouds) then
3353	call write_output('ccnqco2','CCNco2 mmr',
3354	& 'kg.kg-1',zq(:,:,igcm_ccnco2_mass))
3355
3356	call write_output('ccnNco2','CCNco2 number',
3357	& 'part.kg-1',zq(:,:,igcm_ccnco2_number))
3358
3359	call write_output('co2_ice','co2_ice mmr in atm',
3360	& 'kg.kg-1',zq(:,:,igcm_co2_ice))
3361
3362	call write_output("mtotco2","total mass atm of co2",
3363	& "kg.m-2",mtotco2(:))
3364	call write_output("icetotco2","total mass atm of co2 ice",
3365	& "kg.m-2", icetotco2(:))
3366	call write_output("vaptotco2","total mass atm of co2 "//
3367	& "vapor","kg.m-2", vaptotco2(:))
3368	if (co2useh2o) then
3369	call write_output('ccnqco2_h2o_m_ice',
3370	& 'CCNco2_h2o_mass_ice mmr',
3371	& 'kg.kg-1',zq(:,:,igcm_ccnco2_h2o_mass_ice))
3372
3373	call write_output('ccnqco2_h2o_m_ccn',
3374	& 'CCNco2_h2o_mass_ccn mmr',
3375	& 'kg.kg-1',zq(:,:,igcm_ccnco2_h2o_mass_ccn))
3376
3377	call write_output('ccnNco2_h2o','CCNco2_h2o number',
3378	& 'part.kg-1',zq(:,:,igcm_ccnco2_h2o_number))
3379	end if
3380
3381	if (meteo_flux) then
3382	call write_output('ccnqco2_meteor','CCNco2_meteor mmr',
3383	& 'kg.kg-1',zq(:,:,igcm_ccnco2_meteor_mass))
3384
3385	call write_output('ccnNco2_meteor','CCNco2_meteor number',
3386	& 'part.kg-1',zq(:,:,igcm_ccnco2_meteor_number))
3387	end if
3388
3389	end if ! of if (co2clouds)
3390	end if ! of if (igcm_co2.ne.0)
3391
3392	! Output He tracer, if there is one
3393	if (igcm_he.ne.0) then
3394	call write_output("he","helium mass mixing ratio",
3395	& "kg/kg",zq(:,:,igcm_he))
3396	vmr = zq(1:ngrid,1:nlayer,igcm_he)
3397	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_he)
3398	call write_output('vmr_he','helium vol. mixing ratio',
3399	& 'mol/mol',vmr(:,:))
3400	end if
3401
3402	c ----------------------------------------------------------
3403	c Outputs of the water cycle
3404	c ----------------------------------------------------------
3405	if (water) then
3406	#ifdef MESOINI
3407	!!!! waterice = q01, voir readmeteo.F90
3408	call write_output('q01',noms(igcm_h2o_ice),
3409	& 'kg/kg',
3410	& zq(:,:,igcm_h2o_ice))
3411	!!!! watervapor = q02, voir readmeteo.F90
3412	call write_output('q02',noms(igcm_h2o_vap),
3413	& 'kg/kg',
3414	& zq(:,:,igcm_h2o_vap))
3415	!!!! surface waterice qsurf02 (voir readmeteo)
3416	call write_output('qsurf02','surface tracer',
3417	& 'kg.m-2',
3418	& qsurf(:,igcm_h2o_ice,iflat))
3419	do islope=1,nslope
3420	write(str2(1:2),'(i2.2)') islope
3421	call write_output('qsurf02_slope'//str2,
3422	& 'surface tracer','kg.m-2',
3423	& qsurf(:,igcm_h2o_ice,islope))
3424	ENDDO
3425	#endif
3426	call write_output('mtot',
3427	& 'total mass of water vapor',
3428	& 'kg/m2',mtot(:))
3429	call write_output('icetot',
3430	& 'total mass of water ice',
3431	& 'kg/m2',icetot(:))
3432	vmr = zq(1:ngrid,1:nlayer,igcm_h2o_ice)
3433	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
3434	call write_output('vmr_h2oice','h2o ice vmr',
3435	& 'mol/mol',vmr(:,:))
3436	vmr = zq(1:ngrid,1:nlayer,igcm_h2o_vap)
3437	& * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
3438	call write_output('vmr_h2ovap','h2o vap vmr',
3439	& 'mol/mol',vmr(:,:))
3440	call write_output('reffice',
3441	& 'Mean reff',
3442	& 'm',rave(:))
3443	call write_output('h2o_ice','h2o_ice','kg/kg',
3444	& zq(:,:,igcm_h2o_ice))
3445	call write_output('h2o_vap','h2o_vap','kg/kg',
3446	& zq(:,:,igcm_h2o_vap))
3447
3448	if (hdo) then
3449	vmr=zq(1:ngrid,1:nlayer,igcm_hdo_ice)
3450	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_ice)
3451	CALL write_output('vmr_hdoice','hdo ice vmr',
3452	& 'mol/mol',vmr(:,:))
3453	vmr=zq(1:ngrid,1:nlayer,igcm_hdo_vap)
3454	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_vap)
3455	CALL write_output('vmr_hdovap','hdo vap vmr',
3456	& 'mol/mol',vmr(:,:))
3457	call write_output('hdo_ice','hdo_ice','kg/kg',
3458	& zq(:,:,igcm_hdo_ice))
3459	call write_output('hdo_vap','hdo_vap','kg/kg',
3460	& zq(:,:,igcm_hdo_vap))
3461
3462	CALL write_output('mtotD',
3463	& 'total mass of HDO vapor',
3464	& 'kg/m2',mtotD(:))
3465	CALL write_output('icetotD',
3466	& 'total mass of HDO ice',
3467	& 'kg/m2',icetotD(:))
3468
3469	C Calculation of the D/H ratio
3470	do l=1,nlayer
3471	do ig=1,ngrid
3472	if (zq(ig,l,igcm_h2o_vap).gt.qperemin) then
3473	DoH_vap(ig,l) = ( zq(ig,l,igcm_hdo_vap)/
3474	& zq(ig,l,igcm_h2o_vap) )1./(2.155.76e-6)
3475	else
3476	DoH_vap(ig,l) = 0.
3477	endif
3478	enddo
3479	enddo
3480
3481	do l=1,nlayer
3482	do ig=1,ngrid
3483	if (zq(ig,l,igcm_h2o_ice).gt.qperemin) then
3484	DoH_ice(ig,l) = ( zq(ig,l,igcm_hdo_ice)/
3485	& zq(ig,l,igcm_h2o_ice) )/(2.*155.76e-6)
3486	else
3487	DoH_ice(ig,l) = 0.
3488	endif
3489	enddo
3490	enddo
3491
3492	CALL write_output('DoH_vap',
3493	& 'D/H ratio in vapor',
3494	& ' ',DoH_vap(:,:))
3495	CALL write_output('DoH_ice',
3496	& 'D/H ratio in ice',
3497	& '',DoH_ice(:,:))
3498
3499	endif !hdo
3500
3501	!A. Pottier
3502	! CALL write_output('rmoym',
3503	! & 'alternative reffice',
3504	! & 'm',rave2(:))
3505	call write_output('h2o_saturation',
3506	& 'h2o vap saturation ratio','',satu(:,:))
3507	if (scavenging) then
3508	CALL write_output("Nccntot",
3509	& "condensation nuclei","Nbr/m2",
3510	& Nccntot(:))
3511	CALL write_output("Mccntot",
3512	& "mass condensation nuclei","kg/m2",
3513	& Mccntot(:))
3514	endif
3515	call write_output('rice','Water ice particle size',
3516	& 'm',rice(:,:))
3517	call write_output('h2o_ice_s',
3518	& 'surface h2o_ice',
3519	& 'kg.m-2',qsurf(:,igcm_h2o_ice,iflat))
3520	do islope=1,nslope
3521	write(str2(1:2),'(i2.2)') islope
3522	call write_output('h2o_ice_s_slope'//str2,
3523	& 'surface h2o_ice',
3524	& 'kg.m-2',qsurf(:,igcm_h2o_ice,islope))
3525	ENDDO
3526	if (hdo) then
3527	call write_output('hdo_ice_s',
3528	& 'surface hdo_ice',
3529	& 'kg.m-2',qsurf(:,igcm_hdo_ice,iflat))
3530	do islope=1,nslope
3531	write(str2(1:2),'(i2.2)') islope
3532	call write_output('hdo_ice_s_slope'//str2,
3533	& 'surface hdo_ice',
3534	& 'kg.m-2',qsurf(:,igcm_hdo_ice,islope))
3535	ENDDO
3536
3537	do ig=1,ngrid
3538	if (qsurf_meshavg(ig,igcm_h2o_ice).gt.qperemin) then
3539	DoH_surf(ig) = 0.5*( qsurf_meshavg(ig,igcm_hdo_ice)/
3540	& qsurf_meshavg(ig,igcm_h2o_ice) )/155.76e-6
3541	else
3542	DoH_surf(ig) = 0.
3543	endif
3544	enddo
3545
3546	call write_output('DoH_surf',
3547	& 'surface D/H',
3548	& '',DoH_surf(:))
3549	endif ! hdo
3550
3551	CALL write_output('albedo',
3552	& 'albedo',
3553	& '',albedo(:,1,iflat))
3554	do islope=1,nslope
3555	write(str2(1:2),'(i2.2)') islope
3556	CALL write_output('albedo_slope'//str2,
3557	& 'albedo',
3558	& '',albedo(:,1,islope))
3559	ENDDO
3560	if (surfaceice_tifeedback.or.poreice_tifeedback) then
3561	call write_output("soiltemp",
3562	& "Soil temperature","K",
3563	& tsoil(:,:,iflat))
3564	call write_output('soilti',
3565	& 'Soil Thermal Inertia',
3566	& 'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,iflat))
3567
3568	do islope=1,nslope
3569	write(str2(1:2),'(i2.2)') islope
3570	call write_output('soilti_slope'//str2,
3571	& 'Soil Thermal Inertia',
3572	& 'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,islope))
3573	ENDDO
3574	endif
3575	!A. Pottier
3576	if (CLFvarying) then !AP14 nebulosity
3577	call write_output('totcloudfrac',
3578	& 'Total cloud fraction',
3579	& ' ',totcloudfrac(:))
3580	end if !clf varying
3581	end if !(water)
3582
3583	c ----------------------------------------------------------
3584	c Outputs of the dust cycle
3585	c ----------------------------------------------------------
3586
3587	call write_output('tau_pref_scenario',
3588	& 'Prescribed visible dust optical depth at 610Pa',
3589	& 'NU',tau_pref_scenario(:))
3590
3591	call write_output('tau_pref_gcm',
3592	& 'Visible dust optical depth at 610Pa in the GCM',
3593	& 'NU',tau_pref_gcm(:))
3594
3595	if (reff_driven_IRtoVIS_scenario) then
3596	call write_output('IRtoVIScoef',
3597	& 'Conversion coeff for dust tau from abs9.3um to ext0.67um',
3598	& '/',IRtoVIScoef(:))
3599	endif
3600
3601	if (dustbin.ne.0) then
3602
3603	#ifndef MESOINI
3604	if (doubleq) then
3605	call write_output('dqsdust',
3606	& 'deposited surface dust mass',
3607	& 'kg.m-2.s-1',dqdustsurf(:))
3608	call write_output('dqndust',
3609	& 'deposited surface dust number',
3610	& 'number.m-2.s-1',dndustsurf(:))
3611	call write_output('reffdust','reffdust',
3612	& 'm',rdust(:,:)*ref_r0)
3613	call write_output('dustq','Dust mass mr',
3614	& 'kg/kg',qdust(:,:))
3615	call write_output('dustN','Dust number',
3616	& 'part/kg',ndust(:,:))
3617
3618	select case (trim(dustiropacity))
3619	case ("tes")
3620	call write_output('dsodust_TES',
3621	& 'density scaled extinction opacity of std dust at 9.3um(TES)',
3622	& 'm2.kg-1',dsodust(:,:))
3623	call write_output('dso_TES',
3624	& 'density scaled extinction opacity of all dust at 9.3um(TES)',
3625	& 'm2.kg-1',dsodust(:,:)+dsords(:,:)+dsotop(:,:))
3626	case ("mcs")
3627	call write_output('dsodust',
3628	& 'density scaled extinction opacity of std dust at 21.6um(MCS)',
3629	& 'm2.kg-1',dsodust(:,:))
3630	call write_output('dso',
3631	& 'density scaled extinction opacity of all dust at 21.6um(MCS)',
3632	& 'm2.kg-1',dsodust(:,:)+dsords(:,:)+dsotop(:,:))
3633	end select
3634	else ! (doubleq=.false.)
3635	do iq=1,dustbin
3636	write(str2(1:2),'(i2.2)') iq
3637	call write_output('q'//str2,'mix. ratio',
3638	& 'kg/kg',zq(:,:,iq))
3639	call write_output('qsurf'//str2,'qsurf',
3640	& 'kg.m-2',qsurf(:,iq,iflat))
3641	do islope=1,nslope
3642	write(str2(1:2),'(i2.2)') islope
3643	call write_output('qsurf_slope'//str2,'qsurf',
3644	& 'kg.m-2',qsurf(:,iq,islope))
3645	ENDDO
3646	end do
3647	endif ! (doubleq)
3648
3649	if (rdstorm) then ! writediagfi tendencies stormdust tracers
3650	call write_output('reffstormdust','reffstormdust',
3651	& 'm',rstormdust(:,:)*ref_r0)
3652	call write_output('mstormdtot',
3653	& 'total mass of stormdust only',
3654	& 'kg.m-2',mstormdtot(:))
3655	call write_output('mdusttot',
3656	& 'total mass of dust only',
3657	& 'kg.m-2',mdusttot(:))
3658	call write_output('rdsdqsdust',
3659	& 'deposited surface stormdust mass',
3660	& 'kg.m-2.s-1',rdsdqdustsurf(:))
3661	call write_output('rdsdustq','storm Dust mass mr',
3662	& 'kg/kg',rdsqdust(:,:))
3663	call write_output('rdsdustqmodel','storm Dust massmr',
3664	& 'kg/kg',pq(:,:,igcm_stormdust_mass))
3665	call write_output('rdsdustN','storm Dust number',
3666	& 'part/kg',rdsndust(:,:))
3667	call write_output("stormfract",
3668	& "fraction of the mesh, with stormdust","none",
3669	& totstormfract(:))
3670	! call write_output('qsurf',
3671	! & 'stormdust injection',
3672	! & 'kg.m-2',qsurf(:,igcm_stormdust_mass,iflat))
3673	! do islope=1,nslope
3674	! write(str2(1:2),'(i2.2)') islope
3675	! call write_output('qsurf_slope'//str2,
3676	! & 'stormdust injection',
3677	! & 'kg.m-2',qsurf(:,igcm_stormdust_mass,islope))
3678	! ENDDO
3679	! call write_output('pdqsurf',
3680	! & 'tendancy stormdust mass at surface',
3681	! & 'kg.m-2',dqsurf(:,igcm_stormdust_mass,iflat))
3682	! do islope=1,nslope
3683	! write(str2(1:2),'(i2.2)') islope
3684	! call write_output('pdqsurf_slope'//str2,
3685	! & 'tendancy stormdust mass at surface',
3686	! & 'kg.m-2',dqsurf(:,igcm_stormdust_mass,islope))
3687	! ENDDO
3688	call write_output('wspeed_stormdust',
3689	& 'vertical velocity of stormdust',
3690	& 'm/s',wspeed(:,:))
3691	call write_output('zdqsed_dust_mass'
3692	& ,'sedimentation tendency of background dust mmr'
3693	& ,'kg/kg.s-1',
3694	& zdqsed(:,:,igcm_dust_mass))
3695	call write_output('zdqssed_dust_mass'
3696	& ,'sedimentation tendency of background dust on surface'
3697	& ,'kg.m-2.s-1',
3698	& zdqssed(:,igcm_dust_mass))
3699	call write_output('zdqsed_stormdust_mass'
3700	& ,'sedimentation tendency of stormdust dust mmr'
3701	& ,'kg/kg.s-1',
3702	& zdqsed(:,:,igcm_stormdust_mass))
3703	call write_output('zdqsed_dust_number'
3704	& ,'sedimentation tendency of background dust number'
3705	& ,'nbr/kg.s-1',
3706	& zdqsed(:,:,igcm_dust_number))
3707	call write_output('rdust','rdust',
3708	& 'm',rdust(:,:))
3709	call write_output('rstormdust','rstormdust',
3710	& 'm',rstormdust(:,:))
3711
3712	select case (trim(dustiropacity))
3713	case ("tes")
3714	call write_output('dsords_TES',
3715	& 'density scaled extinction opacity of stormdust at 9.3um(TES)',
3716	& 'm2.kg-1',dsords(:,:))
3717	case ("mcs")
3718	call write_output('dsords',
3719	& 'density scaled extinction opacity of stormdust at 21.6um(MCS)',
3720	& 'm2.kg-1',dsords(:,:))
3721	end select
3722	endif ! (rdstorm)
3723
3724	if (topflows) then
3725	call write_output('refftopdust',
3726	& 'Topdust dust effective radius',
3727	& 'm',rtopdust(:,:)*ref_r0)
3728	call write_output('topdustq','top Dust mass mr',
3729	& 'kg/kg',pq(:,:,igcm_topdust_mass))
3730	call write_output('topdustN','top Dust number',
3731	& 'part/kg',pq(:,:,igcm_topdust_number))
3732	select case (trim(dustiropacity))
3733	case ("tes")
3734	call write_output('dsotop_TES',
3735	& 'density scaled extinction opacity of topdust at 9.3um(TES)',
3736	& 'm2.kg-1',dsotop(:,:))
3737	case ("mcs")
3738	call write_output('dsotop',
3739	& 'density scaled extinction opacity of topdust at 21.6um(MCS)',
3740	& 'm2.kg-1',dsotop(:,:))
3741	end select
3742	endif ! (topflows)
3743
3744	if (dustscaling_mode==2) then
3745	call write_output("dust_rad_adjust",
3746	& "radiative adjustment coefficient for dust",
3747	& "",dust_rad_adjust(:))
3748	endif
3749
3750	! if (scavenging) then ! these outputs should be in the scavenging routine
3751	! call write_output('ccnq','CCN mass mr',
3752	! & 'kg/kg',qccn(:,:))
3753	! call write_output('ccnN','CCN number',
3754	! & 'part/kg',nccn(:,:))
3755	! call write_output('surfccnq','Surf nuclei mass mr',
3756	! & 'kg.m-2',qsurf(:,igcm_ccn_mass,iflat))
3757	! do islope=1,nslope
3758	! write(str2(1:2),'(i2.2)') islope
3759	! call write_output('surfccnq_slope'//str2,
3760	! & 'Surf nuclei mass mr',
3761	! & 'kg.m-2',qsurf(:,igcm_ccn_mass,islope))
3762	! ENDDO
3763	! call write_output('surfccnN','Surf nuclei number',
3764	! & 'kg.m-2',qsurf(:,igcm_ccn_number,iflat))
3765	! do islope=1,nslope
3766	! write(str2(1:2),'(i2.2)') islope
3767	! call write_output('surfccnN_slope'//str2,
3768	! & 'Surf nuclei number',
3769	! & 'kg.m-2',qsurf(:,igcm_ccn_number,islope))
3770	! ENDDO
3771	! endif ! (scavenging)
3772
3773	#else
3774	! !!! to initialize mesoscale we need scaled variables
3775	! !!! because this must correspond to starting point for tracers
3776	! call write_output('dustq','Dust mass mr',
3777	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass))
3778	! call write_output('dustN','Dust number',
3779	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number))
3780	! call write_output('ccn','Nuclei mass mr',
3781	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass))
3782	! call write_output('ccnN','Nuclei number',
3783	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number))
3784	if (freedust) then
3785	call write_output('dustq','Dust mass mr',
3786	& 'kg/kg',qdust)
3787	call write_output('dustN','Dust number',
3788	& 'part/kg',ndust)
3789	call write_output('ccn','CCN mass mr',
3790	& 'kg/kg',qccn)
3791	call write_output('ccnN','CCN number',
3792	& 'part/kg',nccn)
3793	else
3794	call write_output('dustq','Dust mass mr',
3795	& 'kg/kg',pq(:,:,igcm_dust_mass))
3796	call write_output('dustN','Dust number',
3797	& 'part/kg',pq(:,:,igcm_dust_number))
3798	call write_output('ccn','Nuclei mass mr',
3799	& 'kg/kg',pq(:,:,igcm_ccn_mass))
3800	call write_output('ccnN','Nuclei number',
3801	& 'part/kg',pq(:,:,igcm_ccn_number))
3802	endif
3803	#endif
3804
3805	end if ! (dustbin.ne.0)
3806
3807	c ----------------------------------------------------------
3808	c Thermospheric outputs
3809	c ----------------------------------------------------------
3810	if(callthermos) then
3811
3812	call write_output("quv","UV heating","K/s",
3813	$ zdteuv(:,:))
3814	call write_output("cond","Thermal conduction","K/s",
3815	$ zdtconduc(:,:))
3816
3817	!H, H2 and D escape fluxes
3818
3819	call write_output("PhiH","H escape flux","s-1",
3820	$ PhiEscH)
3821	call write_output("PhiH2","H2 escape flux","s-1",
3822	$ PhiEscH2)
3823	call write_output("PhiD","D escape flux","s-1",
3824	$ PhiEscD)
3825
3826	endif !(callthermos)
3827
3828	call write_output("q15um","15 um cooling","K/s",
3829	$ zdtnlte(:,:))
3830	call write_output("qnir","NIR heating","K/s",
3831	$ zdtnirco2(:,:))
3832
3833	c ----------------------------------------------------------
3834	c ----------------------------------------------------------
3835	c PBL OUTPUS
3836	c ----------------------------------------------------------
3837	c ----------------------------------------------------------
3838
3839	c ----------------------------------------------------------
3840	c Outputs of thermals
3841	c ----------------------------------------------------------
3842	if (calltherm) then
3843	call write_output('lmax_th',
3844	& 'index of vertical extension of thermals',
3845	& 'grid level',lmax_th_out(:))
3846	call write_output('zmax_th',
3847	& 'vertical extension of thermals','m',
3848	& zmax_th(:))
3849	call write_output('hfmax_th',
3850	& 'maximum heat flux in thermals','K.m/s',
3851	& hfmax_th(:))
3852	call write_output('wstar',
3853	& 'maximum thermals vertical velocity','m/s',
3854	& wstar(:))
3855	end if
3856
3857	c ----------------------------------------------------------
3858	c ----------------------------------------------------------
3859	c END OF PBL OUTPUS
3860	c ----------------------------------------------------------
3861	c ----------------------------------------------------------
3862
3863	c ----------------------------------------------------------
3864	c Output in netcdf file "diagsoil.nc" for subterranean
3865	c variables (output every "ecritphy", as for writediagfi)
3866	c ----------------------------------------------------------
3867	! Write soil temperature
3868	call write_output("soiltemp","Soil temperature","K",
3869	& tsoil(:,:,iflat))
3870	do islope=1,nslope
3871	write(str2(1:2),'(i2.2)') islope
3872	call write_output("soiltemp_slope"//str2,
3873	& "Soil temperature","K",
3874	& tsoil(:,:,islope))
3875	ENDDO
3876
3877	!PREVIOUSLY IN 1D ONLY
3878	call write_output("dtrad","rad. heat. rate",
3879	& "K.s-1",dtrad(:,:))
3880
3881	if (rdstorm) then
3882	call write_output('aerosol_dust','opacity of env. dust',''
3883	& ,aerosol(:,:,iaer_dust_doubleq))
3884	call write_output('aerosol_stormdust',
3885	& 'opacity of storm dust',''
3886	& ,aerosol(:,:,iaer_stormdust_doubleq))
3887	call write_output('dqsdifdustq',
3888	&'tendency due to vertical diffusion of background dust on surface'
3889	& ,'kg.m-2.s-1',zdqsdif(:,igcm_dust_mass,iflat))
3890	do islope=1,nslope
3891	write(str2(1:2),'(i2.2)') islope
3892	call write_output('dqsdifdustq_slope'//str2,
3893	&'tendency due to vertical diffusion of background dust on surface'
3894	& ,'kg.m-2.s-1',zdqsdif(:,igcm_dust_mass,islope))
3895	ENDDO
3896	call write_output('dqsdifrdsq',
3897	& 'tendency due to vertical diffusion of stormdust on surface',
3898	& 'kg.m-2.s-1',zdqsdif(:,igcm_stormdust_mass,iflat))
3899	do islope=1,nslope
3900	write(str2(1:2),'(i2.2)') islope
3901	call write_output('dqsdifrdsq_slope'//str2,
3902	& 'tendency due to vertical diffusion of stormdust on surface',
3903	& 'kg.m-2.s-1',zdqsdif(:,igcm_stormdust_mass,islope))
3904	ENDDO
3905	endif !(rdstorm)
3906
3907	if(water) then
3908	if (.not.scavenging) then
3909	call write_output('zdqcloud_ice','cloud ice',
3910	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_ice))
3911	call write_output('zdqcloud_vap','cloud vap',
3912	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_vap))
3913	call write_output('zdqcloud','cloud',
3914	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_ice)
3915	& +zdqcloud(:,:,igcm_h2o_vap))
3916	IF (hdo) THEN
3917	call write_output('zdqcloud_iceD','cloud ice hdo',
3918	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_hdo_ice))
3919	call write_output('zdqcloud_vapD','cloud vap hdo',
3920	& 'kg.m-2.s-1',zdqcloud(:,:,igcm_hdo_vap))
3921	ENDIF ! hdo
3922	endif !not.scavenging
3923
3924	! Output needed by the PEM
3925	DO ig = 1,ngrid
3926	ztmp1 =(1/m_co2 - 1/m_noco2)
3927	ztmp2=1/m_noco2
3928	pvap_surf(ig) = 1/(ztmp1*zq(ig,1,igcm_co2)+ztmp2)
3929	& * zq(ig,1,igcm_h2o_vap)/(mmol(igcm_h2o_vap)1.e-3)ps(ig)
3930
3931	DO islope = 1,nslope
3932	! Clapeyron law for psat (psat = exp(beta/Th2o+alpha)),following Murphy and Koop 2005
3933	rhowater_surf_sat(ig,islope) =
3934	& exp(beta_clap_h2o/tsurf(ig,islope)+alpha_clap_h2o)
3935	& / tsurf(ig,islope)
3936	& * mmol(igcm_h2o_vap)/(mugaz*r)
3937
3938	if(qsurf(ig,igcm_h2o_ice,islope).gt.(1.e-4)) then
3939	! we consider to be at saturation above 1.e-4 kg.m-2
3940	rhowater_surf(ig,islope) = rhowater_surf_sat(ig,islope)
3941	else
3942	! otherwise, use vapor partial pressure
3943	rhowater_surf(ig,islope) = pvap_surf(ig)
3944	& / tsurf(ig,islope)
3945	& * mmol(igcm_h2o_vap)/(mugaz*r)
3946	endif
3947	DO isoil = 1,nsoilmx
3948	rhowater_soil(ig,isoil,islope) =
3949	& exp(beta_clap_h2o/tsoil(ig,isoil,islope)+alpha_clap_h2o)
3950	& / tsoil(ig,isoil,islope)
3951	& * mmol(igcm_h2o_vap)/(mugaz*r)
3952	ENDDO
3953	ENDDO
3954	ENDDO
3955
3956	CALL write_output("waterdensity_soil",
3957	& "rhowater_soil",'kg.m-3',
3958	& rhowater_soil(:,:,iflat))
3959	CALL write_output("waterdensity_surface",
3960	& "rhowater_surface",'kg.m-3',
3961	& rhowater_surf(:,iflat))
3962	DO islope = 1,nslope
3963	write(str2(1:2),'(i2.2)') islope
3964	CALL write_output("waterdensity_soil_slope"//str2,
3965	& "rhowater_soil_slope"//str2,'kg.m-3',
3966	& rhowater_soil(:,:,islope))
3967	CALL write_output("waterdensity_surface"//str2,
3968	& "rhowater_surface"//str2,'kg.m-3',
3969	& rhowater_surf(:,islope))
3970	ENDDO
3971
3972	CALL write_output("h2o_layer1","h2o mass mr in the first layer",
3973	& 'kg/kg',zq(:,1,igcm_h2o_vap))
3974	CALL write_output("co2_layer1","co2 mass mr in the first layer",
3975	& 'kg/kg',zq(:,1,igcm_co2))
3976	ENDIF ! of IF (water)
3977
3978	!PREVIOUSLY IN 1D ONLY
3979
3980	c ==========================================================
3981	c END OF WRITEDIAGFI
3982	c ==========================================================
3983	#endif
3984	! of ifdef MESOSCALE
3985
3986	c ELSE ! if(ngrid.eq.1)
3987
3988	c#ifndef MESOSCALE
3989	c write(*,
3990	c & '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)')
3991	c & zls*180./pi,odpref,tau_pref_scenario
3992	c#endif
3993
3994	c END IF ! if(ngrid.ne.1)
3995
3996	! test for co2 conservation with co2 microphysics
3997	if (igcm_co2_ice.ne.0) then
3998	co2totB = 0. ! added by C.M.
3999	do ig=1,ngrid
4000	do l=1,nlayer
4001	co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
4002	& (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice)
4003	& +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
4004	enddo
4005	co2totB = co2totB + qsurf(ig,igcm_co2,iflat)
4006	enddo
4007	else
4008	co2totB = 0. ! added by C.M.
4009	do ig=1,ngrid
4010	do l=1,nlayer
4011	co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
4012	& (pq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2)*ptimestep)
4013	enddo
4014	co2totB = co2totB + qsurf(ig,igcm_co2,iflat)
4015	enddo
4016	endif ! of if (igcm_co2_ice.ne.0)
4017	co2conservation = (co2totA-co2totB)/co2totA
4018	call write_output( 'co2conservation',
4019	& 'Total CO2 mass conservation in physic',
4020	& 'kg', co2conservation)
4021	! XIOS outputs
4022	#ifdef CPP_XIOS
4023	! Send fields to XIOS: (NB these fields must also be defined as
4024	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
4025
4026	CALL send_xios_field("controle",tab_cntrl_mod,1)
4027
4028	CALL send_xios_field("ap",ap,1)
4029	CALL send_xios_field("bp",bp,1)
4030	CALL send_xios_field("aps",aps,1)
4031	CALL send_xios_field("bps",bps,1)
4032
4033	if (lastcall.and.is_omp_master) then
4034	write(,) "physiq lastcall: call xios_context_finalize"
4035	call xios_context_finalize
4036	endif
4037	#endif
4038
4039	if (check_physics_outputs) then
4040	! Check the validity of updated fields at the end of the physics step
4041	call check_physics_fields("end of physiq:",zt,zu,zv,zplev,zq)
4042	endif
4043
4044	icount=icount+1
4045
4046	END SUBROUTINE physiq
4047
4048	END MODULE physiq_mod

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

Download in other formats:

Original Format