Context Navigation

physiq_mod.F @ 2596

Last change on this file since 2596 was 2594, checked in by emillour, 4 years ago

Mars GCM:
Fixes and improvements in the Non-orographic GW scheme, namely:

increments are not tendencies
missing rho factor in EP flux computation
missing rho at launch altitude
changed inputs, because R and Cp are needed to compute rho and BV

File size: 160.4 KB

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

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

Download in other formats:

Original Format