Context Navigation

physiq.F @ 1038

Last change on this file since 1038 was 1038, checked in by aslmd, 11 years ago

MESOSCALE
LMDZ.MARS

--> Performed the necessary modifications for dynamic tracers

to work with the mesoscale model (new physics).

--> Added precompiling flag MESOSCALE around pressure modifications

done in revision 883. This makes the mesoscale model become crazy.

--> Added an option -e in makemeso to erase a configuration and start over.

NOTE
--> not sure recent versions (rev>1000) are compliant with nesting compilation.
--> use mesoscale model + new physics with caution. still not stabilized.

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

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

Download in other formats:

Original Format