Context Navigation

physiq_mod.F @ 2220

Last change on this file since 2220 was 2220, checked in by dbardet, 5 years ago

Update the nonorographic gravity waves drag parametrization (from the r3599 of Earth s model LMDZ6): 1) add east_ and west_gwdstress variables, 2) delete aleas function: random waves are producted by the MOD function, 3) reproductibility concerning the number of procs is now validated, 4) changing name of some variables like RUW, RVW, ZOP, ZOM,..., 5) tendency of winds due to GW drag are now module variables and written in restart files.

File size: 137.2 KB

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

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

Download in other formats:

Original Format