Context Navigation

physiq_mod.F @ 2240

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

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

Download in other formats:

Original Format