Context Navigation

physiq_mod.F @ 2156

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

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

Download in other formats:

Original Format