Context Navigation

physiq.F @ 1032

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

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

Download in other formats:

Original Format