Context Navigation

physiq.F @ 708

Last change on this file since 708 was 708, checked in by aslmd, 12 years ago
MESOSCALE: minor modifications: added surface dust diagnostics. added an indication of SVN version in makemeso.
File size: 82.8 KB

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

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

Download in other formats:

Original Format