Context Navigation

physiq_mod.F @ 2522

Last change on this file since 2522 was 2514, checked in by romain.vande, 5 years ago

Mars Dynamico : _ update of r2507: correct ztime_fin for dynamico

_ Test for no writting Diagfi for Dynamico
_ Correction for Xios output : Change of calendar : start_date=0 as a diagfi, timestep for output physic=1s, day_length=ndtphys, it is needed to have output in days in the xml to have the correct time

File size: 154.1 KB

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

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

Download in other formats:

Original Format