Context Navigation

physiq.F @ 850

Last change on this file since 850 was 835, checked in by jbmadeleine, 13 years ago
Mars GCM: Cleaned physiq.F (some debugging lines were introduced at revision 833, sorry about that!);
File size: 89.8 KB

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

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

Download in other formats:

Original Format