Context Navigation

physiq_mod.F @ 2530

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

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

Download in other formats:

Original Format