Context Navigation

physiq_mod.F @ 2609

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

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

Download in other formats:

Original Format