Context Navigation

physiq_mod.F @ 2449

Last change on this file since 2449 was 2447, checked in by cmathe, 4 years ago
CO2 microphysics: correction r_sedim (co2cloud.F); add: co2 cloud radiatively active
File size: 151.9 KB

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

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

Download in other formats:

Original Format