Context Navigation

physiq_mod.F @ 2082

Last change on this file since 2082 was 2081, checked in by mvals, 6 years ago

Mars GCM:

follow-up of the last change in physiq_mod.F: put iaer indices (defined in dimradmars_mod.F) for aerosol(:,:,:) in writediagfi calls.

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

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

Download in other formats:

Original Format