Context Navigation

physiq_mod.F @ 1921

Last change on this file since 1921 was 1921, checked in by emillour, 7 years ago

Mars GCM:
CO2 code update:

nucleaCO2.F: code cleanup and use co2useh2o flag to handle cases where

h2o ccns have to be tracked and accounted for.

co2cloud.F & improvedco2clouds.F: code cleanup and use co2useh2o flag to

control updates on water variables if necessary.

physiq.F : cleanup on outputs & compute mtotco2 and icetotco2

File size: 115.4 KB

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

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

Download in other formats:

Original Format