Context Navigation

physiq_mod.F @ 2214

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

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

Download in other formats:

Original Format