Context Navigation

physiq.F @ 758

Last change on this file since 758 was 758, checked in by emillour, 12 years ago

Mars GCM:

Some minor changes and adaptations for the MCDv5 runs:

removed the 'nice' instruction from run0 script
update dissipation factors in inidissip.F (fac_mid=3,fac_up=30,startalt=70)
add some outputs in physiq.F
put example def files used for MCDv5 simulations in deftank: callphys.def.MCD5 run.def.64x48x49.MCD5 traceur.def.MCD5

File size: 87.8 KB

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

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

Download in other formats:

Original Format