Context Navigation

physiq_mod.F @ 2224

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

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

Download in other formats:

Original Format