Context Navigation

physiq_mod.F @ 2616

Last change on this file since 2616 was 2616, checked in by romain.vande, 3 years ago
LMDZ_MARS RV : Open_MP; Put all the "save" variables as "!$OMP THREADPRIVATE" in phymars. The code can now be tested, see README for more info
File size: 160.9 KB

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

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

Download in other formats:

Original Format