Context Navigation

physiq.F @ 1047

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

Mars GCM:

IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
Also updated EOF routines, everything is now in eofdump_mod.F90
Removed unused routine lectfux.F (in dyn3d)

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

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

Download in other formats:

Original Format