Context Navigation

physiq_mod.F @ 3001

Last change on this file since 3001 was 3000, checked in by jbclement, 3 years ago
Mars PCM: Some adaptations to make the 1D model run with XIOS and MPICH for use with 1D PEM. JBC
File size: 158.1 KB

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

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

Download in other formats:

Original Format