Context Navigation

physiq.F @ 986

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

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

Download in other formats:

Original Format