Context Navigation

physiq_mod.F @ 2225

Last change on this file since 2225 was 2225, checked in by emillour, 5 years ago
Mars GCM: Further cleanup of nonorographic gravity waves drag parametrization east_gwstress and west_gwstress need be initialized and are better as module variables. 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	c input sub-grid scale cloud
1066	& clearsky,totcloudfrac,
1067	c input sub-grid scale topography
1068	& nohmons,alpha_hmons,
1069	c output
1070	& pdqrds,wspeed,dsodust,dsords,
1071	& tauref)
1072
1073	c update the tendencies of both dust after vertical transport
1074	DO l=1,nlayer
1075	DO ig=1,ngrid
1076	pdq(ig,l,igcm_stormdust_mass)=
1077	& pdq(ig,l,igcm_stormdust_mass)+
1078	& pdqrds(ig,l,igcm_stormdust_mass)
1079	pdq(ig,l,igcm_stormdust_number)=
1080	& pdq(ig,l,igcm_stormdust_number)+
1081	& pdqrds(ig,l,igcm_stormdust_number)
1082
1083	pdq(ig,l,igcm_dust_mass)=
1084	& pdq(ig,l,igcm_dust_mass)+ pdqrds(ig,l,igcm_dust_mass)
1085	pdq(ig,l,igcm_dust_number)=
1086	& pdq(ig,l,igcm_dust_number)+
1087	& pdqrds(ig,l,igcm_dust_number)
1088
1089	ENDDO
1090	ENDDO
1091	do l=1,nlayer
1092	do ig=1,ngrid
1093	qdustrds1(ig,l)=pq(ig,l,igcm_dust_mass)+
1094	& pdq(ig,l,igcm_dust_mass)*ptimestep
1095	qstormrds1(ig,l)=pq(ig,l,igcm_stormdust_mass)+
1096	& pdq(ig,l,igcm_stormdust_mass)*ptimestep
1097	qdusttotal1(ig)=qdusttotal1(ig)+(qdustrds1(ig,l)+
1098	& qstormrds1(ig,l))*(zplev(ig,l)-
1099	& zplev(ig,l+1))/g
1100	enddo
1101	enddo
1102
1103	c for diagnostics
1104	call writediagfi(ngrid,'qdustrds1','qdust after rds',
1105	& 'kg/kg ',3,qdustrds1)
1106	call writediagfi(ngrid,'qstormrds1','qstorm after rds',
1107	& 'kg/kg ',3,qstormrds1)
1108
1109	call writediagfi(ngrid,'qdusttotal0','q sum before rds',
1110	& 'kg/kg ',2,qdusttotal0)
1111	call writediagfi(ngrid,'qdusttotal1','q sum after rds',
1112	& 'kg/kg ',2,qdusttotal1)
1113
1114	ENDIF ! end of if(rdstorm)
1115
1116	c 3.2 Dust entrained from the PBL up to the top of sub-grid scale topography
1117	c -------------------------------------------
1118	IF (slpwind) THEN
1119	if (ngrid.gt.1) then
1120	hsummit(:)=summit(:)-phisfi(:)/g
1121	else
1122	hsummit(:)=14000.
1123	endif
1124	nohmons=.false.
1125	pdqtop(:,:,:)=0.
1126	CALL topmons(ngrid,nlayer,nq,ptime,ptimestep,
1127	& pq,pdq,pt,pdt,zplev,zplay,zzlev,
1128	& zzlay,zdtsw,zdtlw,
1129	& icount,zday,zls,tsurf,igout,aerosol,
1130	& totstormfract,clearatm,dsords,
1131	& clearsky,totcloudfrac,
1132	& nohmons,hsummit,
1133	& pdqtop,wtop,dsodust,
1134	& tauref)
1135
1136
1137	c update the tendencies of both dust after vertical transport
1138	DO l=1,nlayer
1139	DO ig=1,ngrid
1140	pdq(ig,l,igcm_topdust_mass)=
1141	& pdq(ig,l,igcm_topdust_mass)+
1142	& pdqtop(ig,l,igcm_topdust_mass)
1143	pdq(ig,l,igcm_topdust_number)=
1144	& pdq(ig,l,igcm_topdust_number)+
1145	& pdqtop(ig,l,igcm_topdust_number)
1146	pdq(ig,l,igcm_dust_mass)=
1147	& pdq(ig,l,igcm_dust_mass)+ pdqtop(ig,l,igcm_dust_mass)
1148	pdq(ig,l,igcm_dust_number)=
1149	& pdq(ig,l,igcm_dust_number)+pdqtop(ig,l,igcm_dust_number)
1150
1151	ENDDO
1152	ENDDO
1153
1154	ENDIF ! end of if (slpwind)
1155
1156	c 3.3 Dust injection from the surface
1157	c -------------------------------------------
1158	if (dustinjection.gt.0) then
1159	CALL compute_dtau(ngrid,nlayer,
1160	& zday,pplev,tauref,
1161	& ptimestep,dustliftday,local_time)
1162	endif ! end of if (dustinjection.gt.0)
1163
1164	c-----------------------------------------------------------------------
1165	c 4. Gravity wave and subgrid scale topography drag :
1166	c -------------------------------------------------
1167
1168
1169	IF(calllott)THEN
1170
1171	CALL calldrag_noro(ngrid,nlayer,ptimestep,
1172	& zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
1173
1174	DO l=1,nlayer
1175	DO ig=1,ngrid
1176	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
1177	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
1178	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
1179	ENDDO
1180	ENDDO
1181	ENDIF
1182
1183	c-----------------------------------------------------------------------
1184	c 5. Vertical diffusion (turbulent mixing):
1185	c -----------------------------------------
1186
1187	IF (calldifv) THEN
1188
1189	DO ig=1,ngrid
1190	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
1191	ENDDO
1192
1193	zdum1(:,:)=0
1194	zdum2(:,:)=0
1195	do l=1,nlayer
1196	do ig=1,ngrid
1197	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
1198	enddo
1199	enddo
1200
1201	c ----------------------
1202	c Treatment of a special case : using new surface layer (Richardson based)
1203	c without using the thermals in gcm and mesoscale can yield problems in
1204	c weakly unstable situations when winds are near to 0. For those cases, we add
1205	c a unit subgrid gustiness. Remember that thermals should be used we using the
1206	c Richardson based surface layer model.
1207	IF ( .not.calltherm
1208	. .and. callrichsl
1209	. .and. .not.turb_resolved) THEN
1210	DO ig=1, ngrid
1211	IF (zh(ig,1) .lt. tsurf(ig)) THEN
1212	wstar(ig)=1.
1213	hfmax_th(ig)=0.2
1214	ELSE
1215	wstar(ig)=0.
1216	hfmax_th(ig)=0.
1217	ENDIF
1218	ENDDO
1219	ENDIF
1220	c ----------------------
1221
1222	IF (tke_heat_flux .ne. 0.) THEN
1223	zz1(:)=(pt(:,1)+pdt(:,1)ptimestep)(r/g)*
1224	& (-alog(zplay(:,1)/zplev(:,1)))
1225	pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
1226	ENDIF
1227
1228	c Calling vdif (Martian version WITH CO2 condensation)
1229	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
1230	$ ptimestep,capcal,lwrite,
1231	$ zplay,zplev,zzlay,zzlev,z0,
1232	$ pu,pv,zh,pq,tsurf,emis,qsurf,
1233	$ zdum1,zdum2,zdh,pdq,zflubid,
1234	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
1235	& zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,
1236	& zcondicea_co2microp,sensibFlux,
1237	& dustliftday,local_time)
1238
1239	DO ig=1,ngrid
1240	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
1241	ENDDO
1242
1243	IF (.not.turb_resolved) THEN
1244	DO l=1,nlayer
1245	DO ig=1,ngrid
1246	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
1247	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
1248	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
1249
1250	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
1251	ENDDO
1252	ENDDO
1253
1254	if (tracer) then
1255	DO iq=1, nq
1256	DO l=1,nlayer
1257	DO ig=1,ngrid
1258	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
1259	ENDDO
1260	ENDDO
1261	ENDDO
1262	DO iq=1, nq
1263	DO ig=1,ngrid
1264	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
1265	ENDDO
1266	ENDDO
1267	end if ! of if (tracer)
1268	ELSE
1269	write (,) '******************************************'
1270	write (,) '** LES mode: the difv part is only used to'
1271	write (,) '** - provide HFX and UST to the dynamics'
1272	write (,) '** - update TSURF'
1273	write (,) '******************************************'
1274	!! Specific treatment for lifting in turbulent-resolving mode (AC)
1275	IF (lifting .and. doubleq) THEN
1276	!! lifted dust is injected in the first layer.
1277	!! Sedimentation must be called after turbulent mixing, i.e. on next step, after WRF.
1278	!! => lifted dust is not incremented before the sedimentation step.
1279	zdqdif(1:ngrid,1,1:nq)=0.
1280	zdqdif(1:ngrid,1,igcm_dust_number) =
1281	. -zdqsdif(1:ngrid,igcm_dust_number)
1282	zdqdif(1:ngrid,1,igcm_dust_mass) =
1283	. -zdqsdif(1:ngrid,igcm_dust_mass)
1284	zdqdif(1:ngrid,2:nlayer,1:nq) = 0.
1285	DO iq=1, nq
1286	IF ((iq .ne. igcm_dust_mass)
1287	& .and. (iq .ne. igcm_dust_number)) THEN
1288	zdqsdif(:,iq)=0.
1289	ENDIF
1290	ENDDO
1291	ELSE
1292	zdqdif(1:ngrid,1:nlayer,1:nq) = 0.
1293	zdqsdif(1:ngrid,1:nq) = 0.
1294	ENDIF
1295	ENDIF
1296	ELSE
1297	DO ig=1,ngrid
1298	zdtsurf(ig)=zdtsurf(ig)+
1299	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
1300	ENDDO
1301	IF (turb_resolved) THEN
1302	write(,) 'Turbulent-resolving mode !'
1303	write(,) 'Please set calldifv to T in callphys.def'
1304	call abort_physic("physiq","turbulent-resolving mode",1)
1305	ENDIF
1306	ENDIF ! of IF (calldifv)
1307
1308	c-----------------------------------------------------------------------
1309	c 6. Thermals :
1310	c -----------------------------
1311
1312	if(calltherm .and. .not.turb_resolved) then
1313
1314	call calltherm_interface(ngrid,nlayer,nq,
1315	$ tracer,igcm_co2,
1316	$ zzlev,zzlay,
1317	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
1318	$ zplay,zplev,pphi,zpopsk,
1319	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
1320	$ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
1321
1322	DO l=1,nlayer
1323	DO ig=1,ngrid
1324	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
1325	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
1326	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
1327	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
1328	ENDDO
1329	ENDDO
1330
1331	DO ig=1,ngrid
1332	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
1333	ENDDO
1334
1335	if (tracer) then
1336	DO iq=1,nq
1337	DO l=1,nlayer
1338	DO ig=1,ngrid
1339	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
1340	ENDDO
1341	ENDDO
1342	ENDDO
1343	endif
1344
1345	lmax_th_out(:)=real(lmax_th(:))
1346
1347	else !of if calltherm
1348	lmax_th(:)=0
1349	wstar(:)=0.
1350	hfmax_th(:)=0.
1351	lmax_th_out(:)=0.
1352	end if
1353	c-----------------------------------------------------------------------
1354	c 7. Dry convective adjustment:
1355	c -----------------------------
1356
1357	IF(calladj) THEN
1358
1359	DO l=1,nlayer
1360	DO ig=1,ngrid
1361	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
1362	ENDDO
1363	ENDDO
1364	zduadj(:,:)=0
1365	zdvadj(:,:)=0
1366	zdhadj(:,:)=0
1367
1368	CALL convadj(ngrid,nlayer,nq,ptimestep,
1369	$ zplay,zplev,zpopsk,lmax_th,
1370	$ pu,pv,zh,pq,
1371	$ pdu,pdv,zdh,pdq,
1372	$ zduadj,zdvadj,zdhadj,
1373	$ zdqadj)
1374
1375	DO l=1,nlayer
1376	DO ig=1,ngrid
1377	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
1378	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
1379	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
1380
1381	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
1382	ENDDO
1383	ENDDO
1384
1385	if(tracer) then
1386	DO iq=1, nq
1387	DO l=1,nlayer
1388	DO ig=1,ngrid
1389	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
1390	ENDDO
1391	ENDDO
1392	ENDDO
1393	end if
1394	ENDIF ! of IF(calladj)
1395
1396	c-----------------------------------------------------
1397	c 8. Non orographic Gravity waves :
1398	c -------------------------------------------------
1399
1400	IF (calllott_nonoro) THEN
1401
1402	CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep,zplay,
1403	& zmax_th, ! max altitude reached by thermals (m)
1404	& pt, pu, pv,
1405	& pdt, pdu, pdv,
1406	& zustrhi,zvstrhi,
1407	& d_t_hin, d_u_hin, d_v_hin)
1408
1409	! Update tendencies
1410	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)
1411	& +d_t_hin(1:ngrid,1:nlayer)
1412	! d_t_hin(:,:)= d_t_hin(:,:)/ptimestep ! K/s (for outputs?)
1413	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)
1414	& +d_u_hin(1:ngrid,1:nlayer)
1415	! d_u_hin(:,:)= d_u_hin(:,:)/ptimestep ! (m/s)/s (for outputs?)
1416	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)
1417	& +d_v_hin(1:ngrid,1:nlayer)
1418	! d_v_hin(:,:)= d_v_hin(:,:)/ptimestep ! (m/s)/s (for outputs?)
1419
1420	ENDIF ! of IF (calllott_nonoro)
1421
1422	c-----------------------------------------------------------------------
1423	c 9. Specific parameterizations for tracers
1424	c: -----------------------------------------
1425
1426	if (tracer) then
1427
1428	c 9a. Water and ice
1429	c ---------------
1430
1431	c ---------------------------------------
1432	c Water ice condensation in the atmosphere
1433	c ----------------------------------------
1434	IF (water) THEN
1435
1436	call watercloud(ngrid,nlayer,ptimestep,
1437	& zplev,zplay,pdpsrf,zzlay, pt,pdt,
1438	& pq,pdq,zdqcloud,zdtcloud,
1439	& nq,tau,tauscaling,rdust,rice,nuice,
1440	& rsedcloud,rhocloud,totcloudfrac)
1441
1442	c Temperature variation due to latent heat release
1443	if (activice) then
1444	pdt(1:ngrid,1:nlayer) =
1445	& pdt(1:ngrid,1:nlayer) +
1446	& zdtcloud(1:ngrid,1:nlayer)
1447	endif
1448
1449	! increment water vapour and ice atmospheric tracers tendencies
1450	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) =
1451	& pdq(1:ngrid,1:nlayer,igcm_h2o_vap) +
1452	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_vap)
1453	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
1454	& pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
1455	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice)
1456
1457	! increment dust and ccn masses and numbers
1458	! We need to check that we have Nccn & Ndust > 0
1459	! This is due to single precision rounding problems
1460	if (microphys) then
1461	pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
1462	& pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
1463	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
1464	pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
1465	& pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
1466	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
1467	where (pq(:,:,igcm_ccn_mass) +
1468	& ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
1469	pdq(:,:,igcm_ccn_mass) =
1470	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1471	pdq(:,:,igcm_ccn_number) =
1472	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1473	end where
1474	where (pq(:,:,igcm_ccn_number) +
1475	& ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
1476	pdq(:,:,igcm_ccn_mass) =
1477	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1478	pdq(:,:,igcm_ccn_number) =
1479	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1480	end where
1481	endif
1482
1483	if (scavenging) then
1484	pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
1485	& pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
1486	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_mass)
1487	pdq(1:ngrid,1:nlayer,igcm_dust_number) =
1488	& pdq(1:ngrid,1:nlayer,igcm_dust_number) +
1489	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_number)
1490	where (pq(:,:,igcm_dust_mass) +
1491	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
1492	pdq(:,:,igcm_dust_mass) =
1493	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1494	pdq(:,:,igcm_dust_number) =
1495	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1496	end where
1497	where (pq(:,:,igcm_dust_number) +
1498	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
1499	pdq(:,:,igcm_dust_mass) =
1500	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1501	pdq(:,:,igcm_dust_number) =
1502	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1503	end where
1504	endif ! of if scavenging
1505
1506	END IF ! of IF (water)
1507
1508	c 9a bis. CO2 clouds (CL & JA)
1509	c ---------------------------------------
1510	c CO2 ice cloud condensation in the atmosphere
1511	c ----------------------------------------
1512	c flag needed in callphys.def:
1513	c co2clouds=.true. is mandatory (default is .false.)
1514	c co2useh2o=.true. if you want to allow co2 condensation
1515	c on water ice particles
1516	c meteo_flux=.true. if you want to add a meteoritic
1517	c supply of CCN
1518	c CLFvaryingCO2=.true. if you want to have a sub-grid
1519	c temperature distribution
1520	c spantCO2=integer (i.e. 3) amplitude of the sub-grid T disti
1521	c nuiceco2_sed=0.2 variance of the size distribution for the
1522	c sedimentation
1523	c nuiceco2_ref=0.2 variance of the size distribution for the
1524	c nucleation
1525	c imicroco2=50 micro-timestep is 1/50 of physical timestep
1526
1527	IF (co2clouds ) THEN
1528
1529
1530	call co2cloud(ngrid,nlayer,ptimestep,
1531	& zplev,zplay,pdpsrf,zzlay,pt,pdt,
1532	& pq,pdq,zdqcloudco2,zdtcloudco2,
1533	& nq,tau,tauscaling,rdust,rice,riceco2,nuice,
1534	& rsedcloudco2,rhocloudco2,
1535	& rsedcloud,rhocloud,zzlev,zdqssed_co2,
1536	& pdu,pu,zcondicea_co2microp)
1537
1538
1539	c Temperature variation due to latent heat release
1540	c if (activice) then !Maybe create activice_co2 ?
1541	pdt(1:ngrid,1:nlayer) =
1542	& pdt(1:ngrid,1:nlayer) +
1543	& zdtcloudco2(1:ngrid,1:nlayer)! --> in co2condens
1544	c endif
1545
1546
1547	! increment dust and ccn masses and numbers
1548	! We need to check that we have Nccn & Ndust > 0
1549	! This is due to single precision rounding problems
1550
1551	! increment dust tracers tendancies
1552	pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
1553	& pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
1554	& zdqcloudco2(1:ngrid,1:nlayer,igcm_dust_mass)
1555	pdq(1:ngrid,1:nlayer,igcm_dust_number) =
1556	& pdq(1:ngrid,1:nlayer,igcm_dust_number) +
1557	& zdqcloudco2(1:ngrid,1:nlayer,igcm_dust_number)
1558	pdq(1:ngrid,1:nlayer,igcm_co2) =
1559	& pdq(1:ngrid,1:nlayer,igcm_co2) +
1560	& zdqcloudco2(1:ngrid,1:nlayer,igcm_co2)
1561	pdq(1:ngrid,1:nlayer,igcm_co2_ice) =
1562	& pdq(1:ngrid,1:nlayer,igcm_co2_ice) +
1563	& zdqcloudco2(1:ngrid,1:nlayer,igcm_co2_ice)
1564	pdq(1:ngrid,1:nlayer,igcm_ccnco2_mass) =
1565	& pdq(1:ngrid,1:nlayer,igcm_ccnco2_mass) +
1566	& zdqcloudco2(1:ngrid,1:nlayer,igcm_ccnco2_mass)
1567	pdq(1:ngrid,1:nlayer,igcm_ccnco2_number) =
1568	& pdq(1:ngrid,1:nlayer,igcm_ccnco2_number) +
1569	& zdqcloudco2(1:ngrid,1:nlayer,igcm_ccnco2_number)
1570	!Update water ice clouds values as well
1571	if (co2useh2o) then
1572	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
1573	& pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
1574	& zdqcloudco2(1:ngrid,1:nlayer,igcm_h2o_ice)
1575	pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
1576	& pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
1577	& zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_mass)
1578	pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
1579	& pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
1580	& zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
1581	where (pq(:,:,igcm_ccn_mass) +
1582	& ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
1583	pdq(:,:,igcm_ccn_mass) =
1584	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1585	pdq(:,:,igcm_ccn_number) =
1586	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1587	end where
1588	where (pq(:,:,igcm_ccn_number) +
1589	& ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
1590	pdq(:,:,igcm_ccn_mass) =
1591	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1592	pdq(:,:,igcm_ccn_number) =
1593	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1594	end where
1595	endif ! of if (co2useh2o)
1596	c Negative values?
1597	where (pq(:,:,igcm_ccnco2_mass) +
1598	& ptimestep*pdq(:,:,igcm_ccnco2_mass) < 0.)
1599	pdq(:,:,igcm_ccnco2_mass) =
1600	& - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
1601	pdq(:,:,igcm_ccnco2_number) =
1602	& - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
1603	end where
1604	where (pq(:,:,igcm_ccnco2_number) +
1605	& ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.)
1606	pdq(:,:,igcm_ccnco2_mass) =
1607	& - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
1608	pdq(:,:,igcm_ccnco2_number) =
1609	& - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
1610	end where
1611
1612	c Negative values?
1613	where (pq(:,:,igcm_dust_mass) +
1614	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
1615	pdq(:,:,igcm_dust_mass) =
1616	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1617	pdq(:,:,igcm_dust_number) =
1618	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1619	end where
1620	where (pq(:,:,igcm_dust_number) +
1621	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
1622	pdq(:,:,igcm_dust_mass) =
1623	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1624	pdq(:,:,igcm_dust_number) =
1625	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1626	end where
1627
1628	END IF ! of IF (co2clouds)
1629
1630	c 9b. Aerosol particles
1631	c -------------------
1632	c ----------
1633	c Dust devil :
1634	c ----------
1635	IF(callddevil) then
1636	call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
1637	& zdqdev,zdqsdev)
1638
1639	if (dustbin.ge.1) then
1640	do iq=1,nq
1641	DO l=1,nlayer
1642	DO ig=1,ngrid
1643	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
1644	ENDDO
1645	ENDDO
1646	enddo
1647	do iq=1,nq
1648	DO ig=1,ngrid
1649	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
1650	ENDDO
1651	enddo
1652	endif ! of if (dustbin.ge.1)
1653
1654	END IF ! of IF (callddevil)
1655
1656	c -------------
1657	c Sedimentation : acts also on water ice
1658	c -------------
1659	IF (sedimentation) THEN
1660	zdqsed(1:ngrid,1:nlayer,1:nq)=0
1661	zdqssed(1:ngrid,1:nq)=0
1662
1663	c Sedimentation for co2 clouds tracers are inside co2cloud microtimestep
1664	c Zdqssed isn't
1665	call callsedim(ngrid,nlayer,ptimestep,
1666	& zplev,zzlev,zzlay,pt,pdt,
1667	& rdust,rstormdust,rtopdust,
1668	& rice,rsedcloud,rhocloud,
1669	& pq,pdq,zdqsed,zdqssed,nq,
1670	& tau,tauscaling)
1671	c Flux at the surface of co2 ice computed in co2cloud microtimestep
1672	IF (co2clouds) THEN
1673	zdqssed(1:ngrid,igcm_co2_ice)=zdqssed_co2(1:ngrid)
1674	ENDIF
1675
1676	IF (rdstorm) THEN
1677	c Storm dust cannot sediment to the surface
1678	DO ig=1,ngrid
1679	zdqsed(ig,1,igcm_stormdust_mass)=
1680	& zdqsed(ig,1,igcm_stormdust_mass)+
1681	& zdqssed(ig,igcm_stormdust_mass) /
1682	& ((pplev(ig,1)-pplev(ig,2))/g)
1683	zdqsed(ig,1,igcm_stormdust_number)=
1684	& zdqsed(ig,1,igcm_stormdust_number)+
1685	& zdqssed(ig,igcm_stormdust_number) /
1686	& ((pplev(ig,1)-pplev(ig,2))/g)
1687	zdqssed(ig,igcm_stormdust_mass)=0.
1688	zdqssed(ig,igcm_stormdust_number)=0.
1689	ENDDO
1690	ENDIF !rdstorm
1691
1692	DO iq=1, nq
1693	DO l=1,nlayer
1694	DO ig=1,ngrid
1695	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
1696	ENDDO
1697	ENDDO
1698	ENDDO
1699	DO iq=1, nq
1700	DO ig=1,ngrid
1701	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
1702	ENDDO
1703	ENDDO
1704
1705	END IF ! of IF (sedimentation)
1706
1707	c Add lifted dust to tendancies after sedimentation in the LES (AC)
1708	IF (turb_resolved) THEN
1709	DO iq=1, nq
1710	DO l=1,nlayer
1711	DO ig=1,ngrid
1712	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
1713	ENDDO
1714	ENDDO
1715	ENDDO
1716	DO iq=1, nq
1717	DO ig=1,ngrid
1718	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
1719	ENDDO
1720	ENDDO
1721	ENDIF
1722	c
1723	c 9c. Chemical species
1724	c ------------------
1725
1726	#ifndef MESOSCALE
1727	c --------------
1728	c photochemistry :
1729	c --------------
1730	IF (photochem) then
1731
1732	if (modulo(icount-1,ichemistry).eq.0) then
1733	! compute chemistry every ichemistry physics step
1734
1735	! dust and ice surface area
1736	call surfacearea(ngrid, nlayer, naerkind,
1737	$ ptimestep, zplay, zzlay,
1738	$ pt, pq, pdq, nq,
1739	$ rdust, rice, tau, tauscaling,
1740	$ surfdust, surfice)
1741	! call photochemistry
1742	call calchim(ngrid,nlayer,nq,
1743	& ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
1744	$ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
1745	$ zdqcloud,zdqscloud,tau(:,1),co2ice,
1746	$ pu,pdu,pv,pdv,surfdust,surfice)
1747
1748	endif ! of if (modulo(icount-1,ichemistry).eq.0)
1749
1750	! increment values of tracers:
1751	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1752	! tracers is zero anyways
1753	DO l=1,nlayer
1754	DO ig=1,ngrid
1755	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
1756	ENDDO
1757	ENDDO
1758	ENDDO ! of DO iq=1,nq
1759
1760	! add condensation tendency for H2O2
1761	if (igcm_h2o2.ne.0) then
1762	DO l=1,nlayer
1763	DO ig=1,ngrid
1764	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
1765	& +zdqcloud(ig,l,igcm_h2o2)
1766	ENDDO
1767	ENDDO
1768	endif
1769
1770	! increment surface values of tracers:
1771	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1772	! tracers is zero anyways
1773	DO ig=1,ngrid
1774	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
1775	ENDDO
1776	ENDDO ! of DO iq=1,nq
1777
1778	! add condensation tendency for H2O2
1779	if (igcm_h2o2.ne.0) then
1780	DO ig=1,ngrid
1781	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
1782	& +zdqscloud(ig,igcm_h2o2)
1783	ENDDO
1784	endif
1785
1786	END IF ! of IF (photochem)
1787	#endif
1788
1789	endif ! of if (tracer)
1790
1791	#ifndef MESOSCALE
1792	c-----------------------------------------------------------------------
1793	c 10. THERMOSPHERE CALCULATION
1794	c-----------------------------------------------------------------------
1795
1796	if (callthermos) then
1797	call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol,
1798	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
1799	& pt,pq,pu,pv,pdt,pdq,
1800	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
1801
1802	DO l=1,nlayer
1803	DO ig=1,ngrid
1804	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
1805	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
1806	& +zdteuv(ig,l)
1807	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
1808	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
1809	DO iq=1, nq
1810	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
1811	ENDDO
1812	ENDDO
1813	ENDDO
1814
1815	endif ! of if (callthermos)
1816	#endif
1817	c-----------------------------------------------------------------------
1818	c 11. Carbon dioxide condensation-sublimation:
1819	c (should be the last atmospherical physical process to be computed)
1820	c -------------------------------------------
1821
1822	IF (tituscap) THEN
1823	!!! get the actual co2 seasonal cap from Titus observations
1824	CALL geticecover( ngrid, 180.*zls/pi,
1825	. 180.longitude/pi, 180.latitude/pi, co2ice )
1826	co2ice = co2ice * 10000.
1827	ENDIF
1828
1829
1830	pdpsrf(:) = 0
1831
1832	IF (callcond) THEN
1833	CALL co2condens(ngrid,nlayer,nq,ptimestep,
1834	$ capcal,zplay,zplev,tsurf,pt,
1835	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
1836	$ co2ice,albedo,emis,rdust,
1837	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
1838	$ fluxsurf_sw,zls,
1839	$ zdqssed_co2,zcondicea_co2microp,
1840	& zdtcloudco2,zdqsc)
1841
1842	DO l=1,nlayer
1843	DO ig=1,ngrid
1844	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
1845	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
1846	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
1847	ENDDO
1848	ENDDO
1849	DO ig=1,ngrid
1850	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
1851	ENDDO
1852
1853	IF (tracer) THEN
1854	DO iq=1, nq
1855	DO l=1,nlayer
1856	DO ig=1,ngrid
1857	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
1858	ENDDO
1859	ENDDO
1860	ENDDO
1861
1862	DO iq=1, nq
1863	DO ig=1,ngrid
1864	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqsc(ig,iq)
1865	ENDDO ! (ig)
1866	ENDDO ! (iq)
1867
1868	ENDIF ! of IF (tracer)
1869
1870	#ifndef MESOSCALE
1871	! update surface pressure
1872	DO ig=1,ngrid
1873	ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
1874	ENDDO
1875
1876	! update pressure levels
1877	DO l=1,nlayer
1878	DO ig=1,ngrid
1879	zplay(ig,l) = aps(l) + bps(l)*ps(ig)
1880	zplev(ig,l) = ap(l) + bp(l)*ps(ig)
1881	ENDDO
1882	ENDDO
1883	zplev(:,nlayer+1) = 0.
1884
1885	! update layers altitude
1886	DO l=2,nlayer
1887	DO ig=1,ngrid
1888	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
1889	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
1890	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
1891	ENDDO
1892	ENDDO
1893	#endif
1894
1895	ENDIF ! of IF (callcond)
1896
1897	c-----------------------------------------------------------------------
1898	c Updating tracer budget on surface
1899	c-----------------------------------------------------------------------
1900	IF (tracer) THEN
1901	DO iq=1, nq
1902	DO ig=1,ngrid
1903
1904	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1905
1906	ENDDO ! (ig)
1907	ENDDO ! (iq)
1908	ENDIF
1909
1910	c-----------------------------------------------------------------------
1911	c 12. Surface and sub-surface soil temperature
1912	c-----------------------------------------------------------------------
1913	c
1914	c
1915	c 12.1 Increment Surface temperature:
1916	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1917
1918	DO ig=1,ngrid
1919	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1920	ENDDO
1921
1922	c Prescribe a cold trap at south pole (except at high obliquity !!)
1923	c Temperature at the surface is set there to be the temperature
1924	c corresponding to equilibrium temperature between phases of CO2
1925
1926
1927	IF (tracer.AND.water.AND.(ngrid.NE.1)) THEN
1928	!#ifndef MESOSCALE
1929	! if (caps.and.(obliquit.lt.27.)) then => now done in co2condens
1930	! NB: Updated surface pressure, at grid point 'ngrid', is
1931	! ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
1932	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
1933	! & (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
1934	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095ps(ngrid)))
1935	! endif
1936	!#endif
1937	c -------------------------------------------------------------
1938	c Change of surface albedo in case of ground frost
1939	c everywhere except on the north permanent cap and in regions
1940	c covered by dry ice.
1941	c ALWAYS PLACE these lines after co2condens !!!
1942	c -------------------------------------------------------------
1943	do ig=1,ngrid
1944	if ((co2ice(ig).eq.0).and.
1945	& (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
1946	albedo(ig,1) = albedo_h2o_ice
1947	albedo(ig,2) = albedo_h2o_ice
1948	c write(,) "frost thickness", qsurf(ig,igcm_h2o_ice)
1949	c write(,) "physiq.F frost :"
1950	c & ,latitude(ig)180./pi, longitude(ig)180./pi
1951	endif
1952	enddo ! of do ig=1,ngrid
1953	ENDIF ! of IF (tracer.AND.water.AND.(ngrid.NE.1))
1954
1955	c
1956	c 12.2 Compute soil temperatures and subsurface heat flux:
1957	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1958	IF (callsoil) THEN
1959	c Thermal inertia feedback
1960	IF (tifeedback) THEN
1961	CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
1962	CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
1963	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
1964	ELSE
1965	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
1966	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
1967	ENDIF
1968	ENDIF
1969
1970	c To avoid negative values
1971	IF (rdstorm) THEN
1972	where (pq(:,:,igcm_stormdust_mass) +
1973	& ptimestep*pdq(:,:,igcm_stormdust_mass) < 0.)
1974	pdq(:,:,igcm_stormdust_mass) =
1975	& - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30
1976	pdq(:,:,igcm_stormdust_number) =
1977	& - pq(:,:,igcm_stormdust_number)/ptimestep + 1.e-30
1978	end where
1979	where (pq(:,:,igcm_stormdust_number) +
1980	& ptimestep*pdq(:,:,igcm_stormdust_number) < 0.)
1981	pdq(:,:,igcm_stormdust_mass) =
1982	& - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30
1983	pdq(:,:,igcm_stormdust_number) =
1984	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1985	end where
1986
1987	where (pq(:,:,igcm_dust_mass) +
1988	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
1989	pdq(:,:,igcm_dust_mass) =
1990	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1991	pdq(:,:,igcm_dust_number) =
1992	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1993	end where
1994	where (pq(:,:,igcm_dust_number) +
1995	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
1996	pdq(:,:,igcm_dust_mass) =
1997	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1998	pdq(:,:,igcm_dust_number) =
1999	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
2000	end where
2001	ENDIF !(rdstorm)
2002
2003	c-----------------------------------------------------------------------
2004	c 13. Write output files
2005	c ----------------------
2006
2007	c -------------------------------
2008	c Dynamical fields incrementation
2009	c -------------------------------
2010	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
2011	! temperature, zonal and meridional wind
2012	DO l=1,nlayer
2013	DO ig=1,ngrid
2014	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
2015	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
2016	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
2017	ENDDO
2018	ENDDO
2019
2020	! tracers
2021	DO iq=1, nq
2022	DO l=1,nlayer
2023	DO ig=1,ngrid
2024	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
2025	ENDDO
2026	ENDDO
2027	ENDDO
2028
2029	! Density
2030	DO l=1,nlayer
2031	DO ig=1,ngrid
2032	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
2033	ENDDO
2034	ENDDO
2035
2036	! Potential Temperature
2037
2038	DO ig=1,ngrid
2039	DO l=1,nlayer
2040	zh(ig,l) = zt(ig,l)(zplev(ig,1)/zplay(ig,l))*rcp
2041	ENDDO
2042	ENDDO
2043
2044
2045	c Compute surface stress : (NB: z0 is a common in surfdat.h)
2046	c DO ig=1,ngrid
2047	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
2048	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
2049	c ENDDO
2050
2051	c Sum of fluxes in solar spectral bands (for output only)
2052	DO ig=1,ngrid
2053	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
2054	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
2055	ENDDO
2056	c ***** TEST ****************************************************
2057	ztim1 = 999
2058	DO l=1,nlayer
2059	DO ig=1,ngrid
2060	if (pt(ig,l).lt.ztim1) then
2061	ztim1 = pt(ig,l)
2062	igmin = ig
2063	lmin = l
2064	end if
2065	ENDDO
2066	ENDDO
2067	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
2068	write(,) 'PHYSIQ: stability WARNING :'
2069	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
2070	& 'ig l =', igmin, lmin
2071	end if
2072	c *******************************************************************
2073
2074	c ---------------------
2075	c Outputs to the screen
2076	c ---------------------
2077
2078	IF (lwrite) THEN
2079	PRINT*,'Global diagnostics for the physics'
2080	PRINT,'Variables and their increments x and dx/dt dt'
2081	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
2082	WRITE(*,'(2f10.5,2f15.5)')
2083	s tsurf(igout),zdtsurf(igout)*ptimestep,
2084	s zplev(igout,1),pdpsrf(igout)*ptimestep
2085	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
2086	WRITE(*,'(i4,6f10.5)') (l,
2087	s pu(igout,l),pdu(igout,l)*ptimestep,
2088	s pv(igout,l),pdv(igout,l)*ptimestep,
2089	s pt(igout,l),pdt(igout,l)*ptimestep,
2090	s l=1,nlayer)
2091	ENDIF ! of IF (lwrite)
2092
2093	c ----------------------------------------------------------
2094	c ----------------------------------------------------------
2095	c INTERPOLATIONS IN THE SURFACE-LAYER
2096	c ----------------------------------------------------------
2097	c ----------------------------------------------------------
2098
2099	n_out=0 ! number of elements in the z_out array.
2100	! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set
2101	! to 5
2102	IF (n_out .ne. 0) THEN
2103
2104	IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out))
2105	IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out))
2106	IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out))
2107
2108	z_out(:)=[3.,2.,1.,0.5,0.1]
2109	u_out(:,:)=0.
2110	T_out(:,:)=0.
2111
2112	call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
2113	& g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,n_out,
2114	& T_out,u_out,ustar,tstar,L_mo,vhf,vvv)
2115	! pourquoi ustar recalcule ici? fait dans vdifc.
2116
2117	#ifndef MESOSCALE
2118	IF (ngrid .eq. 1) THEN
2119	dimout=0
2120	ELSE
2121	dimout=2
2122	ENDIF
2123	DO n=1,n_out
2124	write(zstring, '(F8.6)') z_out(n)
2125	call WRITEDIAGFI(ngrid,'T_out_'//trim(zstring),
2126	& 'potential temperature at z_out','K',dimout,T_out(:,n))
2127	call WRITEDIAGFI(ngrid,'u_out_'//trim(zstring),
2128	& 'horizontal velocity norm at z_out','m/s',dimout,u_out(:,n))
2129	ENDDO
2130	call WRITEDIAGFI(ngrid,'u_star',
2131	& 'friction velocity','m/s',dimout,ustar)
2132	call WRITEDIAGFI(ngrid,'teta_star',
2133	& 'friction potential temperature','K',dimout,tstar)
2134	! call WRITEDIAGFI(ngrid,'L',
2135	! & 'Monin Obukhov length','m',dimout,L_mo)
2136	call WRITEDIAGFI(ngrid,'vvv',
2137	& 'Vertical velocity variance at zout','m',dimout,vvv)
2138	call WRITEDIAGFI(ngrid,'vhf',
2139	& 'Vertical heat flux at zout','m',dimout,vhf)
2140	#else
2141	T_out1(:)=T_out(:,1)
2142	u_out1(:)=u_out(:,1)
2143	#endif
2144
2145	ENDIF
2146
2147	c ----------------------------------------------------------
2148	c ----------------------------------------------------------
2149	c END OF SURFACE LAYER INTERPOLATIONS
2150	c ----------------------------------------------------------
2151	c ----------------------------------------------------------
2152
2153	IF (ngrid.NE.1) THEN
2154
2155	#ifndef MESOSCALE
2156	c -------------------------------------------------------------------
2157	c Writing NetCDF file "RESTARTFI" at the end of the run
2158	c -------------------------------------------------------------------
2159	c Note: 'restartfi' is stored just before dynamics are stored
2160	c in 'restart'. Between now and the writting of 'restart',
2161	c there will have been the itau=itau+1 instruction and
2162	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
2163	c thus we store for time=time+dtvr
2164
2165	IF( ((ecritstart.GT.0) .and.
2166	. (MOD(icount*iphysiq,ecritstart).EQ.0))
2167	. .or. lastcall ) THEN
2168
2169	ztime_fin = pday + ptime + ptimestep/(float(iphysiq)*daysec)
2170	. - day_ini - time_phys
2171	print*, pday,ptime,day_ini, time_phys
2172	write(*,'(A,I7,A,F12.5)')
2173	. 'PHYSIQ: Ecriture du fichier restartfi ; icount=',
2174	. icount,' date=',ztime_fin
2175
2176
2177	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
2178	. ptimestep,ztime_fin,
2179	. tsurf,tsoil,co2ice,albedo,emis,
2180	. q2,qsurf,tauscaling,totcloudfrac,wstar,
2181	. mem_Mccn_co2,mem_Nccn_co2,mem_Mh2o_co2)
2182
2183	ENDIF
2184	#endif
2185
2186	c -------------------------------------------------------------------
2187	c Calculation of diagnostic variables written in both stats and
2188	c diagfi files
2189	c -------------------------------------------------------------------
2190
2191	if (tracer) then
2192
2193	if(doubleq) then
2194	do ig=1,ngrid
2195	dqdustsurf(ig) =
2196	& zdqssed(ig,igcm_dust_mass)*tauscaling(ig)
2197	dndustsurf(ig) =
2198	& zdqssed(ig,igcm_dust_number)*tauscaling(ig)
2199	ndust(ig,:) =
2200	& zq(ig,:,igcm_dust_number)*tauscaling(ig)
2201	qdust(ig,:) =
2202	& zq(ig,:,igcm_dust_mass)*tauscaling(ig)
2203	enddo
2204	if (scavenging) then
2205	do ig=1,ngrid
2206	dqdustsurf(ig) = dqdustsurf(ig) +
2207	& zdqssed(ig,igcm_ccn_mass)*tauscaling(ig)
2208	dndustsurf(ig) = dndustsurf(ig) +
2209	& zdqssed(ig,igcm_ccn_number)*tauscaling(ig)
2210	nccn(ig,:) =
2211	& zq(ig,:,igcm_ccn_number)*tauscaling(ig)
2212	qccn(ig,:) =
2213	& zq(ig,:,igcm_ccn_mass)*tauscaling(ig)
2214	enddo
2215	endif
2216	endif ! of (doubleq)
2217
2218	if (rdstorm) then ! diagnostics of stormdust tendancies for 1D and 3D
2219	mstormdtot(:)=0
2220	mdusttot(:)=0
2221	qdusttotal(:,:)=0
2222	do ig=1,ngrid
2223	rdsdqdustsurf(ig) =
2224	& zdqssed(ig,igcm_stormdust_mass)*tauscaling(ig)
2225	rdsdndustsurf(ig) =
2226	& zdqssed(ig,igcm_stormdust_number)*tauscaling(ig)
2227	rdsndust(ig,:) =
2228	& pq(ig,:,igcm_stormdust_number)*tauscaling(ig)
2229	rdsqdust(ig,:) =
2230	& pq(ig,:,igcm_stormdust_mass)*tauscaling(ig)
2231	do l=1,nlayer
2232	mstormdtot(ig) = mstormdtot(ig) +
2233	& zq(ig,l,igcm_stormdust_mass) *
2234	& (zplev(ig,l) - zplev(ig,l+1)) / g
2235	mdusttot(ig) = mdusttot(ig) +
2236	& zq(ig,l,igcm_dust_mass) *
2237	& (zplev(ig,l) - zplev(ig,l+1)) / g
2238	qdusttotal(ig,l) = qdust(ig,l)+rdsqdust(ig,l) !calculate total dust
2239	enddo
2240	enddo
2241	endif !(rdstorm)
2242
2243	if (co2clouds) then
2244	do ig=1,ngrid
2245	nccnco2(ig,:) =
2246	& zq(ig,:,igcm_ccnco2_number)*tauscaling(ig)
2247	qccnco2(ig,:) =
2248	& zq(ig,:,igcm_ccnco2_mass)*tauscaling(ig)
2249	enddo
2250	c D. BARDET compute integrated CO2 vapor and ice content
2251	mtotco2(:)=0
2252	icetotco2(:)=0
2253	do ig=1,ngrid
2254	do l=1,nlayer
2255	mtotco2(ig) = mtotco2(ig) +
2256	& zq(ig,l,igcm_co2) *
2257	& (zplev(ig,l) - zplev(ig,l+1)) / g
2258	icetotco2(ig) = icetotco2(ig) +
2259	& zq(ig,l,igcm_co2_ice) *
2260	& (zplev(ig,l) - zplev(ig,l+1)) / g
2261	enddo
2262	enddo
2263	endif ! of if (co2clouds)
2264
2265	if (water) then
2266	mtot(:)=0
2267	icetot(:)=0
2268	rave(:)=0
2269	tauTES(:)=0
2270	do ig=1,ngrid
2271	do l=1,nlayer
2272	mtot(ig) = mtot(ig) +
2273	& zq(ig,l,igcm_h2o_vap) *
2274	& (zplev(ig,l) - zplev(ig,l+1)) / g
2275	icetot(ig) = icetot(ig) +
2276	& zq(ig,l,igcm_h2o_ice) *
2277	& (zplev(ig,l) - zplev(ig,l+1)) / g
2278	c Computing abs optical depth at 825 cm-1 in each
2279	c layer to simulate NEW TES retrieval
2280	Qabsice = min(
2281	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
2282	& )
2283	opTES(ig,l)= 0.75 * Qabsice *
2284	& zq(ig,l,igcm_h2o_ice) *
2285	& (zplev(ig,l) - zplev(ig,l+1)) / g
2286	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
2287	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
2288	enddo
2289	c rave(ig)=rave(ig)/max(icetot(ig),1.e-30) ! mass weight
2290	c if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
2291	enddo
2292
2293	call watersat(ngrid*nlayer,zt,zplay,zqsat)
2294	satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
2295
2296	if (scavenging) then
2297	Nccntot(:)= 0
2298	Mccntot(:)= 0
2299	rave(:)=0
2300	do ig=1,ngrid
2301	do l=1,nlayer
2302	Nccntot(ig) = Nccntot(ig) +
2303	& zq(ig,l,igcm_ccn_number)*tauscaling(ig)
2304	& *(zplev(ig,l) - zplev(ig,l+1)) / g
2305	Mccntot(ig) = Mccntot(ig) +
2306	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
2307	& *(zplev(ig,l) - zplev(ig,l+1)) / g
2308	cccc Column integrated effective ice radius
2309	cccc is weighted by total ice surface area (BETTER than total ice mass)
2310	rave(ig) = rave(ig) +
2311	& tauscaling(ig) *
2312	& zq(ig,l,igcm_ccn_number) *
2313	& (zplev(ig,l) - zplev(ig,l+1)) / g *
2314	& rice(ig,l) * rice(ig,l)* (1.+nuice_ref)
2315	enddo
2316	rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75
2317	& /max(pi*rave(ig),1.e-30) ! surface weight
2318	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
2319	enddo
2320	else ! of if (scavenging)
2321	rave(:)=0
2322	do ig=1,ngrid
2323	do l=1,nlayer
2324	rave(ig) = rave(ig) +
2325	& zq(ig,l,igcm_h2o_ice) *
2326	& (zplev(ig,l) - zplev(ig,l+1)) / g *
2327	& rice(ig,l) * (1.+nuice_ref)
2328	enddo
2329	rave(ig) = max(rave(ig) /
2330	& max(icetot(ig),1.e-30),1.e-30) ! mass weight
2331	enddo
2332	endif ! of if (scavenging)
2333
2334	!Alternative A. Pottier weighting
2335	rave2(:) = 0.
2336	totrave2(:) = 0.
2337	do ig=1,ngrid
2338	do l=1,nlayer
2339	rave2(ig) =rave2(ig)+ zq(ig,l,igcm_h2o_ice)*rice(ig,l)
2340	totrave2(ig) = totrave2(ig) + zq(ig,l,igcm_h2o_ice)
2341	end do
2342	rave2(ig)=max(rave2(ig)/max(totrave2(ig),1.e-30),1.e-30)
2343	end do
2344
2345	endif ! of if (water)
2346
2347
2348	endif ! of if (tracer)
2349
2350	#ifndef MESOSCALE
2351	c -----------------------------------------------------------------
2352	c WSTATS: Saving statistics
2353	c -----------------------------------------------------------------
2354	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
2355	c which can later be used to make the statistic files of the run:
2356	c "stats") only possible in 3D runs !
2357
2358	IF (callstats) THEN
2359
2360	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
2361	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
2362	call wstats(ngrid,"co2ice","CO2 ice cover",
2363	& "kg.m-2",2,co2ice)
2364	call wstats(ngrid,"tauref","reference dod at 610 Pa","NU",
2365	& 2,tauref)
2366	call wstats(ngrid,"fluxsurf_lw",
2367	& "Thermal IR radiative flux to surface","W.m-2",2,
2368	& fluxsurf_lw)
2369	call wstats(ngrid,"fluxsurf_sw",
2370	& "Solar radiative flux to surface","W.m-2",2,
2371	& fluxsurf_sw_tot)
2372	call wstats(ngrid,"fluxtop_lw",
2373	& "Thermal IR radiative flux to space","W.m-2",2,
2374	& fluxtop_lw)
2375	call wstats(ngrid,"fluxtop_sw",
2376	& "Solar radiative flux to space","W.m-2",2,
2377	& fluxtop_sw_tot)
2378	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
2379	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
2380	call wstats(ngrid,"v","Meridional (North-South) wind",
2381	& "m.s-1",3,zv)
2382	call wstats(ngrid,"w","Vertical (down-up) wind",
2383	& "m.s-1",3,pw)
2384	call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
2385	call wstats(ngrid,"pressure","Pressure","Pa",3,zplay)
2386	call wstats(ngrid,"q2",
2387	& "Boundary layer eddy kinetic energy",
2388	& "m2.s-2",3,q2)
2389	call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
2390	& emis)
2391	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
2392	c & 2,zstress)
2393	c call wstats(ngrid,"sw_htrt","sw heat.rate",
2394	c & "W.m-2",3,zdtsw)
2395	c call wstats(ngrid,"lw_htrt","lw heat.rate",
2396	c & "W.m-2",3,zdtlw)
2397
2398	if (calltherm) then
2399	call wstats(ngrid,"zmax_th","Height of thermals",
2400	& "m",2,zmax_th)
2401	call wstats(ngrid,"hfmax_th","Max thermals heat flux",
2402	& "K.m/s",2,hfmax_th)
2403	call wstats(ngrid,"wstar",
2404	& "Max vertical velocity in thermals",
2405	& "m/s",2,wstar)
2406	endif
2407
2408	if (tracer) then
2409	if (water) then
2410	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
2411	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2412	call wstats(ngrid,"vmr_h2ovap",
2413	& "H2O vapor volume mixing ratio","mol/mol",
2414	& 3,vmr)
2415	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
2416	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
2417	call wstats(ngrid,"vmr_h2oice",
2418	& "H2O ice volume mixing ratio","mol/mol",
2419	& 3,vmr)
2420	! also store vmr_ice*rice for better diagnostics of rice
2421	vmr(1:ngrid,1:nlayer)=vmr(1:ngrid,1:nlayer)*
2422	& rice(1:ngrid,1:nlayer)
2423	call wstats(ngrid,"vmr_h2oice_rice",
2424	& "H2O ice mixing ratio times ice particule size",
2425	& "(mol/mol)*m",
2426	& 3,vmr)
2427	vmr=zqsat(1:ngrid,1:nlayer)
2428	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2429	call wstats(ngrid,"vmr_h2osat",
2430	& "saturation volume mixing ratio","mol/mol",
2431	& 3,vmr)
2432	call wstats(ngrid,"h2o_ice_s",
2433	& "surface h2o_ice","kg/m2",
2434	& 2,qsurf(1,igcm_h2o_ice))
2435	call wstats(ngrid,'albedo',
2436	& 'albedo',
2437	& '',2,albedo(1,1))
2438	call wstats(ngrid,"mtot",
2439	& "total mass of water vapor","kg/m2",
2440	& 2,mtot)
2441	call wstats(ngrid,"icetot",
2442	& "total mass of water ice","kg/m2",
2443	& 2,icetot)
2444	call wstats(ngrid,"reffice",
2445	& "Mean reff","m",
2446	& 2,rave)
2447	call wstats(ngrid,"Nccntot",
2448	& "condensation nuclei","Nbr/m2",
2449	& 2,Nccntot)
2450	call wstats(ngrid,"Mccntot",
2451	& "condensation nuclei mass","kg/m2",
2452	& 2,Mccntot)
2453	call wstats(ngrid,"rice",
2454	& "Ice particle size","m",
2455	& 3,rice)
2456	if (.not.activice) then
2457	call wstats(ngrid,"tauTESap",
2458	& "tau abs 825 cm-1","",
2459	& 2,tauTES)
2460	else
2461	call wstats(ngrid,'tauTES',
2462	& 'tau abs 825 cm-1',
2463	& '',2,taucloudtes)
2464	endif
2465
2466	endif ! of if (water)
2467
2468
2469	if (dustbin.ne.0) then
2470
2471	call wstats(ngrid,'tau','taudust','SI',2,tau(1,1))
2472
2473	if (doubleq) then
2474	c call wstats(ngrid,'qsurf','qsurf',
2475	c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))
2476	c call wstats(ngrid,'Nsurf','N particles',
2477	c & 'N.m-2',2,qsurf(1,igcm_dust_number))
2478	c call wstats(ngrid,'dqsdev','ddevil lift',
2479	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
2480	c call wstats(ngrid,'dqssed','sedimentation',
2481	c & 'kg.m-2.s-1',2,zdqssed(1,1))
2482	c call wstats(ngrid,'dqsdif','diffusion',
2483	c & 'kg.m-2.s-1',2,zdqsdif(1,1))
2484	call wstats(ngrid,'dqsdust',
2485	& 'deposited surface dust mass',
2486	& 'kg.m-2.s-1',2,dqdustsurf)
2487	call wstats(ngrid,'dqndust',
2488	& 'deposited surface dust number',
2489	& 'number.m-2.s-1',2,dndustsurf)
2490	call wstats(ngrid,'reffdust','reffdust',
2491	& 'm',3,rdust*ref_r0)
2492	call wstats(ngrid,'dustq','Dust mass mr',
2493	& 'kg/kg',3,qdust)
2494	call wstats(ngrid,'dustN','Dust number',
2495	& 'part/kg',3,ndust)
2496	if (rdstorm) then
2497	call wstats(ngrid,'reffstormdust','reffdust',
2498	& 'm',3,rstormdust*ref_r0)
2499	call wstats(ngrid,'rdsdustq','Dust mass mr',
2500	& 'kg/kg',3,rdsqdust)
2501	call wstats(ngrid,'rdsdustN','Dust number',
2502	& 'part/kg',3,rdsndust)
2503	end if
2504	else
2505	do iq=1,dustbin
2506	write(str2(1:2),'(i2.2)') iq
2507	call wstats(ngrid,'q'//str2,'mix. ratio',
2508	& 'kg/kg',3,zq(1,1,iq))
2509	call wstats(ngrid,'qsurf'//str2,'qsurf',
2510	& 'kg.m-2',2,qsurf(1,iq))
2511	end do
2512	endif ! (doubleq)
2513
2514	if (scavenging) then
2515	call wstats(ngrid,'ccnq','CCN mass mr',
2516	& 'kg/kg',3,qccn)
2517	call wstats(ngrid,'ccnN','CCN number',
2518	& 'part/kg',3,nccn)
2519	endif ! (scavenging)
2520
2521	endif ! (dustbin.ne.0)
2522
2523	if (photochem) then
2524	do iq=1,nq
2525	if (noms(iq) .ne. "dust_mass" .and.
2526	$ noms(iq) .ne. "dust_number" .and.
2527	$ noms(iq) .ne. "ccn_mass" .and.
2528	$ noms(iq) .ne. "ccn_number" .and.
2529	$ noms(iq) .ne. "ccnco2_mass" .and.
2530	$ noms(iq) .ne. "ccnco2_number") then
2531
2532	! volume mixing ratio
2533
2534	vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
2535	& *mmean(1:ngrid,1:nlayer)/mmol(iq)
2536
2537	call wstats(ngrid,"vmr_"//trim(noms(iq)),
2538	$ "Volume mixing ratio","mol/mol",3,vmr)
2539	if ((noms(iq).eq."o")
2540	$ .or. (noms(iq).eq."co2")
2541	$ .or. (noms(iq).eq."o3")
2542	$ .or. (noms(iq).eq."ar")
2543	$ .or. (noms(iq).eq."o2")
2544	$ .or. (noms(iq).eq."h2o_vap") ) then
2545	call writediagfi(ngrid,"vmr_"//trim(noms(iq)),
2546	$ "Volume mixing ratio","mol/mol",3,vmr)
2547	end if
2548
2549	! number density (molecule.cm-3)
2550
2551	rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
2552	& rho(1:ngrid,1:nlayer)n_avog/
2553	& (1000*mmol(iq))
2554
2555	call wstats(ngrid,"rho_"//trim(noms(iq)),
2556	$ "Number density","cm-3",3,rhopart)
2557	call writediagfi(ngrid,"rho_"//trim(noms(iq)),
2558	$ "Number density","cm-3",3,rhopart)
2559
2560	! vertical column (molecule.cm-2)
2561
2562	do ig = 1,ngrid
2563	colden(ig,iq) = 0.
2564	end do
2565	do l=1,nlayer
2566	do ig=1,ngrid
2567	colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
2568	$ *(zplev(ig,l)-zplev(ig,l+1))
2569	$ 6.022e22/(mmol(iq)g)
2570	end do
2571	end do
2572
2573	call wstats(ngrid,"c_"//trim(noms(iq)),
2574	$ "column","mol cm-2",2,colden(1,iq))
2575	call writediagfi(ngrid,"c_"//trim(noms(iq)),
2576	$ "column","mol cm-2",2,colden(1,iq))
2577
2578	! global mass (g)
2579
2580	call planetwide_sumval(colden(:,iq)/6.022e23
2581	$ mmol(iq)1.e4*cell_area(:),mass(iq))
2582
2583	call writediagfi(ngrid,"mass_"//trim(noms(iq)),
2584	$ "global mass","g",0,mass(iq))
2585
2586	end if ! of if (noms(iq) .ne. "dust_mass" ...)
2587	end do ! of do iq=1,nq
2588	end if ! of if (photochem)
2589
2590	end if ! of if (tracer)
2591
2592	IF(lastcall) THEN
2593	write (,) "Writing stats..."
2594	call mkstats(ierr)
2595	ENDIF
2596
2597	ENDIF !if callstats
2598
2599	c (Store EOF for Mars Climate database software)
2600	IF (calleofdump) THEN
2601	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
2602	ENDIF
2603	#endif
2604	!endif of ifndef MESOSCALE
2605
2606	#ifdef MESOSCALE
2607
2608	!! see comm_wrf.
2609	!! not needed when an array is already in a shared module.
2610	!! --> example : hfmax_th, zmax_th
2611
2612	CALL allocate_comm_wrf(ngrid,nlayer)
2613
2614	!state real HR_SW ikj misc 1 - h "HR_SW" "HEATING RATE SW" "K/s"
2615	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
2616	!state real HR_LW ikj misc 1 - h "HR_LW" "HEATING RATE LW" "K/s"
2617	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
2618	!state real SWDOWNZ ij misc 1 - h "SWDOWNZ" "DOWNWARD SW FLUX AT SURFACE" "W m-2"
2619	comm_SWDOWNZ(1:ngrid) = fluxsurf_sw_tot(1:ngrid)
2620	!state real TAU_DUST ij misc 1 - h "TAU_DUST" "REFERENCE VISIBLE DUST OPACITY" ""
2621	comm_TAU_DUST(1:ngrid) = tauref(1:ngrid)
2622	!state real RDUST ikj misc 1 - h "RDUST" "DUST RADIUS" "m"
2623	comm_RDUST(1:ngrid,1:nlayer) = rdust(1:ngrid,1:nlayer)
2624	!state real QSURFDUST ij misc 1 - h "QSURFDUST" "DUST MASS AT SURFACE" "kg m-2"
2625	IF (igcm_dust_mass .ne. 0) THEN
2626	comm_QSURFDUST(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
2627	ELSE
2628	comm_QSURFDUST(1:ngrid) = 0.
2629	ENDIF
2630	!state real MTOT ij misc 1 - h "MTOT" "TOTAL MASS WATER VAPOR in pmic" "pmic"
2631	comm_MTOT(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
2632	!state real ICETOT ij misc 1 - h "ICETOT" "TOTAL MASS WATER ICE" "kg m-2"
2633	comm_ICETOT(1:ngrid) = icetot(1:ngrid) * 1.e6 / rho_ice
2634	!state real VMR_ICE ikj misc 1 - h "VMR_ICE" "VOL. MIXING RATIO ICE" "ppm"
2635	IF (igcm_h2o_ice .ne. 0) THEN
2636	comm_VMR_ICE(1:ngrid,1:nlayer) = 1.e6
2637	. * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
2638	. * mmean(1:ngrid,1:nlayer) / mmol(igcm_h2o_ice)
2639	ELSE
2640	comm_VMR_ICE(1:ngrid,1:nlayer) = 0.
2641	ENDIF
2642	!state real TAU_ICE ij misc 1 - h "TAU_ICE" "CLOUD OD at 825 cm-1 TES" ""
2643	if (activice) then
2644	comm_TAU_ICE(1:ngrid) = taucloudtes(1:ngrid)
2645	else
2646	comm_TAU_ICE(1:ngrid) = tauTES(1:ngrid)
2647	endif
2648	!state real RICE ikj misc 1 - h "RICE" "ICE RADIUS" "m"
2649	comm_RICE(1:ngrid,1:nlayer) = rice(1:ngrid,1:nlayer)
2650
2651	!! calculate sensible heat flux in W/m2 for outputs
2652	!! -- the one computed in vdifc is not the real one
2653	!! -- vdifc must have been called
2654	if (.not.callrichsl) then
2655	sensibFlux(1:ngrid) = zflubid(1:ngrid)
2656	. - capcal(1:ngrid)*zdtsdif(1:ngrid)
2657	else
2658	sensibFlux(1:ngrid) =
2659	& (pplay(1:ngrid,1)/(rpt(1:ngrid,1)))cpp
2660	& sqrt(pu(1:ngrid,1)pu(1:ngrid,1)+pv(1:ngrid,1)*pv(1:ngrid,1)
2661	& +(log(1.+0.7wstar(1:ngrid) + 2.3wstar(1:ngrid)2))2)
2662	& zcdh(1:ngrid)(tsurf(1:ngrid)-zh(1:ngrid,1))
2663	endif
2664
2665	#else
2666	#ifndef MESOINI
2667
2668	c ==========================================================
2669	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
2670	c any variable for diagnostic (output with period
2671	c "ecritphy", set in "run.def")
2672	c ==========================================================
2673	c WRITEDIAGFI can ALSO be called from any other subroutines
2674	c for any variables !!
2675	c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
2676	c & emis)
2677	c call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
2678	c call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
2679	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
2680	& tsurf)
2681	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
2682	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
2683	& ,"kg.m-2",2,co2ice)
2684
2685	call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
2686	& "K",2,zt(1,7))
2687	call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
2688	& fluxsurf_lw)
2689	call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
2690	& fluxsurf_sw_tot)
2691	call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
2692	& fluxtop_lw)
2693	call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
2694	& fluxtop_sw_tot)
2695	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
2696	call WRITEDIAGFI(ngrid,"Sols","Time","sols",0,zday)
2697	call WRITEDIAGFI(ngrid,"Ls","Solar longitude","deg",
2698	& 0,zls*180./pi)
2699	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
2700	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
2701	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
2702	call WRITEDIAGFI(ngrid,"rho","density","kg.m-3",3,rho)
2703	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
2704	c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
2705	call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,zplay)
2706	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
2707	c & zstress)
2708	c call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
2709	c & 'w.m-2',3,zdtsw)
2710	c call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
2711	c & 'w.m-2',3,zdtlw)
2712
2713	if (.not.activice) then
2714	CALL WRITEDIAGFI(ngrid,'tauTESap',
2715	& 'tau abs 825 cm-1',
2716	& '',2,tauTES)
2717	else
2718	CALL WRITEDIAGFI(ngrid,'tauTES',
2719	& 'tau abs 825 cm-1',
2720	& '',2,taucloudtes)
2721	endif
2722
2723	#else
2724	!!! this is to ensure correct initialisation of mesoscale model
2725	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
2726	& tsurf)
2727	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
2728	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
2729	& co2ice)
2730	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
2731	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
2732	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
2733	call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
2734	& emis)
2735	call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
2736	& "K",3,tsoil)
2737	call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
2738	& "K",3,inertiedat)
2739	#endif
2740
2741
2742	c ----------------------------------------------------------
2743	c Outputs of the CO2 cycle
2744	c ----------------------------------------------------------
2745
2746	if (tracer.and.(igcm_co2.ne.0)) then
2747	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
2748	! & "kg/kg",2,zq(1,1,igcm_co2))
2749	call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
2750	& "kg/kg",3,zq(1,1,igcm_co2))
2751	if (co2clouds) then
2752	CALL WRITEDIAGFI(ngrid,'mtotco2',
2753	& 'total mass of CO2 vapor',
2754	& 'kg/m2',2,mtotco2)
2755	CALL WRITEDIAGFI(ngrid,'zdtcloudco2',
2756	& 'temperature variation of CO2 latent heat',
2757	& 'K/s',3,zdtcloudco2)
2758
2759	CALL WRITEDIAGFI(ngrid,'icetotco2',
2760	& 'total mass of CO2 ice',
2761	& 'kg/m2',2,icetotco2)
2762
2763	call WRITEDIAGFI(ngrid,'ccnqco2','CCNco2 mass mr',
2764	& 'kg/kg',3,qccnco2)
2765	call WRITEDIAGFI(ngrid,'ccnNco2','CCNco2 number',
2766	& 'part/kg',3,nccnco2)
2767	call WRITEDIAGFI(ngrid,'co2_ice','co2_ice','kg/kg',
2768	& 3,zq(:,:,igcm_co2_ice))
2769	call WRITEDIAGFI(ngrid,'precip_co2_ice',
2770	& 'surface deposition of co2 ice',
2771	& 'kg.m-2.s-1',2,
2772	& zdqssed(1:ngrid,igcm_co2_ice))
2773	endif ! of if (co2clouds)
2774	endif ! of if (tracer.and.(igcm_co2.ne.0))
2775	! Output He tracer, if there is one
2776	if (tracer.and.(igcm_he.ne.0)) then
2777	call WRITEDIAGFI(ngrid,"he","helium mass mixing ratio",
2778	& "kg/kg",3,zq(1,1,igcm_he))
2779	vmr=zq(1:ngrid,1:nlayer,igcm_he)
2780	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_he)
2781	call WRITEDIAGFI(ngrid,'vmr_he','helium vol. mixing ratio',
2782	& 'mol/mol',3,vmr)
2783	endif
2784
2785	c ----------------------------------------------------------
2786	c Outputs of the water cycle
2787	c ----------------------------------------------------------
2788	if (tracer) then
2789	if (water) then
2790
2791	#ifdef MESOINI
2792	!!!! waterice = q01, voir readmeteo.F90
2793	call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice),
2794	& 'kg/kg',3,
2795	& zq(1:ngrid,1:nlayer,igcm_h2o_ice))
2796	!!!! watervapor = q02, voir readmeteo.F90
2797	call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap),
2798	& 'kg/kg',3,
2799	& zq(1:ngrid,1:nlayer,igcm_h2o_vap))
2800	!!!! surface waterice qsurf02 (voir readmeteo)
2801	call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
2802	& 'kg.m-2',2,
2803	& qsurf(1:ngrid,igcm_h2o_ice))
2804	#endif
2805	CALL WRITEDIAGFI(ngrid,'mtot',
2806	& 'total mass of water vapor',
2807	& 'kg/m2',2,mtot)
2808	CALL WRITEDIAGFI(ngrid,'icetot',
2809	& 'total mass of water ice',
2810	& 'kg/m2',2,icetot)
2811	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
2812	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
2813	CALL WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
2814	& 'mol/mol',3,vmr)
2815	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
2816	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2817	CALL WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
2818	& 'mol/mol',3,vmr)
2819	CALL WRITEDIAGFI(ngrid,'reffice',
2820	& 'Mean reff',
2821	& 'm',2,rave)
2822
2823	call WRITEDIAGFI(ngrid,'h2o_ice','h2o_ice','kg/kg',
2824	& 3,zq(:,:,igcm_h2o_ice))
2825	call WRITEDIAGFI(ngrid,'h2o_vap','h2o_vap','kg/kg',
2826	& 3,zq(:,:,igcm_h2o_vap))
2827
2828
2829	!A. Pottier
2830	CALL WRITEDIAGFI(ngrid,'rmoym',
2831	& 'alternative reffice',
2832	& 'm',2,rave2)
2833	call WRITEDIAGFI(ngrid,'saturation',
2834	& 'h2o vap saturation ratio','dimless',3,satu)
2835	if (scavenging) then
2836	CALL WRITEDIAGFI(ngrid,"Nccntot",
2837	& "condensation nuclei","Nbr/m2",
2838	& 2,Nccntot)
2839	CALL WRITEDIAGFI(ngrid,"Mccntot",
2840	& "mass condensation nuclei","kg/m2",
2841	& 2,Mccntot)
2842	endif
2843	call WRITEDIAGFI(ngrid,'rice','Ice particle size',
2844	& 'm',3,rice)
2845	call WRITEDIAGFI(ngrid,'h2o_ice_s',
2846	& 'surface h2o_ice',
2847	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
2848	CALL WRITEDIAGFI(ngrid,'albedo',
2849	& 'albedo',
2850	& '',2,albedo(1,1))
2851	if (tifeedback) then
2852	call WRITEDIAGSOIL(ngrid,"soiltemp",
2853	& "Soil temperature","K",
2854	& 3,tsoil)
2855	call WRITEDIAGSOIL(ngrid,'soilti',
2856	& 'Soil Thermal Inertia',
2857	& 'J.s-1/2.m-2.K-1',3,inertiesoil)
2858	endif
2859	!A. Pottier
2860	if (CLFvarying) then !AP14 nebulosity
2861	call WRITEDIAGFI(ngrid,'totcloudfrac',
2862	& 'Total cloud fraction',
2863	& ' ',2,totcloudfrac)
2864	endif !clf varying
2865
2866	endif !(water)
2867
2868
2869	if (water.and..not.photochem) then
2870	iq=nq
2871	c write(str2(1:2),'(i2.2)') iq
2872	c call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
2873	c & 'kg.m-2',2,zdqscloud(1,iq))
2874	c call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
2875	c & 'kg/kg',3,zdqchim(1,1,iq))
2876	c call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
2877	c & 'kg/kg',3,zdqdif(1,1,iq))
2878	c call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
2879	c & 'kg/kg',3,zdqadj(1,1,iq))
2880	c call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
2881	c & 'kg/kg',3,zdqc(1,1,iq))
2882	endif !(water.and..not.photochem)
2883	endif
2884
2885	c ----------------------------------------------------------
2886	c Outputs of the dust cycle
2887	c ----------------------------------------------------------
2888
2889	call WRITEDIAGFI(ngrid,'tauref',
2890	& 'Dust ref opt depth','NU',2,tauref)
2891
2892	if (tracer.and.(dustbin.ne.0)) then
2893
2894	call WRITEDIAGFI(ngrid,'tau','taudust','SI',2,tau(1,1))
2895
2896	#ifndef MESOINI
2897	if (doubleq) then
2898	c call WRITEDIAGFI(ngrid,'qsurf','qsurf',
2899	c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))
2900	c call WRITEDIAGFI(ngrid,'Nsurf','N particles',
2901	c & 'N.m-2',2,qsurf(1,igcm_dust_number))
2902	c call WRITEDIAGFI(ngrid,'dqsdev','ddevil lift',
2903	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
2904	c call WRITEDIAGFI(ngrid,'dqssed','sedimentation',
2905	c & 'kg.m-2.s-1',2,zdqssed(1,1))
2906	c call WRITEDIAGFI(ngrid,'dqsdif','diffusion',
2907	c & 'kg.m-2.s-1',2,zdqsdif(1,1))
2908	c call WRITEDIAGFI(ngrid,'sedice','sedimented ice',
2909	c & 'kg.m-2.s-1',2,zdqssed(:,igcm_h2o_ice))
2910	c call WRITEDIAGFI(ngrid,'subice','sublimated ice',
2911	c & 'kg.m-2.s-1',2,zdqsdif(:,igcm_h2o_ice))
2912	call WRITEDIAGFI(ngrid,'dqsdust',
2913	& 'deposited surface dust mass',
2914	& 'kg.m-2.s-1',2,dqdustsurf)
2915	call WRITEDIAGFI(ngrid,'dqndust',
2916	& 'deposited surface dust number',
2917	& 'number.m-2.s-1',2,dndustsurf)
2918	call WRITEDIAGFI(ngrid,'reffdust','reffdust',
2919	& 'm',3,rdust*ref_r0)
2920	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2921	& 'kg/kg',3,qdust)
2922	call WRITEDIAGFI(ngrid,'dustN','Dust number',
2923	& 'part/kg',3,ndust)
2924	call WRITEDIAGFI(ngrid,'dsodust',
2925	& 'density scaled optical depth in the IR',
2926	& 'm2.kg-1',3,dsodust)
2927	call WRITEDIAGFI(ngrid,'dso',
2928	& 'density scaled optical depth in the IR',
2929	& 'm2.kg-1',3,dsodust+dsords)
2930	else
2931	do iq=1,dustbin
2932	write(str2(1:2),'(i2.2)') iq
2933	call WRITEDIAGFI(ngrid,'q'//str2,'mix. ratio',
2934	& 'kg/kg',3,zq(1,1,iq))
2935	call WRITEDIAGFI(ngrid,'qsurf'//str2,'qsurf',
2936	& 'kg.m-2',2,qsurf(1,iq))
2937	end do
2938	endif ! (doubleq)
2939
2940	if (rdstorm) then ! writediagfi tendencies stormdust tracers
2941	call WRITEDIAGFI(ngrid,'reffdust','reffdust',
2942	& 'm',3,rdust*ref_r0)
2943	call WRITEDIAGFI(ngrid,'reffstormdust','reffstormdust',
2944	& 'm',3,rstormdust*ref_r0)
2945	call WRITEDIAGFI(ngrid,'mstormdtot',
2946	& 'total mass of stormdust only',
2947	& 'kg.m-2',2,mstormdtot)
2948	call WRITEDIAGFI(ngrid,'mdusttot',
2949	& 'total mass of dust only',
2950	& 'kg.m-2',2,mdusttot)
2951	call WRITEDIAGFI(ngrid,'rdsdqsdust',
2952	& 'deposited surface stormdust mass',
2953	& 'kg.m-2.s-1',2,rdsdqdustsurf)
2954	call WRITEDIAGFI(ngrid,'rdsdustq','storm Dust mass mr',
2955	& 'kg/kg',3,rdsqdust)
2956	call WRITEDIAGFI(ngrid,'rdsdustqmodel','storm Dust massmr',
2957	& 'kg/kg',3,pq(:,:,igcm_stormdust_mass))
2958	call WRITEDIAGFI(ngrid,'rdsdustN','storm Dust number',
2959	& 'part/kg',3,rdsndust)
2960	call WRITEDIAGFI(ngrid,"stormfract",
2961	& "fraction of the mesh, with stormdust","none",
2962	& 2,totstormfract)
2963	call WRITEDIAGFI(ngrid,'qsurf',
2964	& 'stormdust injection',
2965	& 'kg.m-2',2,qsurf(:,igcm_stormdust_mass))
2966	call WRITEDIAGFI(ngrid,'pdqsurf',
2967	& 'tendancy stormdust mass at surface',
2968	& 'kg.m-2',2,dqsurf(:,igcm_stormdust_mass))
2969	call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust',
2970	& 'm/s',3,wspeed(:,1:nlayer))
2971	call WRITEDIAGFI(ngrid,'zdqsed_dustq'
2972	& ,'sedimentation q','kg.m-2.s-1',3,
2973	& zdqsed(:,:,igcm_dust_mass))
2974	call WRITEDIAGFI(ngrid,'zdqssed_dustq'
2975	& ,'sedimentation q','kg.m-2.s-1',2,
2976	& zdqssed(:,igcm_dust_mass))
2977	call WRITEDIAGFI(ngrid,'zdqsed_rdsq'
2978	& ,'sedimentation q','kg.m-2.s-1',3,
2979	& zdqsed(:,:,igcm_stormdust_mass))
2980	call WRITEDIAGFI(ngrid,'rdust','rdust',
2981	& 'm',3,rdust)
2982	call WRITEDIAGFI(ngrid,'rstormdust','rstormdust',
2983	& 'm',3,rstormdust)
2984	call WRITEDIAGFI(ngrid,'totaldustq','total dust mass',
2985	& 'kg/kg',3,qdusttotal)
2986	call WRITEDIAGFI(ngrid,'dsords',
2987	& 'density scaled opacity of stormdust',
2988	& 'm2.kg-1',3,dsords)
2989	endif ! (rdstorm)
2990
2991	if (scavenging) then
2992	call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr',
2993	& 'kg/kg',3,qccn)
2994	call WRITEDIAGFI(ngrid,'ccnN','CCN number',
2995	& 'part/kg',3,nccn)
2996	call WRITEDIAGFI(ngrid,'surfccnq','Surf nuclei mass mr',
2997	& 'kg.m-2',2,qsurf(1,igcm_ccn_mass))
2998	call WRITEDIAGFI(ngrid,'surfccnN','Surf nuclei number',
2999	& 'kg.m-2',2,qsurf(1,igcm_ccn_number))
3000	endif ! (scavenging)
3001
3002	c if (submicron) then
3003	c call WRITEDIAGFI(ngrid,'dustsubm','subm mass mr',
3004	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
3005	c endif ! (submicron)
3006
3007	#else
3008	! !!! to initialize mesoscale we need scaled variables
3009	! !!! because this must correspond to starting point for tracers
3010	! call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
3011	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass))
3012	! call WRITEDIAGFI(ngrid,'dustN','Dust number',
3013	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number))
3014	! call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
3015	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass))
3016	! call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
3017	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number))
3018	if (freedust) then
3019	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
3020	& 'kg/kg',3,qdust)
3021	call WRITEDIAGFI(ngrid,'dustN','Dust number',
3022	& 'part/kg',3,ndust)
3023	call WRITEDIAGFI(ngrid,'ccn','CCN mass mr',
3024	& 'kg/kg',3,qccn)
3025	call WRITEDIAGFI(ngrid,'ccnN','CCN number',
3026	& 'part/kg',3,nccn)
3027	else
3028	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
3029	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
3030	call WRITEDIAGFI(ngrid,'dustN','Dust number',
3031	& 'part/kg',3,pq(1,1,igcm_dust_number))
3032	call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
3033	& 'kg/kg',3,pq(1,1,igcm_ccn_mass))
3034	call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
3035	& 'part/kg',3,pq(1,1,igcm_ccn_number))
3036	endif
3037	#endif
3038
3039	end if ! (tracer.and.(dustbin.ne.0))
3040
3041	c ----------------------------------------------------------
3042	c GW non-oro outputs
3043	c ----------------------------------------------------------
3044
3045	if(calllott_nonoro) then
3046	call WRITEDIAGFI(ngrid,"dugwno","GW non-oro dU","m/s2",
3047	$ 3,d_u_hin/ptimestep)
3048	call WRITEDIAGFI(ngrid,"dvgwno","GW non-oro dV","m/s2",
3049	$ 3,d_v_hin/ptimestep)
3050	endif !(calllott_nonoro)
3051
3052	c ----------------------------------------------------------
3053	c Thermospheric outputs
3054	c ----------------------------------------------------------
3055
3056	if(callthermos) then
3057
3058	call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s",
3059	$ 3,zdtnlte)
3060	call WRITEDIAGFI(ngrid,"quv","UV heating","K/s",
3061	$ 3,zdteuv)
3062	call WRITEDIAGFI(ngrid,"cond","Thermal conduction","K/s",
3063	$ 3,zdtconduc)
3064	call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s",
3065	$ 3,zdtnirco2)
3066
3067	! call wstats(ngrid,"q15um","15 um cooling","K/s",
3068	! $ 3,zdtnlte)
3069	! call wstats(ngrid,"quv","UV heating","K/s",
3070	! $ 3,zdteuv)
3071	! call wstats(ngrid,"cond","Thermal conduction","K/s",
3072	! $ 3,zdtconduc)
3073	! call wstats(ngrid,"qnir","NIR heating","K/s",
3074	! $ 3,zdtnirco2)
3075
3076	endif !(callthermos)
3077
3078	c ----------------------------------------------------------
3079	c ----------------------------------------------------------
3080	c PBL OUTPUS
3081	c ----------------------------------------------------------
3082	c ----------------------------------------------------------
3083
3084	c ----------------------------------------------------------
3085	c Outputs of thermals
3086	c ----------------------------------------------------------
3087	if (calltherm) then
3088
3089	! call WRITEDIAGFI(ngrid,'dtke',
3090	! & 'tendance tke thermiques','m2/s2',
3091	! & 3,dtke_th)
3092	! call WRITEDIAGFI(ngrid,'d_u_ajs',
3093	! & 'tendance u thermiques','m/s',
3094	! & 3,pdu_th*ptimestep)
3095	! call WRITEDIAGFI(ngrid,'d_v_ajs',
3096	! & 'tendance v thermiques','m/s',
3097	! & 3,pdv_th*ptimestep)
3098	! if (tracer) then
3099	! if (nq .eq. 2) then
3100	! call WRITEDIAGFI(ngrid,'deltaq_th',
3101	! & 'delta q thermiques','kg/kg',
3102	! & 3,ptimestep*pdq_th(:,:,2))
3103	! endif
3104	! endif
3105
3106	call WRITEDIAGFI(ngrid,'zmax_th',
3107	& 'hauteur du thermique','m',
3108	& 2,zmax_th)
3109	call WRITEDIAGFI(ngrid,'hfmax_th',
3110	& 'maximum TH heat flux','K.m/s',
3111	& 2,hfmax_th)
3112	call WRITEDIAGFI(ngrid,'wstar',
3113	& 'maximum TH vertical velocity','m/s',
3114	& 2,wstar)
3115
3116	endif
3117
3118	c ----------------------------------------------------------
3119	c ----------------------------------------------------------
3120	c END OF PBL OUTPUS
3121	c ----------------------------------------------------------
3122	c ----------------------------------------------------------
3123
3124
3125	c ----------------------------------------------------------
3126	c Output in netcdf file "diagsoil.nc" for subterranean
3127	c variables (output every "ecritphy", as for writediagfi)
3128	c ----------------------------------------------------------
3129
3130	! Write soil temperature
3131	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
3132	! & 3,tsoil)
3133	! Write surface temperature
3134	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
3135	! & 2,tsurf)
3136
3137	c ==========================================================
3138	c END OF WRITEDIAGFI
3139	c ==========================================================
3140	#endif
3141	! of ifdef MESOSCALE
3142
3143	ELSE ! if(ngrid.eq.1)
3144
3145	#ifndef MESOSCALE
3146	write(*,'("Ls =",f11.6," tauref(",f4.0," Pa) =",f9.6)')
3147	& zls*180./pi,odpref,tauref
3148	c ----------------------------------------------------------------------
3149	c Output in grads file "g1d" (ONLY when using testphys1d)
3150	c (output at every X physical timestep)
3151	c ----------------------------------------------------------------------
3152	c
3153	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
3154	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
3155	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
3156
3157	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
3158	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
3159	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
3160	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
3161
3162	! THERMALS STUFF 1D
3163	if(calltherm) then
3164
3165	call WRITEDIAGFI(ngrid,'lmax_th',
3166	& 'hauteur du thermique','point',
3167	& 0,lmax_th_out)
3168	call WRITEDIAGFI(ngrid,'zmax_th',
3169	& 'hauteur du thermique','m',
3170	& 0,zmax_th)
3171	call WRITEDIAGFI(ngrid,'hfmax_th',
3172	& 'maximum TH heat flux','K.m/s',
3173	& 0,hfmax_th)
3174	call WRITEDIAGFI(ngrid,'wstar',
3175	& 'maximum TH vertical velocity','m/s',
3176	& 0,wstar)
3177
3178	co2col(:)=0.
3179	if (tracer) then
3180	do l=1,nlayer
3181	do ig=1,ngrid
3182	co2col(ig)=co2col(ig)+
3183	& zq(ig,l,1)*(zplev(ig,l)-zplev(ig,l+1))/g
3184	enddo
3185	enddo
3186
3187	end if
3188	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
3189	& ,'kg/m-2',0,co2col)
3190	endif ! of if (calltherm)
3191
3192	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
3193	& ,'m/s',1,pw)
3194	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
3195	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
3196	& tsurf)
3197	call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
3198	call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
3199
3200	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
3201	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
3202	call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
3203	call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate", &
3204	& "K.s-1",1,dtrad)
3205	call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
3206	& 'w.m-2',1,zdtsw)
3207	call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
3208	& 'w.m-2',1,zdtlw)
3209	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
3210	& ,"kg.m-2",0,co2ice)
3211
3212	if (igcm_co2.ne.0) then
3213	call co2sat(ngrid*nlayer,zt,zplay,zqsatco2)
3214	do ig=1,ngrid
3215	do l=1,nlayer
3216	satuco2(ig,l) = zq(ig,l,igcm_co2)*
3217	& (mmean(ig,l)/44.01)*zplay(ig,l)/zqsatco2(ig,l)
3218
3219	c write(,) "In PHYSIQMOD, pt,zt,time ",pt(ig,l)
3220	c & ,zt(ig,l),ptime
3221	enddo
3222	enddo
3223	endif
3224
3225	c CALL writeg1d(ngrid,nlayer,zt,'temp','K')
3226	c CALL writeg1d(ngrid,nlayer,riceco2,'riceco2','m')
3227	c CALL writeg1d(ngrid,nlayer,satuco2,'satuco2','satu')
3228
3229
3230	c call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",1,
3231	c & satuco2)
3232	c call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
3233	c & ,1,riceco2)
3234	! or output in diagfi.nc (for testphys1d)
3235	call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
3236	call WRITEDIAGFI(ngrid,'temp','Temperature ',
3237	& 'K JA',1,zt)
3238	c call WRITEDIAGFI(ngrid,'temp2','Temperature ',
3239	c & 'K JA2',1,pt)
3240
3241	if(tracer) then
3242	c CALL writeg1d(ngrid,1,tau,'tau','SI')
3243	do iq=1,nq
3244	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
3245	call WRITEDIAGFI(ngrid,trim(noms(iq)),
3246	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
3247	end do
3248	if (doubleq) then
3249	call WRITEDIAGFI(ngrid,'rdust','rdust',
3250	& 'm',1,rdust)
3251	endif ! doubleq 1D
3252	if (rdstorm) then
3253	call writediagfi(1,'aerosol_dust','opacity of env. dust',''
3254	& ,1,aerosol(:,:,iaer_dust_doubleq))
3255	call writediagfi(1,'aerosol_stormdust',
3256	& 'opacity of storm dust',''
3257	& ,1,aerosol(:,:,iaer_stormdust_doubleq))
3258	call WRITEDIAGFI(ngrid,'dqsdifdustq','diffusion',
3259	& 'kg.m-2.s-1',0,zdqsdif(1,igcm_dust_mass))
3260	call WRITEDIAGFI(ngrid,'dqsdifrdsq','diffusion',
3261	& 'kg.m-2.s-1',0,zdqsdif(1,igcm_stormdust_mass))
3262	call WRITEDIAGFI(ngrid,'mstormdtot',
3263	& 'total mass of stormdust only',
3264	& 'kg.m-2',0,mstormdtot)
3265	call WRITEDIAGFI(ngrid,'mdusttot',
3266	& 'total mass of dust only',
3267	& 'kg.m-2',0,mdusttot)
3268	call WRITEDIAGFI(ngrid,'tauref',
3269	& 'Dust ref opt depth','NU',0,tauref)
3270	call WRITEDIAGFI(ngrid,'rdsdqsdust',
3271	& 'deposited surface stormdust mass',
3272	& 'kg.m-2.s-1',0,rdsdqdustsurf)
3273	call WRITEDIAGFI(ngrid,'rdsdustq','storm Dust mass mr',
3274	& 'kg/kg',1,rdsqdust)
3275	call WRITEDIAGFI(ngrid,"stormfract",
3276	& "fraction of the mesh,with stormdust",
3277	& "none",0,totstormfract)
3278	call WRITEDIAGFI(ngrid,'rdsqsurf',
3279	& 'stormdust at surface',
3280	& 'kg.m-2',0,qsurf(:,igcm_stormdust_mass))
3281	call WRITEDIAGFI(ngrid,'qsurf',
3282	& 'dust mass at surface',
3283	& 'kg.m-2',0,qsurf(:,igcm_dust_mass))
3284	call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust',
3285	& 'm/s',1,wspeed)
3286	call WRITEDIAGFI(ngrid,'totaldustq','total dust mass',
3287	& 'kg/kg',1,qdusttotal)
3288	call WRITEDIAGFI(ngrid,'dsords',
3289	& 'density scaled opacity of stormdust',
3290	& 'm2.kg-1',1,dsords)
3291	call WRITEDIAGFI(ngrid,'zdqsed_dustq'
3292	& ,'sedimentation q','kg.m-2.s-1',1,
3293	& zdqsed(1,:,igcm_dust_mass))
3294	call WRITEDIAGFI(ngrid,'zdqsed_rdsq'
3295	& ,'sedimentation q','kg.m-2.s-1',1,
3296	& zdqsed(1,:,igcm_stormdust_mass))
3297	endif !(rdstorm 1D)
3298
3299	if (water.AND.tifeedback) then
3300	call WRITEDIAGFI(ngrid,"soiltemp",
3301	& "Soil temperature","K",
3302	& 1,tsoil)
3303	call WRITEDIAGFI(ngrid,'soilti',
3304	& 'Soil Thermal Inertia',
3305	& 'J.s-1/2.m-2.K-1',1,inertiesoil)
3306	endif
3307	end if
3308
3309	cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc
3310	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3311	IF (water) THEN
3312
3313	if (.not.activice) then
3314
3315	tauTES=0
3316	do l=1,nlayer
3317	Qabsice = min(
3318	& max(0.4e6rice(1,l)(1.+nuice_ref)-0.05 ,0.),1.2
3319	& )
3320	opTES(1,l)= 0.75 * Qabsice *
3321	& zq(1,l,igcm_h2o_ice) *
3322	& (zplev(1,l) - zplev(1,l+1)) / g
3323	& / (rho_ice * rice(1,l) * (1.+nuice_ref))
3324	tauTES=tauTES+ opTES(1,l)
3325	enddo
3326	CALL WRITEDIAGFI(ngrid,'tauTESap',
3327	& 'tau abs 825 cm-1',
3328	& '',0,tauTES)
3329	else
3330
3331	CALL WRITEDIAGFI(ngrid,'tauTES',
3332	& 'tau abs 825 cm-1',
3333	& '',0,taucloudtes)
3334	endif
3335
3336	mtot = 0
3337	icetot = 0
3338	h2otot = qsurf(1,igcm_h2o_ice)
3339
3340	do l=1,nlayer
3341	mtot = mtot + zq(1,l,igcm_h2o_vap)
3342	& * (zplev(1,l) - zplev(1,l+1)) / g
3343	icetot = icetot + zq(1,l,igcm_h2o_ice)
3344	& * (zplev(1,l) - zplev(1,l+1)) / g
3345	end do
3346	h2otot = h2otot+mtot+icetot
3347
3348	CALL WRITEDIAGFI(ngrid,'h2otot',
3349	& 'h2otot',
3350	& 'kg/m2',0,h2otot)
3351	CALL WRITEDIAGFI(ngrid,'mtot',
3352	& 'mtot',
3353	& 'kg/m2',0,mtot)
3354	CALL WRITEDIAGFI(ngrid,'icetot',
3355	& 'icetot',
3356	& 'kg/m2',0,icetot)
3357
3358	if (scavenging) then
3359
3360	rave = 0
3361	do l=1,nlayer
3362	cccc Column integrated effective ice radius
3363	cccc is weighted by total ice surface area (BETTER)
3364	rave = rave + tauscaling(1) *
3365	& zq(1,l,igcm_ccn_number) *
3366	& (zplev(1,l) - zplev(1,l+1)) / g *
3367	& rice(1,l) * rice(1,l)* (1.+nuice_ref)
3368	enddo
3369	rave=icetot0.75/max(ravepi*rho_ice,1.e-30) ! surface weight
3370
3371	Nccntot= 0
3372	call watersat(ngrid*nlayer,zt,zplay,zqsat)
3373	do l=1,nlayer
3374	Nccntot = Nccntot +
3375	& zq(1,l,igcm_ccn_number)*tauscaling(1)
3376	& *(zplev(1,l) - zplev(1,l+1)) / g
3377	satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l)
3378	satu(1,l) = (max(satu(1,l),float(1))-1)
3379	! & * zq(1,l,igcm_h2o_vap) *
3380	! & (zplev(1,l) - zplev(1,l+1)) / g
3381	enddo
3382	call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1,
3383	& satu)
3384	CALL WRITEDIAGFI(ngrid,'Nccntot',
3385	& 'Nccntot',
3386	& 'nbr/m2',0,Nccntot)
3387
3388	call WRITEDIAGFI(ngrid,'zdqsed_dustq'
3389	& ,'sedimentation q','kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass))
3390	call WRITEDIAGFI(ngrid,'zdqsed_dustN'
3391	&,'sedimentation N','Nbr.m-2.s-1',1,
3392	& zdqsed(1,:,igcm_dust_number))
3393
3394	else ! of if (scavenging)
3395
3396	cccc Column integrated effective ice radius
3397	cccc is weighted by total ice mass (LESS GOOD)
3398	rave = 0
3399	do l=1,nlayer
3400	rave = rave + zq(1,l,igcm_h2o_ice)
3401	& * (zplev(1,l) - zplev(1,l+1)) / g
3402	& * rice(1,l) * (1.+nuice_ref)
3403	enddo
3404	rave=max(rave/max(icetot,1.e-30),1.e-30) ! mass weight
3405	endif ! of if (scavenging)
3406
3407
3408	CALL WRITEDIAGFI(ngrid,'reffice',
3409	& 'reffice',
3410	& 'm',0,rave)
3411
3412	!Alternative A. Pottier weighting
3413	rave2 = 0.
3414	totrave2 = 0.
3415	do l=1,nlayer
3416	rave2 =rave2+ zq(1,l,igcm_h2o_ice)*rice(1,l)
3417	totrave2 = totrave2 + zq(1,l,igcm_h2o_ice)
3418	end do
3419	rave2=max(rave2/max(totrave2,1.e-30),1.e-30)
3420	CALL WRITEDIAGFI(ngrid,'rmoym',
3421	& 'reffice',
3422	& 'm',0,rave2)
3423
3424	do iq=1,nq
3425	call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s',
3426	& trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq))
3427	end do
3428
3429	call WRITEDIAGFI(ngrid,'zdqcloud_ice','cloud ice',
3430	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))
3431	call WRITEDIAGFI(ngrid,'zdqcloud_vap','cloud vap',
3432	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap))
3433	call WRITEDIAGFI(ngrid,'zdqcloud','cloud ice',
3434	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)
3435	& +zdqcloud(1,:,igcm_h2o_vap))
3436
3437	call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
3438	& rice)
3439
3440	if (CLFvarying) then
3441	call WRITEDIAGFI(ngrid,'totcloudfrac',
3442	& 'Total cloud fraction',
3443	& ' ',0,totcloudfrac)
3444	endif !clfvarying
3445
3446	ENDIF ! of IF (water)
3447
3448	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3449	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3450
3451
3452	zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g
3453
3454	do l=2,nlayer-1
3455	tmean=zt(1,l)
3456	if(zt(1,l).ne.zt(1,l-1))
3457	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
3458	zlocal(l)= zlocal(l-1)
3459	& -log(zplay(1,l)/zplay(1,l-1))rnew(1,l)tmean/g
3460	enddo
3461	zlocal(nlayer)= zlocal(nlayer-1)-
3462	& log(zplay(1,nlayer)/zplay(1,nlayer-1))*
3463	& rnew(1,nlayer)*tmean/g
3464	#endif
3465
3466	END IF ! if(ngrid.ne.1)
3467
3468	! XIOS outputs
3469	#ifdef CPP_XIOS
3470	! Send fields to XIOS: (NB these fields must also be defined as
3471	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
3472	CALL send_xios_field("ls",zls*180./pi)
3473
3474	CALL send_xios_field("ps",ps)
3475	CALL send_xios_field("area",cell_area)
3476
3477	CALL send_xios_field("temperature",zt)
3478	CALL send_xios_field("u",zu)
3479	CALL send_xios_field("v",zv)
3480
3481	if (lastcall.and.is_omp_master) then
3482	write(,) "physiq lastcall: call xios_context_finalize"
3483	call xios_context_finalize
3484	endif
3485	#endif
3486
3487
3488	icount=icount+1
3489
3490
3491	END SUBROUTINE physiq
3492
3493	END MODULE physiq_mod

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

Download in other formats:

Original Format