Context Navigation

physiq.F @ 1009

Last change on this file since 1009 was 1005, checked in by emillour, 11 years ago

Mars GCM:

Bug fix: when running with photochemistry, ccns did not sediment! Fixed initracer.F. Also added that callsedim/newsedim use updated temperatures.
Converted run0 and run_mcd scripts to bash.

File size: 91.0 KB

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

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

Download in other formats:

Original Format