Context Navigation

physiq.F @ 1124

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

Mars GCM:

Improved 15um cooling routines, with also better handling of errors.

FGG

File size: 93.1 KB

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

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

Download in other formats:

Original Format