Context Navigation

physiq_mod.F @ 1660

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

Mars GCM:

Added possibility to run with an Helium "he" tracer (to be initialized at constant value of 3.6e-7 kg/kg_air, i.e. the 4ppm of Krasnopolsky 1996 EUVE satellite, using newstart).
corrected case for CH4 in aeronomars/photochemistry.F (was using undefined indexes when there is no CH4).
updated aki/cpi coefficients for Argon used to compute mean atmospheric Cp and thermal conductivity in aeronomars/concentrations.F

JYC+EM

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

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

Download in other formats:

Original Format