Context Navigation

physiq.F @ 1312

Last change on this file since 1312 was 1312, checked in by emillour, 10 years ago
Generic and Mars GCM: LMDZ.COMMON dynamics sends mass flux to physics and not vertical velocity. Add the computation of vertical velocity from input mass flux in the physics, and also modify "old" LMDZ.GENERIC and LMDZ.MARS dynamics to be consistent. EM
File size: 99.8 KB

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

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

Download in other formats:

Original Format