Context Navigation

physiq_mod.F @ 1621

Last change on this file since 1621 was 1621, checked in by emillour, 9 years ago

Further work on full dynamics/physics separation.

LMDZ.COMMON:

added phy_common/vertical_layers_mod.F90 to store information on vertical grid. This is where routines in the physics should get the information.
The contents of vertical_layers_mod intialized via dynphy_lonlat/inigeomphy_mod.F90.

LMDZ.MARS:

physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
created an "ini_tracer_mod" routine in module "tracer_mod" for a cleaner initialization of the later.
removed some purely dynamics-related outputs (etot0, zoom parameters, etc.) from diagfi.nc and stats.nc outputs as these informations are not available in the physics.

LMDZ.GENERIC:

physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
added nqtot to tracer_h.F90.
removed some purely dynamics-related outputs (etot0, zoom parameters, etc.) from diagfi.nc and stats.nc outputs as these informations are not available in the physics.

LMDZ.VENUS:

physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
added infotrac_phy.F90 to store information on tracers in the physics. Initialized via iniphysiq.
added cpdet_phy_mod.F90 to store t2tpot etc. functions to be used in the physics. Initialized via iniphysiq. IMPORTANT: there are some hard-coded constants! These should match what is in cpdet_mod.F90 in the dynamics.
got rid of references to moyzon_mod module within the physics. The required variables (tmoy, plevmoy) are passed to the physics as arguments to physiq.

LMDZ.TITAN:

added infotrac_phy.F90 to store information on tracers in the physics. Initialized via iniphysiq.
added cpdet_phy_mod.F90 to store t2tpot etc. functions to be used in the physics.
Extra work required to completely decouple physics and dynamics: moyzon_mod should be cleaned up and information passed from dynamics to physics as as arguments. Likewise moyzon_ch and moyzon_mu should not be queried from logic_mod (which is in the dynamics).

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

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

Download in other formats:

Original Format