Context Navigation

physiq.F @ 890

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

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

Download in other formats:

Original Format