Context Navigation

physiq.F @ 756

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

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

Download in other formats:

Original Format