Context Navigation

physiq.F @ 804

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

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

Download in other formats:

Original Format