Context Navigation

physiq_mod.F @ 1930

Last change on this file since 1930 was 1922, checked in by emillour, 8 years ago

Mars GCM:
CO2 code updates:

make co2cloud a module and save mem_* variables (initialized via phys_state_var_init)
make improvedCO2cloud a module
read/write mem_* variables in phyetat0.F and phyredem.F

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

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

Download in other formats:

Original Format