Context Navigation

physiq_mod.F @ 2275

Last change on this file since 2275 was 2265, checked in by emillour, 6 years ago
Mars GCM: Save "dtau", the opacity difference between model and target dust scenario in the restartfi.nc file so that we have 1+1=2 when running with dust injection schemes. EM
File size: 140.7 KB

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

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

Download in other formats:

Original Format