Context Navigation

physiq_mod.F @ 2642

Last change on this file since 2642 was 2628, checked in by abierjon, 4 years ago

Mars GCM:
Big changes on mountain top dust flows for GCM6:

the scheme now activates only in grid meshes that contain a mountain among a hard-written list, instead of every meshes. This is done to prevent strong artificial reinjections of dust in places that don't present huge converging slopes enabling the concentration of dust (ex: Valles Marineris, Hellas). Topdust is now also detrained as soon as it leaves the column it originated from.
the list of the 19 allowed mountains is used by the subroutine topmons_setup in module topmons_mod, to compute a logical array contains_mons(ngrid). alpha_hmons and hsummit are also set up once and for all by this subroutine, which is called in physiq_mod's firstcall.
contains_mons, alpha_hmons and hsummit are now saved variables of the module surfdat_h, and are called as such and not as arguments in the subroutines using them.
the logical flag "slpwind" and the comments in the code have also been updated to the new terminology "mountain top dust flows", accordingly to ticket #71. The new flag read in callphys.def is "topflows".

File size: 160.1 KB

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

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

Download in other formats:

Original Format