Context Navigation

physiq_mod.F @ 2429

Last change on this file since 2429 was 2417, checked in by emillour, 5 years ago

Mars GCM:
Add a new scheme to handle the dust and its radiative impact. Triggered using
a new flag dustscaling_mode=2 (dustscaling_mod=0: no rescaling at all, and
dustscaling_mode=1: full rescaling using tauscaling, GCMv5.3 style). Rescaling
is then only done on the radiative impact (see dust_scaling_mod.F90) of dust.
Moreover the scaling factor "dust_rad_adjust" is evaluated using the target dust
scenario opacity for the next sol and left to evolve linearly until then to not
impose the diurnal evolution of dust.
In practice, main changes or additions in the code are:

renamed flag "tauscaling_mode" as "dustscaling_mode"
moved parameter "t_scenario_sol" to "dust_param_mod"
adapted "compute_dustscaling" routine in "dust_scaling_mod"
added module "dust_rad_adjust_mod"
2D fields "dust_rad_adjust_prev" and "dust_rad_adjust_next" required to compute coefficient "dust_rad_adjust" need to be stored in (re)startfi files

File size: 152.0 KB

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

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

Download in other formats:

Original Format