Context Navigation

physiq.F @ 3400

Last change on this file since 3400 was 3373, checked in by afalco, 20 months ago
Pluto.old PCM: Corrected some runtime issues with gfortran. AF
File size: 82.7 KB

Line
1	subroutine physiq(ngrid,nlayer,nq,
2	* firstcall,lastcall,
3	* pday,ptime,ptimestep,
4	* pplev,pplay,pphi,
5	* pu,pv,pt,pq,
6	* pw,
7	* pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
8
9	use radinc_h, only : naerkind
10	use planet_h
11	USE ioipsl_getincom
12	use comgeomfi_h
13	use comsoil_h
14	use aerosol_mod
15	implicit none
16
17	!==================================================================
18	! Purpose
19	! -------
20	! Central subroutine for all the physics parameterisations in the
21	! universal model. Originally adapted from the Mars LMDZ model.
22	!
23	! The model can be run without or with tracer transport
24	! depending on the value of "tracer" in file "callphys.def".
25	!
26	! It includes:
27	!
28	! 1. Initialization:
29	! 1.1 Firstcall initializations
30	! 1.2 Initialization for every call to physiq
31	! 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
32	! 2. Compute radiative transfer tendencies
33	! (longwave and shortwave).
34	! 4. Vertical diffusion (turbulent mixing):
35	! 5. Convective adjustment
36	! 6. Condensation and sublimation of gases (currently just CO2).
37	! 7. TRACERS :
38	! 7c. other schemes for tracer transport (lifting, sedimentation)
39	! 7d. updates (pressure variations, surface budget)
40	! 9. Surface and sub-surface temperature calculations
41	! 10. Write outputs :
42	! - "startfi", "histfi" if it's time
43	! - Saving statistics if "callstats = .true."
44	! - Output any needed variables in "diagfi"
45	! 10. Diagnostics: mass conservation of tracers, radiative energy balance etc.
46	!
47	! arguments
48	! ---------
49	!
50	! input
51	! -----
52	! ecri period (in dynamical timestep) to write output
53	! ngrid Size of the horizontal grid.
54	! All internal loops are performed on that grid.
55	! nlayer Number of vertical layers.
56	! nq Number of advected fields
57	! firstcall True at the first call
58	! lastcall True at the last call
59	! pday Number of days counted from the North. Spring
60	! equinoxe.
61	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
62	! ptimestep timestep (s)
63	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
64	! pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
65	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
66	! pu(ngrid,nlayer) u component of the wind (ms-1)
67	! pv(ngrid,nlayer) v component of the wind (ms-1)
68	! pt(ngrid,nlayer) Temperature (K)
69	! pq(ngrid,nlayer,nq) Advected fields
70	! pudyn(ngrid,nlayer) \
71	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
72	! ptdyn(ngrid,nlayer) / corresponding variables
73	! pqdyn(ngrid,nlayer,nq) /
74	! pw(ngrid,?) vertical velocity
75	!
76	! output
77	! ------
78	!
79	! pdu(ngrid,nlayermx) \
80	! pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding
81	! pdt(ngrid,nlayermx) / variables due to physical processes.
82	! pdq(ngrid,nlayermx) /
83	! pdpsrf(ngrid) /
84	! tracerdyn call tracer in dynamical part of GCM ?
85	!
86	! Authors
87	! -------
88	! Frederic Hourdin 15/10/93
89	! Francois Forget 1994
90	! Christophe Hourdin 02/1997
91	! Subroutine completly rewritten by F.Forget (01/2000)
92	! Water ice clouds: Franck Montmessin (update 06/2003)
93	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
94	! New correlated-k radiative scheme: R. Wordsworth (2009)
95	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
96	! Pluto version improved :T. Bertrand (2015,2020)
97	!==================================================================
98
99	c 0. Declarations :
100	c ------------------
101	#include "dimensions.h"
102	#include "dimphys.h"
103	#include "surfdat.h"
104	#include "comdiurn.h"
105	#include "callkeys.h"
106	#include "comcstfi.h"
107	#include "comsaison.h"
108	#include "control.h"
109	#include "comg1d.h"
110	#include "tracer.h"
111	#include "netcdf.inc"
112
113	c Arguments :
114	c -----------
115
116	c inputs:
117	c -------
118	INTEGER ngrid,nlayer,nq
119	REAL ptimestep
120	REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
121	REAL pphi(ngridmx,nlayer)
122	REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
123	REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
124	REAL pw(ngridmx,nlayer) ! pvervel transmitted by dyn3d
125	REAL zh(ngridmx,nlayermx) ! potential temperature (K)
126	LOGICAL firstcall,lastcall
127	REAL pday
128	REAL ptime
129	logical tracerdyn
130
131	c outputs:
132	c --------
133	c physical tendencies
134	REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
135	REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
136	REAL pdpsrf(ngridmx) ! surface pressure tendency
137	REAL picen2(ngrid)
138
139	c Local saved variables:
140	c ----------------------
141	INTEGER day_ini ! Initial date of the run (sol since Ls=0)
142	INTEGER icount ! counter of calls to physiq during the run.
143	REAL tsurf(ngridmx) ! Surface temperature (K)
144	REAL tsoil(ngridmx,nsoilmx) ! sub-surface temperatures (K)
145	REAL,SAVE :: tidat(ngridmx,nsoilmx) ! thermal inertia soil
146	REAL tidat_out(ngridmx,nlayermx) ! thermal inertia soil but output on vertical grid
147	REAL tsub(ngridmx) ! temperature of the deepest layer
148	REAL albedo(ngridmx) ! Surface albedo
149
150	REAL emis(ngridmx) ! Thermal IR surface emissivity
151	REAL,SAVE :: dtrad(ngridmx,nlayermx) ! Net atm. radiative heating rate (K.s-1)
152	REAL,SAVE :: fluxrad_sky(ngridmx) ! rad. flux from sky absorbed by surface (W.m-2)
153	REAL,SAVE :: fluxrad(ngridmx) ! Net radiative surface flux (W.m-2)
154	REAL,SAVE :: capcal(ngridmx) ! surface heat capacity (J m-2 K-1)
155	REAL dplanck(ngridmx) ! for output (W.s/m2/K)
156	REAL,SAVE :: fluxgrd(ngridmx) ! surface conduction flux (W.m-2)
157	REAL,SAVE :: qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
158	REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy
159	REAL qsurf1(ngridmx,nqmx) ! saving qsurf to calculate flux over long timescales kg.m-2
160	REAL flusurf(ngridmx,nqmx) ! flux cond/sub kg.m-2.s-1
161	REAL,SAVE :: flusurfold(ngridmx,nqmx) ! old flux cond/sub kg.m-2.s-1
162	REAL totmass ! total mass of n2 (atm + surf)
163	REAL globave ! globalaverage 2D ps
164	REAL globaveice(nqmx) ! globalaverage 2D ice
165	REAL globavenewice(nqmx) ! globalaverage 2D ice
166
167	REAL,SAVE :: ptime0 ! store the first time
168
169	REAL dstep
170	REAL,SAVE :: glastep=20 ! step en pluto day pour etaler le glacier
171
172	SAVE day_ini,icount
173	SAVE tsurf,tsoil,tsub
174	SAVE albedo,emis,q2
175
176
177	REAL stephan
178	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
179	SAVE stephan
180
181	REAL acond,bcond
182	SAVE acond,bcond
183	REAL tcond1p4Pa
184	DATA tcond1p4Pa/38/
185
186	c Local variables :
187	c -----------------
188	! Tendencies for the paleoclimate mode
189	REAL,SAVE :: qsurfyear(ngridmx,nqmx) ! kg.m-2 averaged mass of ice lost/gained in the last Pluto year of the run
190	REAL,SAVE :: phisfinew(ngridmx) ! geopotential of the bedrock (= phisfi-qsurf/1000*g)
191	REAL qsurfpal(ngridmx,nqmx) ! qsurf after a paleoclimate step : for physdem1 and restartfi
192	REAL phisfipal(ngridmx) ! geopotential after a paleoclimate step : for physdem1 and restartfi
193	REAL oblipal ! change of obliquity
194	REAL peri_daypal ! new periday
195	REAL eccpal ! change of eccentricity
196	REAL tpalnew ! change of time
197	REAL adjustnew ! change in N2 ice albedo
198	REAL pdaypal ! new pday = day_ini + step
199	REAL zdt_tot ! time range corresponding to the flux of qsurfyear
200	REAL massacc(nqmx) ! accumulated mass
201	REAL masslost(nqmx) ! accumulated mass
202
203	! aerosol (ice or haze) extinction optical depth at reference wavelength
204	! for the "naerkind" optically active aerosols:
205	REAL aerosol(ngridmx,nlayermx,naerkind)
206
207	CHARACTER*80 fichier
208	INTEGER l,ig,ierr,igout,iq,i, tapphys
209	INTEGER lecttsoil ! lecture of tsoil from proftsoil
210	! INTEGER iqmin ! Used if iceparty engaged
211
212	REAL fluxsurf_lw(ngridmx) ! incident LW (IR) surface flux (W.m-2)
213	REAL fluxsurf_sw(ngridmx) ! incident SW (solar) surface flux (W.m-2)
214	REAL,SAVE :: fluxtop_lw(ngridmx) ! Outgoing LW (IR) flux to space (W.m-2)
215	REAL,SAVE :: fluxtop_sw(ngridmx) ! Outgoing SW (solar) flux to space (W.m-2)
216
217	! included by RW for equilibration test
218	real,save :: fluxtop_dn(ngridmx)
219	real fluxabs_sw(ngridmx) ! absorbed shortwave radiation
220	real fluxdyn(ngridmx) ! horizontal heat transport by dynamics
221
222	REAL zls ! solar longitude (rad)
223	REAL zfluxuv ! Lyman alpha flux at 1AU
224	! ph/cm2/s
225	REAL zday ! date (time since Ls=0, in martian days)
226	REAL saveday ! saved date
227	SAVE saveday
228	REAL savedeclin ! saved declin
229	SAVE savedeclin
230	!REAL adjust ! adjustment for surface albedo
231	REAL zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
232	REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
233
234	c Aerosol effective radius used for radiative transfer (units=meters)
235	REAL reffrad(ngridmx,nlayermx,naerkind)
236
237	c Tendencies due to various processes:
238	REAL dqsurf(ngridmx,nqmx)
239	REAL zdtlw(ngridmx,nlayermx) ! (K/s)
240	REAL zdtsw(ngridmx,nlayermx) ! (K/s)
241	REAL cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear areas
242	REAL cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear areas
243	REAL zdtsurf(ngridmx) ! (K/s)
244	REAL zdtlscale(ngridmx,nlayermx)
245	REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2)
246	REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s)
247	REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2)
248	REAL zdhadj(ngridmx,nlayermx) ! (K/s)
249	REAL zdtgw(ngridmx,nlayermx) ! (K/s)
250	REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2)
251	REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
252	REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
253	REAL tconds(ngridmx,nlayermx)
254
255	REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
256	REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
257	REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
258	REAL zdqadj(ngridmx,nlayermx,nqmx)
259	REAL zdqc(ngridmx,nlayermx,nqmx)
260	REAL zdqlscale(ngridmx,nlayermx,nqmx)
261	REAL zdqslscale(ngridmx,nqmx), zdqsc(ngridmx,nqmx)
262	REAL zdqchim(ngridmx,nlayermx,nqmx)
263	REAL zdqschim(ngridmx,nqmx)
264	REAL zdqch4cloud(ngridmx,nlayermx,nqmx)
265	REAL zdqsch4cloud(ngridmx,nqmx)
266	REAL zdtch4cloud(ngridmx,nlayermx)
267	REAL zdqcocloud(ngridmx,nlayermx,nqmx)
268	REAL zdqscocloud(ngridmx,nqmx)
269	REAL zdtcocloud(ngridmx,nlayermx)
270	REAL rice_ch4(ngridmx,nlayermx) ! Methane ice geometric mean radius (m)
271	REAL rice_co(ngridmx,nlayermx) ! CO ice geometric mean radius (m)
272
273	REAL zdqsch4fast(ngridmx) ! used only for fast model nogcm
274	REAL zdqch4fast(ngridmx) ! used only for fast model nogcm
275	REAL zdqscofast(ngridmx) ! used only for fast model nogcm
276	REAL zdqcofast(ngridmx) ! used only for fast model nogcm
277	REAL zdqflow(ngridmx,nqmx)
278
279	REAL zdteuv(ngridmx,nlayermx) ! (K/s)
280	REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
281	REAL zdumolvis(ngridmx,nlayermx)
282	REAL zdvmolvis(ngridmx,nlayermx)
283	real zdqmoldiff(ngridmx,nlayermx,nqmx)
284
285	! Haze relatated tendancies
286	REAL zdqhaze(ngridmx,nlayermx,nqmx)
287	REAL zdqprodhaze(ngridmx,nqmx)
288	REAL zdqsprodhaze(ngridmx)
289	REAL zdqhaze_col(ngridmx)
290	REAL zdqphot_prec(ngrid,nlayer)
291	REAL zdqphot_ch4(ngrid,nlayer)
292	REAL zdqconv_prec(ngrid,nlayer)
293	REAL zdq_source(ngridmx,nlayermx,nqmx)
294	! Fast Haze relatated tendancies
295	REAL fluxbot(ngridmx)
296	REAL gradflux(ngridmx)
297	REAL fluxlym_sol_bot(ngridmx) ! Solar flux Lyman alpha ph.m-2.s-1 reaching the surface
298	REAL fluxlym_ipm_bot(ngridmx) ! IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 reaching the surface
299	REAL flym_sol(ngridmx) ! Incident Solar flux Lyman alpha ph.m-2.s-1
300	REAL flym_ipm(ngridmx) ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1
301
302	real nconsMAX
303	real vdifcncons(ngridmx)
304	real cadjncons(ngridmx)
305	real ch4csMAX
306	real ch4cncons(ngridmx)
307	real ch4sedncons(ngridmx)
308	real ch4sedsMAX
309	real condncons(ngridmx)
310	real dWtot, dWtots
311	real dWtotn2, dWtotsn2
312	real nconsMAXn2
313	real vdifcnconsn2(ngridmx)
314	real cadjnconsn2(ngridmx)
315	real condnconsn2(ngridmx)
316
317	c Local variable for local intermediate calcul:
318	REAL zflubid(ngridmx)
319	REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
320	REAL zdum1(ngridmx,nlayermx)
321	REAL zdum2(ngridmx,nlayermx)
322	REAL ztim1,ztim2,ztim3, z1,z2
323	REAL ztime_fin
324	REAL zdh(ngridmx,nlayermx)
325	INTEGER length
326	PARAMETER (length=100)
327
328	c local variables only used for diagnostic (output in file "diagfi" or "stats")
329	c -----------------------------------------------------------------------------
330	REAL ps(ngridmx), zt(ngridmx,nlayermx)
331	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
332	REAL zq(ngridmx,nlayermx,nqmx)
333	REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
334	character*2 str2
335	character*5 str5
336	real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
337	real zdtdyn(ngridmx,nlayermx),ztprevious(ngridmx,nlayermx)
338	save ztprevious
339	real qtot1,qtot2 ! total aerosol mass
340	integer igmin, lmin
341	logical tdiag
342
343	REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
344	REAL cd
345	real vmr(ngridmx,nlayermx) ! volume mixing ratio
346
347	REAL time_phys
348
349	REAL tau_col(ngrid)
350	real flux1,flux2,flux3,ts1,ts2,ts3
351
352	real qcol(ngridmx,nqmx)
353	! Pluto adding variables
354	real vmr_ch4(ngridmx) ! vmr ch4
355	real vmr_co(ngridmx) ! vmr co
356	real rho(ngridmx,nlayermx) ! density
357	real zrho_ch4(ngridmx,nlayermx) ! density methane kg.m-3
358	real zrho_co(ngridmx,nlayermx) ! density CO kg.m-3
359	real zrho_haze(ngridmx,nlayermx) ! density haze kg.m-3
360	real zdqrho_photprec(ngridmx,nlayermx) !photolysis rate kg.m-3.s-1
361	real zq1temp_ch4(ngridmx) !
362	real qsat_ch4(ngridmx) !
363	real qsat_ch4_l1(ngridmx) !
364
365	! CHARACTER(LEN=20) :: txt ! to temporarly store text for eddy tracers
366	real profmmr(ngridmx,nlayermx) ! fixed profile of haze if haze_proffix
367	real sensiblehf1(ngridmx) ! sensible heat flux
368	real sensiblehf2(ngridmx) ! sensible heat flux
369
370	! included by RW to test energy conservation
371	REAL dEtot, dEtots, masse, vabs, dvabs, ang0
372
373	! included by RW to allow variations in cp with location
374	REAL cpp3D(ngridmx,nlayermx) ! specific heat capacity at const. pressure
375	REAL rcp3D(ngridmx,nlayermx) ! R / specific heat capacity
376	real cppNI, rcpNI, nnu ! last one just for Seb version
377	REAL zpopskNI(ngridmx,nlayermx)
378
379	! included by RW to make 1D saves not every timestep
380	integer countG1D
381	! integer countG3D,saveG3D
382	save countG1D
383	! save countG3D,saveG3D
384	CHARACTER(len=10) :: tname
385
386	c=======================================================================
387
388	c 1. Initialisation:
389	c -----------------
390
391	c 1.1 Initialisation only at first call
392	c ---------------------------------------
393	IF (firstcall) THEN
394	if(ngrid.eq.1)then
395	saveG1D=1
396	countG1D=1
397	endif
398
399	c variables set to 0
400	c ~~~~~~~~~~~~~~~~~~
401	call zerophys(ngrid*nlayer,dtrad)
402	call zerophys(ngrid,fluxrad)
403	call zerophys(ngridnlayernq,zdqsed)
404	call zerophys(ngrid*nq, zdqssed)
405	call zerophys(ngridnlayernq,zdqadj)
406
407	c initialize tracer names, indexes and properties
408	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
409	tracerdyn=tracer ! variable tracer in dynamics
410	IF (tracer) THEN
411	CALL initracer()
412	ENDIF ! end tracer
413
414	! Initialize aerosol indexes: done in initracer
415	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
416	! call iniaerosol()
417
418	c read startfi (initial parameters)
419	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
420	CALL phyetat0 ("startfi.nc",0,0,
421	& nsoilmx,nq,
422	& day_ini,time_phys,
423	& tsurf,tsoil,emis,q2,qsurf)
424
425	if (pday.ne.day_ini) then
426	write(,) "PHYSIQ: ERROR: bad synchronization between ",
427	& "physics and dynamics"
428	write(,) "dynamics day: ",pday
429	write(,) "physics day: ",day_ini
430	stop
431	endif
432
433	write (,) 'In physiq day_ini =', day_ini
434	ptime0=ptime
435	write (,) 'In physiq ptime0 =', ptime
436
437	c Initialize albedo and orbital calculation
438	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
439	CALL surfini(ngrid,albedo)
440	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
441
442	saveday=pday
443	savedeclin=0.
444	!adjust=0. ! albedo adjustment for convergeps
445
446	c Initialize surface inventory and geopotential for paleo mode
447	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
448	if (paleo) then
449	write(,) 'Paleo time tpal = ',tpal
450	DO ig=1,ngrid
451	phisfinew(ig)=phisfi(ig)-qsurf(ig,1)*g/1000. ! topo of bedrock below the ice
452	ENDDO
453	endif
454
455	if (conservn2) then
456	write(,) 'conservn2 initial'
457	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
458	endif
459
460	c initialize soil
461	c ~~~~~~~~~~~~~~~
462	IF (callsoil) THEN
463	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
464	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
465	ELSE
466	PRINT*,'WARNING! Thermal conduction in the soil turned off'
467	DO ig=1,ngrid
468	capcal(ig)=1.e5
469	fluxgrd(ig)=0.
470	ENDDO
471	ENDIF
472	icount=1
473
474	ENDIF ! (end of "if firstcall")
475
476	c ---------------------------------------------------
477	c 1.2 Initializations done at every physical timestep:
478	c ---------------------------------------------------
479	c
480	IF (ngrid.NE.ngridmx) THEN
481	PRINT*,'STOP in PHYSIQ'
482	PRINT*,'Probleme de dimensions :'
483	PRINT*,'ngrid = ',ngrid
484	PRINT*,'ngridmx = ',ngridmx
485	STOP
486	ENDIF
487
488	c Initialize various variables
489	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
490	call zerophys(ngrid*nlayer, pdv)
491	call zerophys(ngrid*nlayer, pdu)
492	call zerophys(ngrid*nlayer,pdt)
493	call zerophys(ngridnlayernq, pdq)
494	call zerophys(ngrid, pdpsrf)
495	call zerophys(ngrid, zflubid)
496	call zerophys(ngrid, zdtsurf)
497	zdtsdif(1:ngrid)=0
498	call zerophys(ngrid*nq, dqsurf)
499
500	if (conservn2) then
501	write(,) 'conservn2 iniloop'
502	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
503	endif
504
505	igout=ngrid/2+1
506
507	zday=pday+ptime ! compute time, in sols (and fraction thereof)
508
509	c Compute Solar Longitude (Ls) :
510	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
511	if (season) then
512	call solarlong(zday,zls)
513	else
514	call solarlong(float(day_ini),zls)
515	end if
516
517	c Get Lyman alpha flux at specific Ls
518	if (haze) then
519	call lymalpha(zls,zfluxuv)
520	print*, 'Haze lyman-alpha zls,zfluxuv=',zls,zfluxuv
521	end if
522
523	c Compute geopotential between layers
524	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
525	DO l=1,nlayer
526	DO ig=1,ngrid
527	zzlay(ig,l)=pphi(ig,l)/g
528	ENDDO
529	ENDDO
530	DO ig=1,ngrid
531	zzlev(ig,1)=0.
532	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
533	ENDDO
534	DO l=2,nlayer
535	DO ig=1,ngrid
536	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
537	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
538	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
539	ENDDO
540	ENDDO
541
542	! Potential temperature calculation not the same in physiq and dynamic...
543	c Compute potential temperature zh
544	DO l=1,nlayer
545	DO ig=1,ngrid
546	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
547	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
548	ENDDO
549	ENDDO
550
551	c --------------------------------------------------------
552	c 1.3 thermosphere
553	c --------------------------------------------------------
554
555	c ajout de la conduction depuis la thermosphere
556	IF (callconduct) THEN
557
558	call conduction (ngrid,nlayer,ptimestep,
559	$ pplay,pt,zzlay,zzlev,zdtconduc,tsurf)
560	DO l=1,nlayer
561	DO ig=1,ngrid
562	pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
563	ENDDO
564	ENDDO
565
566	ENDIF
567
568	! ajout de la viscosite moleculaire
569	IF (callmolvis) THEN
570	call molvis(ngrid,nlayer,ptimestep,
571	$ pplay,pt,zzlay,zzlev,
572	$ zdtconduc,pu,tsurf,zdumolvis)
573	call molvis(ngrid,nlayer,ptimestep,
574	$ pplay,pt,zzlay,zzlev,
575	$ zdtconduc,pv,tsurf,zdvmolvis)
576
577	DO l=1,nlayer
578	DO ig=1,ngrid
579	! pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
580	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
581	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
582	ENDDO
583	ENDDO
584	ENDIF
585
586	IF (callmoldiff) THEN
587	call moldiff_red(ngrid,nlayer,nq,
588	& pplay,pplev,pt,pdt,pq,pdq,ptimestep,
589	& zzlay,zdtconduc,zdqmoldiff)
590
591	DO l=1,nlayer
592	DO ig=1,ngrid
593	DO iq=1, nq
594	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
595	ENDDO
596	ENDDO
597	ENDDO
598	ENDIF
599
600	if (conservn2) then
601	write(,) 'conservn2 thermo'
602	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
603	endif
604
605	c-----------------------------------------------------------------------
606	c 2. Compute radiative tendencies :
607	c--------- --------------------------------------------------------------
608	c Saving qsurf to compute paleo flux condensation/sublimation
609	DO iq=1, nq
610	DO ig=1,ngrid
611	IF (qsurf(ig,iq).lt.0.) then
612	qsurf(ig,iq)=0.
613	ENDIF
614	qsurf1(ig,iq)=qsurf(ig,iq)
615	ENDDO
616	ENDDO
617
618	IF (callrad) THEN
619	IF( MOD(icount-1,iradia).EQ.0) THEN
620
621	c Local stellar zenith angle
622	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
623	IF (triton) then
624	CALL orbitetriton(zls,zday,dist_sol,declin)
625	ELSE
626	CALL orbite(zls,dist_sol,declin)
627	ENDIF
628
629	IF (diurnal) THEN
630	ztim1=SIN(declin)
631	ztim2=COS(declin)COS(2.pi*(zday-.5))
632	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
633
634	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
635	s ztim1,ztim2,ztim3,mu0,fract)
636
637	ELSE
638	CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
639	! WARNING: this function appears not to work in 1D
640	ENDIF
641
642	c Albedo /IT changes depending of surface ices (only in 3D)
643	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
644	if (ngrid.ne.1) then
645
646	!! Specific to change albedo of N2 so that Psurf
647	!! converges toward 1.4 Pa in "1989" seasons for Triton
648	!! converges toward 1.1 Pa in "2015" seasons for Pluto
649	if (convergeps) then
650	if (triton) then
651	! 1989 declination
652	if (declin180./pi.gt.-46..and.declin180./pi.lt.-45.
653	& .and.zday.gt.saveday+1000.
654	& .and.declin.lt.savedeclin) then
655	call globalaverage2d(ngrid,pplev(:,1),globave)
656	if (globave.gt.1.5) then
657	adjust=adjust+0.005
658	else if (globave.lt.1.3) then
659	adjust=adjust-0.005
660	endif
661	saveday=zday
662	endif
663	call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,capcal,
664	& adjust,dist_sol,albedo,emis,flusurfold,ptimestep,zls)
665	savedeclin=declin
666	else
667	! Pluto : 2015 declination current epoch
668	if (declin180./pi.gt.50.and.declin180./pi.lt.51.
669	& .and.zday.gt.saveday+10000.
670	& .and.declin.gt.savedeclin) then
671	call globalaverage2d(ngrid,pplev(:,1),globave)
672	if (globave.gt.1.2) then
673	adjust=adjust+0.005
674	else if (globave.lt.1.) then
675	adjust=adjust-0.005
676	endif
677	saveday=zday
678	endif
679	call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,capcal,
680	& adjust,dist_sol,albedo,emis,flusurfold,ptimestep,zls)
681	savedeclin=declin
682	endif
683	else
684	! Surface properties
685	call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,
686	& capcal,adjust,dist_sol,albedo,emis,flusurfold,ptimestep,
687	& zls)
688	end if
689	!TB22
690	else ! here ngrid=1
691	call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,
692	& capcal,adjust,dist_sol,albedo,emis,flusurfold,ptimestep,
693	& zls)
694
695	end if ! if ngrid ne 1
696
697	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
698	c Fixed zenith angle in 1D
699	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
700	if(ngrid.eq.1.and.globmean1d) then
701	!global mean 1D radiative equiilium
702	ig=1
703	mu0(ig) = 0.5
704	fract(ig)= 1/(4*mu0(ig))
705	!print*,'WARNING GLOBAL MEAN CALCULATION'
706	endif
707
708	c Call main radiative transfer scheme
709	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
710
711	c Radiative transfer
712	c ------------------
713	if (corrk) then
714	c print*,'TB22 start corrk',tsurf ! Bug testphys1d nodebug
715	CALL callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,
716	& albedo,emis,mu0,pplev,pplay,pt,
717	& tsurf,fract,dist_sol,igout,aerosol,cpp3D,
718	& zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,
719	& fluxtop_sw,fluxtop_dn,reffrad,tau_col,ptime,pday,
720	& firstcall,lastcall,zzlay)
721
722	c Radiative flux from the sky absorbed by the surface (W.m-2)
723	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
724	DO ig=1,ngrid
725	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
726	& +fluxsurf_sw(ig)*(1.-albedo(ig))
727	ENDDO
728
729	c Net atmospheric radiative heating/cooling rate (K.s-1):
730	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
731	DO l=1,nlayer
732	DO ig=1,ngrid
733	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
734	ENDDO
735	ENDDO
736
737	else ! corrk = F
738
739	c Atmosphere has no radiative effect
740	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
741	DO ig=1,ngrid
742	fluxtop_dn(ig) = fract(ig)mu0(ig)Fat1AU/dist_sol**2
743	fluxrad_sky(ig) = fluxtop_dn(ig)*(1.-albedo(ig))
744	fluxtop_sw(ig) = fluxtop_dn(ig)*albedo(ig)
745	fluxtop_lw(ig) = emis(ig)stephantsurf(ig)**4
746	ENDDO ! radiation skips the atmosphere entirely
747	DO l=1,nlayer
748	DO ig=1,ngrid
749	dtrad(ig,l)=0.0
750	ENDDO
751	ENDDO ! hence no atmospheric radiative heating
752
753	endif ! if corrk
754
755	ENDIF ! of if(mod(icount-1,iradia).eq.0) ie radiative timestep
756
757	! Transformation of the radiative tendencies:
758	! -------------------------------------------
759
760	! Net radiative surface flux (W.m-2)
761	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
762
763	DO ig=1,ngrid
764	zplanck(ig)=tsurf(ig)*tsurf(ig)
765	zplanck(ig)=emis(ig)*
766	& stephanzplanck(ig)zplanck(ig)
767	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
768	ENDDO
769
770	DO l=1,nlayer
771	DO ig=1,ngrid
772	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
773	ENDDO
774	ENDDO
775
776	ENDIF ! of IF (callrad)
777
778	if (conservn2) then
779	write(,) 'conservn2 radiat'
780	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
781	endif
782
783	!-----------------------------------------------------------------------
784	! 4. Vertical diffusion (turbulent mixing):
785	!-----------------------------------------------------------------------
786
787	IF (calldifv) THEN
788
789	DO ig=1,ngrid
790	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
791	ENDDO
792
793	CALL zerophys(ngrid*nlayer,zdum1)
794	CALL zerophys(ngrid*nlayer,zdum2)
795	do l=1,nlayer
796	do ig=1,ngrid
797	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
798	enddo
799	enddo
800
801	c Calling vdif (Martian version WITH N2 condensation)
802	CALL vdifc(ngrid,nlayer,nq,zpopsk,
803	& ptimestep,capcal,lwrite,
804	& pplay,pplev,zzlay,zzlev,z0,
805	& pu,pv,zh,pq,pt,tsurf,emis,qsurf,
806	& zdum1,zdum2,zdh,pdq,pdt,zflubid,
807	& zdudif,zdvdif,zdhdif,zdtsdif,q2,
808	& zdqdif,zdqsdif,qsat_ch4,qsat_ch4_l1) !,zq1temp_ch4,qsat_ch4)
809
810	DO l=1,nlayer
811	DO ig=1,ngrid
812	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
813	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
814	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
815
816	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
817	ENDDO
818	ENDDO
819
820	DO ig=1,ngrid
821	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
822	ENDDO
823
824	bcond=1./tcond1p4Pa
825	acond=r/lw_n2
826
827	if (tracer) then
828	DO iq=1, nq
829	DO l=1,nlayer
830	DO ig=1,ngrid
831	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
832	ENDDO
833	ENDDO
834	ENDDO
835
836	DO iq=1, nq
837	DO ig=1,ngrid
838	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
839	ENDDO
840	ENDDO
841
842	end if ! of if (tracer)
843
844	!------------------------------------------------------------------
845	! test methane conservation
846	c if(methane)then
847	c call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
848	c & ptimestep,pplev,zdqdif,zdqsdif,'CH4',' vdifc ')
849	c endif ! methane
850	!------------------------------------------------------------------
851	! test CO conservation
852	c if(carbox)then
853	c call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
854	c & ptimestep,pplev,zdqdif,zdqsdif,'CO ',' vdifc ')
855	c endif ! carbox
856	!------------------------------------------------------------------
857
858	ELSE ! case with calldifv=F
859	ztim1=4.stephanptimestep
860	DO ig=1,ngrid
861	ztim2=ztim1emis(ig)tsurf(ig)**3
862	z1=capcal(ig)*tsurf(ig)+
863	s ztim2tsurf(ig)+ (fluxrad(ig)+fluxgrd(ig))ptimestep
864	z2= capcal(ig)+ztim2
865	zdtsurf(ig)=(z1/z2 - tsurf(ig))/ptimestep
866
867	! for output:
868	!dplanck(ig)=4.stephanptimestepemis(ig)tsurf(ig)**3
869	ENDDO
870
871	c ------------------------------------------------------------------
872	c Methane surface sublimation and condensation in fast model (nogcm)
873	c ------------------------------------------------------------------
874	if ((methane).and.(fast).and.condmetsurf) THEN
875
876	call ch4surf(ngrid,nlayer,nq,ptimestep,
877	& tsurf,zdtsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf,
878	& zdqch4fast,zdqsch4fast)
879
880	DO ig=1,ngrid
881	dqsurf(ig,igcm_ch4_ice)= dqsurf(ig,igcm_ch4_ice) +
882	& zdqsch4fast(ig)
883	pdq(ig,1,igcm_ch4_gas)= pdq(ig,1,igcm_ch4_gas) +
884	& zdqch4fast(ig)
885	zdtsurf(ig)=zdtsurf(ig)+lw_ch4*zdqsch4fast(ig)/capcal(ig)
886	ENDDO
887	end if
888	c ------------------------------------------------------------------
889	c CO surface sublimation and condensation in fast model (nogcm)
890	c ------------------------------------------------------------------
891	if ((carbox).and.(fast).and.condcosurf) THEN
892
893	call cosurf(ngrid,nlayer,nq,ptimestep,
894	& tsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf,
895	& zdqcofast,zdqscofast)
896
897	DO ig=1,ngrid
898	dqsurf(ig,igcm_co_ice)= dqsurf(ig,igcm_co_ice) +
899	& zdqscofast(ig)
900	pdq(ig,1,igcm_co_gas)= pdq(ig,1,igcm_co_gas) +
901	& zdqcofast(ig)
902	zdtsurf(ig)=zdtsurf(ig)+lw_co*zdqscofast(ig)/capcal(ig)
903	ENDDO
904	end if
905
906	ENDIF ! of IF (calldifv)
907
908	if (conservn2) then
909	write(,) 'conservn2 diff'
910	call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)+
911	& dqsurf(:,1)*ptimestep)
912	endif
913
914	!------------------------------------------------------------------
915	! test CO conservation
916	! if(carbox)then
917	! call testconservfast(ngrid,nlayer,nq,
918	! & ptimestep,pplev,zdqcofast,zdqscofast,'CO ',' vdifc ')
919	! endif ! carbox
920	!------------------------------------------------------------------
921
922	!-----------------------------------------------------------------------
923	! 5. Dry convective adjustment:
924	! -----------------------------
925
926	IF(calladj) THEN
927
928	DO l=1,nlayer
929	DO ig=1,ngrid
930	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
931	ENDDO
932	ENDDO
933	CALL zerophys(ngrid*nlayer,zduadj)
934	CALL zerophys(ngrid*nlayer,zdvadj)
935	CALL zerophys(ngrid*nlayer,zdhadj)
936
937	CALL convadj(ngrid,nlayer,nq,ptimestep,
938	& pplay,pplev,zpopsk,
939	& pu,pv,zh,pq,
940	& pdu,pdv,zdh,pdq,
941	& zduadj,zdvadj,zdhadj,
942	& zdqadj)
943
944
945	DO l=1,nlayer
946	DO ig=1,ngrid
947	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
948	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
949	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
950	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
951	ENDDO
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)+ zdqadj(ig,l,iq)
959	ENDDO
960	ENDDO
961	ENDDO
962	end if
963
964	ENDIF ! of IF(calladj)
965
966	!------------------------------------------------------------------
967	! test methane conservation
968	c if(methane)then
969	c call testchange(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
970	c & ptimestep,pplev,zdqadj,'CH4','calladj')
971	c endif ! methane
972	!------------------------------------------------------------------
973	! test CO conservation
974	c if(carbox)then
975	c call testchange(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
976	c & ptimestep,pplev,zdqadj,'CO ','calladj')
977	c endif ! carbox
978	!------------------------------------------------------------------
979
980	!-----------------------------------------------------------------------
981	! 6. Nitrogen condensation-sublimation:
982	! -------------------------------------------
983
984	IF (N2cond) THEN
985	if (.not.tracer.or.igcm_n2.eq.0) then
986	print*,'We need a n2 ice tracer to condense n2'
987	call abort
988	endif
989
990	c write(,) 'before N2 condens :'
991	c call massn2(ngrid,nlayer,pplev,qsurf(:,1),dqsurf(:,igcm_n2),
992	c & pdpsrf,ptimestep)
993
994	if (tracer) then
995	zdqc(:,:,:)=0.
996	zdqsc(:,:)=0.
997	call condens_n2(ngrid,nlayer,nq,ptimestep,
998	& capcal,pplay,pplev,tsurf,pt,
999	& pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
1000	& qsurf(1,igcm_n2),albedo,emis,
1001	& zdtc,zdtsurfc,pdpsrf,zduc,zdvc,
1002	& zdqc,zdqsc(1,igcm_n2))
1003	endif
1004
1005	!!! TB temporaire
1006	!zdtc=0.
1007	!zdvc=0.
1008	!zduc=0.
1009	!zdqc=0.
1010	!zdqsc=0.
1011	!zdtsurfc=0.
1012
1013	DO l=1,nlayer
1014	DO ig=1,ngrid
1015	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
1016	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
1017	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
1018	ENDDO
1019	ENDDO
1020	DO ig=1,ngrid
1021	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
1022	ENDDO
1023
1024	DO iq=1, nq
1025	c TB: option eddy not condensed with N2
1026	c txt=noms(iq)
1027	c IF (txt(1:4).ne."eddy") THEN
1028	DO l=1,nlayer
1029	DO ig=1,ngrid
1030	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
1031	ENDDO
1032	ENDDO
1033	c ENDIF
1034	ENDDO
1035	DO ig=1,ngrid
1036	dqsurf(ig,igcm_n2)= dqsurf(ig,igcm_n2) + zdqsc(ig,igcm_n2)
1037	ENDDO
1038
1039	ENDIF ! of IF (N2cond)
1040
1041	if (conservn2) then
1042	write(,) 'conservn2 n2cond'
1043	call testconservmass(ngrid,nlayer,pplev(:,1)+
1044	& pdpsrf(:)ptimestep,qsurf(:,1)+dqsurf(:,1)ptimestep)
1045	endif
1046
1047	!------------------------------------------------------------------
1048	! test n2 conservation for one tendency / pplevis not updated here w pdpsrf
1049	! if(tracer)then
1050	! call testconserv(ngrid,nlayer,nq,igcm_n2,igcm_n2,
1051	! & ptimestep,pplev,zdqc,zdqsc,'N2 ',' n2cond')
1052	! call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
1053	! endif ! n2
1054	!------------------------------------------------------------------
1055	! test methane conservation
1056	! if(methane)then
1057	! call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
1058	! & ptimestep,pplev,zdqc,zdqsc,'CH4',' n2cond')
1059	! endif ! methane
1060	!------------------------------------------------------------------
1061	! test CO conservation
1062	c if(carbox)then
1063	c call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
1064	c & ptimestep,pplev,zdqc,zdqsc,'CO ',' n2cond')
1065	c endif ! carbox
1066	!------------------------------------------------------------------
1067
1068	c-----------------------------------------------------------------------
1069	c 7. Specific parameterizations for tracers
1070	c -----------------------------------------
1071
1072	if (tracer) then
1073
1074	c 7a. Methane, CO, and ice
1075	c ---------------------------------------
1076	c Methane ice condensation in the atmosphere
1077	c ----------------------------------------
1078	zdqch4cloud(:,:,:)=0.
1079	if ((methane).and.(metcloud).and.(.not.fast)) THEN
1080	call ch4cloud(ngrid,nlayer,naerkind,ptimestep,
1081	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
1082	& pq,pdq,zdqch4cloud,zdqsch4cloud,zdtch4cloud,
1083	& nq,rice_ch4)
1084
1085	DO l=1,nlayer
1086	DO ig=1,ngrid
1087	pdq(ig,l,igcm_ch4_gas)=pdq(ig,l,igcm_ch4_gas)+
1088	& zdqch4cloud(ig,l,igcm_ch4_gas)
1089	pdq(ig,l,igcm_ch4_ice)=pdq(ig,l,igcm_ch4_ice)+
1090	& zdqch4cloud(ig,l,igcm_ch4_ice)
1091	ENDDO
1092	ENDDO
1093
1094	! Increment methane ice surface tracer tendency
1095	DO ig=1,ngrid
1096	dqsurf(ig,igcm_ch4_ice)=dqsurf(ig,igcm_ch4_ice)+
1097	& zdqsch4cloud(ig,igcm_ch4_ice)
1098	ENDDO
1099
1100	! update temperature tendancy
1101	DO ig=1,ngrid
1102	DO l=1,nlayer
1103	pdt(ig,l)=pdt(ig,l)+zdtch4cloud(ig,l)
1104	ENDDO
1105	ENDDO
1106	else
1107	rice_ch4(:,:)=0 ! initialization needed for callsedim
1108	end if
1109
1110	c ---------------------------------------
1111	c CO ice condensation in the atmosphere
1112	c ----------------------------------------
1113	zdqcocloud(:,:,:)=0.
1114	IF ((carbox).and.(monoxcloud).and.(.not.fast)) THEN
1115	call cocloud(ngrid,nlayer,naerkind,ptimestep,
1116	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
1117	& pq,pdq,zdqcocloud,zdqscocloud,zdtcocloud,
1118	& nq,rice_co,qsurf(1,igcm_n2),dqsurf(1,igcm_n2))
1119
1120	DO l=1,nlayer
1121	DO ig=1,ngrid
1122	pdq(ig,l,igcm_co_gas)=pdq(ig,l,igcm_co_gas)+
1123	& zdqcocloud(ig,l,igcm_co_gas)
1124	pdq(ig,l,igcm_co_ice)=pdq(ig,l,igcm_co_ice)+
1125	& zdqcocloud(ig,l,igcm_co_ice)
1126	ENDDO
1127	ENDDO
1128
1129	! Increment CO ice surface tracer tendency
1130	DO ig=1,ngrid
1131	dqsurf(ig,igcm_co_ice)=dqsurf(ig,igcm_co_ice)+
1132	& zdqscocloud(ig,igcm_co_ice)
1133	ENDDO
1134
1135	! update temperature tendancy
1136	DO ig=1,ngrid
1137	DO l=1,nlayer
1138	pdt(ig,l)=pdt(ig,l)+zdtcocloud(ig,l)
1139	ENDDO
1140	ENDDO
1141	ELSE
1142	rice_co(:,:)=0 ! initialization needed for callsedim
1143	END IF ! of IF (carbox)
1144
1145	!------------------------------------------------------------------
1146	! test methane conservation
1147	c if(methane)then
1148	c call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
1149	c & ptimestep,pplev,zdqch4cloud,zdqsch4cloud,'CH4','ch4clou')
1150	c endif ! methane
1151	!------------------------------------------------------------------
1152	! test CO conservation
1153	c if(carbox)then
1154	c call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
1155	c & ptimestep,pplev,zdqcocloud,zdqscocloud,'CO ','cocloud')
1156	c endif ! carbox
1157	!------------------------------------------------------------------
1158
1159	c 7b. Haze particle production
1160	c -------------------
1161	IF (haze) THEN
1162
1163	zdqphot_prec(:,:)=0.
1164	zdqphot_ch4(:,:)=0.
1165	zdqhaze(:,:,:)=0.
1166	! Forcing to a fixed haze profile if haze_proffix
1167	if (haze_proffix.and.i_haze.gt.0.) then
1168	call haze_prof(ngrid,nlayer,zzlay,pplay,pt,
1169	& reffrad,profmmr)
1170	zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze))
1171	& /ptimestep
1172	else
1173	call hazecloud(ngrid,nlayer,nq,ptimestep,
1174	& pplay,pplev,pq,pdq,dist_sol,mu0,zfluxuv,zdqhaze,
1175	& zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
1176	endif
1177
1178	DO iq=1, nq ! should be updated
1179	DO l=1,nlayer
1180	DO ig=1,ngrid
1181	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqhaze(ig,l,iq)
1182	ENDDO
1183	ENDDO
1184	ENDDO
1185
1186	ENDIF
1187
1188	IF (fast.and.fasthaze) THEN
1189	call prodhaze(ngrid,nlayer,nq,ptimestep,pplev,pq,pdq,dist_sol,
1190	& mu0,declin,zdqprodhaze,zdqsprodhaze,gradflux,fluxbot,
1191	& fluxlym_sol_bot,fluxlym_ipm_bot,flym_sol,flym_ipm)
1192
1193	DO ig=1,ngrid
1194	pdq(ig,1,igcm_ch4_gas)=pdq(ig,1,igcm_ch4_gas)+
1195	& zdqprodhaze(ig,igcm_ch4_gas)
1196	pdq(ig,1,igcm_prec_haze)=pdq(ig,1,igcm_prec_haze)+
1197	& zdqprodhaze(ig,igcm_prec_haze)
1198	pdq(ig,1,igcm_haze)=abs(pdq(ig,1,igcm_haze)+
1199	& zdqprodhaze(ig,igcm_haze))
1200	qsurf(ig,igcm_haze)= qsurf(ig,igcm_haze)+
1201	& zdqsprodhaze(ig)*ptimestep
1202	ENDDO
1203
1204	ENDIF
1205
1206	c 7c. Aerosol particles
1207	c -------------------
1208
1209	c -------------
1210	c Sedimentation
1211	c -------------
1212	IF (sedimentation) THEN
1213	call zerophys(ngridnlayernq, zdqsed)
1214	call zerophys(ngrid*nq, zdqssed)
1215
1216	call callsedim(ngrid,nlayer, ptimestep,
1217	& pplev,zzlev, pt,rice_ch4,rice_co,
1218	& pq, pdq, zdqsed, zdqssed,nq,pphi)
1219
1220	DO iq=1, nq
1221	DO l=1,nlayer
1222	DO ig=1,ngrid
1223	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
1224	ENDDO
1225	ENDDO
1226	ENDDO
1227	DO iq=2, nq
1228	DO ig=1,ngrid
1229	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
1230	ENDDO
1231	ENDDO
1232
1233	END IF ! of IF (sedimentation)
1234
1235	!------------------------------------------------------------------
1236	! test methane conservation
1237	c if(methane)then
1238	c call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
1239	c & ptimestep,pplev,zdqsed,zdqssed,'CH4',' sedim ')
1240	c endif ! methane
1241	!------------------------------------------------------------------
1242	! test CO conservation
1243	c if(carbox)then
1244	c call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
1245	c & ptimestep,pplev,zdqsed,zdqssed,'CO ',' sedim ')
1246	c endif ! carbox
1247	!------------------------------------------------------------------
1248
1249	c 7d. Updates
1250	c ---------
1251	! write(,) 'before update qsurf:'
1252	! call massn2(ngrid,nlayer,pplev,qsurf(:,1),dqsurf(:,igcm_n2),
1253	! & pdpsrf,ptimestep)
1254
1255	c ---------------------------------
1256	c Updating tracer budget on surface (before spread of glacier)
1257	c ---------------------------------
1258	DO iq=1, nq
1259	DO ig=1,ngrid
1260	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1261	ENDDO
1262	ENDDO
1263
1264	endif ! of if (tracer)
1265
1266	if (conservn2) then
1267	write(,) 'conservn2 tracer'
1268	call testconservmass(ngrid,nlayer,pplev(:,1)+
1269	& pdpsrf(:)*ptimestep,qsurf(:,1))
1270	endif
1271
1272	DO ig=1,ngrid
1273	flusurf(ig,igcm_n2)=(qsurf(ig,igcm_n2)-
1274	& qsurf1(ig,igcm_n2))/ptimestep
1275	flusurfold(ig,igcm_n2)=flusurf(ig,igcm_n2)
1276	if (methane) then
1277	flusurf(ig,igcm_ch4_ice)=(qsurf(ig,igcm_ch4_ice)-
1278	& qsurf1(ig,igcm_ch4_ice))/ptimestep
1279	flusurfold(ig,igcm_ch4_ice)=flusurf(ig,igcm_ch4_ice)
1280	endif
1281	if (carbox) then
1282	flusurf(ig,igcm_co_ice)=(qsurf(ig,igcm_co_ice)-
1283	& qsurf1(ig,igcm_co_ice))/ptimestep
1284	!flusurfold(ig,igcm_co_ice)=flusurf(ig,igcm_co_ice)
1285	endif
1286	ENDDO
1287
1288	!! Special source of haze particle !
1289	! todo: should be placed in haze and use tendency of n2 instead of flusurf
1290	IF (source_haze) THEN
1291	call hazesource(ngrid,nlayer,nq,ptimestep,
1292	& pplev,flusurf,mu0,zdq_source)
1293
1294	DO iq=1, nq
1295	DO l=1,nlayer
1296	DO ig=1,ngrid
1297	pdq(ig,l,iq)=pdq(ig,l,iq)+zdq_source(ig,l,iq)
1298	ENDDO
1299	ENDDO
1300	ENDDO
1301	ENDIF
1302
1303	!-----------------------------------------------------------------------
1304	! 9. Surface and sub-surface soil temperature
1305	!-----------------------------------------------------------------------
1306	! For diagnostic
1307	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1308	DO ig=1,ngrid
1309	rho(ig,1) = pplay(ig,1)/(r*pt(ig,1))
1310	sensiblehf1(ig)=rho(ig,1)cpp(0.4/log(zzlay(ig,1)/z0))*2
1311	& (pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1))*0.5
1312	& (tsurf(ig)-pt(ig,1))
1313	sensiblehf2(ig)=zflubid(ig)-capcal(ig)*zdtsdif(ig)
1314
1315	ENDDO
1316
1317	! 9.1 Increment Surface temperature:
1318	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1319	DO ig=1,ngrid
1320	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1321	ENDDO
1322
1323	! Prevent surface (.e.g. non volatile ch4) to exceed max temperature
1324	! Lellouch et al., 2000,2011
1325	IF (tsurfmax) THEN
1326	DO ig=1,ngrid
1327	if (albedo(ig).gt.albmin_ch4.and.
1328	& qsurf(ig,igcm_n2).eq.0.) then
1329	tsurf(ig)=min(tsurf(ig),54.)
1330	endif
1331	ENDDO
1332	ENDIF
1333	!
1334	! 9.2 Compute soil temperatures and subsurface heat flux:
1335	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1336	IF (callsoil) THEN
1337	CALL soil(ngrid,nsoilmx,.false.,tidat,
1338	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
1339	ENDIF
1340
1341	! For output :
1342	tidat_out(:,:)=0.
1343	DO l=1,min(nlayermx,nsoilmx)
1344	tidat_out(:,l)=tidat(:,l)
1345	ENDDO
1346
1347	! 9.3 multiply tendencies of cond/subli for paleo loop only in the
1348	! last Pluto year of the simulation
1349	! Year day must be adapted in the startfi for each object
1350	! Paleo uses year_day to calculate the annual mean tendancies
1351	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1352	IF (paleo) then
1353	if (zday.gt.day_ini+ptime0+nday-year_day) then
1354	DO iq=1,nq
1355	DO ig=1,ngrid
1356	qsurfyear(ig,iq)=qsurfyear(ig,iq)+
1357	& (qsurf(ig,iq)-qsurf1(ig,iq)) !kg m-2 !ptimestep
1358	ENDDO
1359	ENDDO
1360	endif
1361	endif
1362
1363	! 9.4 Glacial flow at each timestep glastep or at lastcall
1364	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1365	IF (fast.and.glaflow) THEN
1366	if ((mod(zday-day_ini-ptime0,glastep)).lt.1.
1367	& .or.lastcall) then
1368	IF (lastcall) then
1369	dstep=mod(zday-day_ini-ptime0,glastep)*daysec
1370	else
1371	dstep=glastep*daysec
1372	endif
1373	zdqflow(:,:)=qsurf(:,:)
1374	IF (paleo) then
1375	call spreadglacier_paleo(ngrid,nq,qsurf,
1376	& phisfinew,dstep,tsurf)
1377	else
1378	call spreadglacier_simple(ngrid,nq,qsurf,dstep)
1379	endif
1380	zdqflow(:,:)=(zdqflow(:,:)-qsurf(:,:))/dstep
1381
1382	if (conservn2) then
1383	write(,) 'conservn2 glaflow'
1384	call testconservmass(ngrid,nlayer,pplev(:,1)+
1385	& pdpsrf(:)*ptimestep,qsurf(:,1))
1386	endif
1387
1388	endif
1389	ENDIF
1390	!-----------------------------------------------------------------------
1391	! 10. Perform diagnostics and write output files
1392	!-----------------------------------------------------------------------
1393	! ---------------------------------------------------------
1394	! Check the energy balance of the simulation during the run
1395	! ---------------------------------------------------------
1396	flux1 = 0
1397	flux2 = 0
1398	flux3 = 0
1399	do ig=1,ngrid
1400	flux1 = flux1 +
1401	& area(ig)*(fluxtop_dn(ig) - fluxtop_sw(ig))
1402	flux2 = flux2 + area(ig)*fluxtop_lw(ig)
1403	flux3 = flux3 + area(ig)*fluxtop_dn(ig)
1404	fluxabs_sw(ig)=fluxtop_dn(ig) - fluxtop_sw(ig)
1405
1406	end do
1407	print*,'Incident solar flux, absorbed solar flux, OLR (W/m2)'
1408	print*, flux3/totarea, ' ', flux1/totarea ,
1409	& ' = ', flux2/totarea
1410
1411	if(meanOLR)then
1412	! to record global radiative balance
1413	open(92,file="rad_bal.out",form='formatted',access='append')
1414	write(92,*) zday,flux3/totarea,flux1/totarea,flux2/totarea
1415	close(92)
1416	endif
1417
1418	! Surface temperature information
1419	ts1 = 0
1420	ts2 = 999
1421	ts3 = 0
1422	do ig=1,ngrid
1423	ts1 = ts1 + area(ig)*tsurf(ig)
1424	ts2 = min(ts2,tsurf(ig))
1425	ts3 = max(ts3,tsurf(ig))
1426	end do
1427	! print*,'Mean Tsurf =',ts1/totarea , ' Min Tsurf=',ts2,
1428	! & ' Max Tsurf =',ts3
1429	! print*,'Max Tsurf=',ts3,'Min Tsurf=',ts2
1430
1431	! -------------------------------
1432	! Dynamical fields incrementation
1433	! -------------------------------
1434	! (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
1435	! temperature, zonal and meridional wind
1436	if(firstcall) ztprevious(:,:)=pt(:,:)
1437	DO l=1,nlayer
1438	DO ig=1,ngrid
1439	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
1440	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
1441	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
1442	! diagnostic
1443	zdtdyn(ig,l)=ztprevious(ig,l)-pt(ig,l)
1444	ztprevious(ig,l)=zt(ig,l)
1445	ENDDO
1446	ENDDO
1447	! if(firstcall)call zerophys(zdtdyn)
1448	! if(firstcall) zdtdyn(:,:)=0.
1449	! dynamical heating diagnostic
1450	DO ig=1,ngrid
1451	fluxdyn(ig)=0.
1452	DO l=1,nlayer
1453	fluxdyn(ig)=fluxdyn(ig) - (zdtdyn(ig,l)/ptimestep)
1454	& (pplev(ig,l)-pplev(ig,l+1))cpp/g
1455	ENDDO
1456	ENDDO
1457
1458	! -------------------------------
1459	! Other diagnostics
1460	! -------------------------------
1461	!! diagnostic qsat_lev1 (ps, T1)
1462	! CALL methanesat(ngridmx,zt(:,1),pplev(1,1),qsat_lev1(:),
1463	! & qsurf(:,igcm_n2))
1464
1465	!! eddy tracers : decay time
1466	! DO l=1,nlayer
1467	! DO ig=1,ngrid
1468	! pdq(ig,l,igcm_eddy1e6)=pdq(ig,l,igcm_eddy1e6)-
1469	! & pq(ig,l,igcm_eddy1e6)*(1-exp(-ptimestep/1.e6))/ptimestep
1470	! pdq(ig,l,igcm_eddy1e7)=pdq(ig,l,igcm_eddy1e7)-
1471	! & pq(ig,l,igcm_eddy1e7)*(1-exp(-ptimestep/1.e7))/ptimestep
1472	! pdq(ig,l,igcm_eddy5e7)=pdq(ig,l,igcm_eddy5e7)-
1473	! & pq(ig,l,igcm_eddy5e7)*(1-exp(-ptimestep/5.e7))/ptimestep
1474	! pdq(ig,l,igcm_eddy1e8)=pdq(ig,l,igcm_eddy1e8)-
1475	! & pq(ig,l,igcm_eddy1e8)*(1-exp(-ptimestep/1.e8))/ptimestep
1476	! pdq(ig,l,igcm_eddy5e8)=pdq(ig,l,igcm_eddy5e8)-
1477	! & pq(ig,l,igcm_eddy5e8)*(1-exp(-ptimestep/5.e8))/ptimestep
1478	! ENDDO
1479	! ENDDO
1480
1481	! -------------------------------
1482	! Update tracers
1483	! -------------------------------
1484	DO iq=1, nq
1485	DO l=1,nlayer
1486	DO ig=1,ngrid
1487	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
1488	ENDDO
1489	ENDDO
1490	ENDDO
1491
1492	DO ig=1,ngrid
1493	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
1494	ENDDO
1495	call globalaverage2d(ngrid,ps,globave)
1496
1497	! pressure density
1498	IF (.not.fast) then !
1499	DO ig=1,ngrid
1500	DO l=1,nlayer
1501	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1502	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1503	rho(ig,l) = zplay(ig,l)/(r*zt(ig,l))
1504	ENDDO
1505	zplev(ig,nlayer+1)=pplev(ig,nlayer+1)/pplev(ig,1)*ps(ig)
1506	ENDDO
1507	ENDIF
1508
1509	! ---------------------------------------------------------
1510	! Compute column amounts (kg m-2) if tracers are enabled
1511	! ---------------------------------------------------------
1512	call zerophys(ngrid*nq,qcol)
1513	if(tracer.and..not.fast)then
1514	do iq=1,nq
1515	DO ig=1,ngrid
1516	DO l=1,nlayer
1517	qcol(ig,iq) = qcol(ig,iq) + zq(ig,l,iq) *
1518	& (zplev(ig,l) - zplev(ig,l+1)) / g
1519	enddo
1520	enddo
1521	enddo
1522	endif
1523
1524	if (methane) then
1525	IF (fast) then ! zq is the mixing ratio supposingly mixed in all atmosphere
1526	DO ig=1,ngrid
1527	vmr_ch4(ig)=zq(ig,1,igcm_ch4_gas)*
1528	& mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.
1529	ENDDO
1530	ELSE
1531	DO ig=1,ngrid
1532	! compute vmr methane
1533	vmr_ch4(ig)=qcol(ig,igcm_ch4_gas)*
1534	& g/ps(ig)mmol(igcm_n2)/mmol(igcm_ch4_gas)100.
1535	! compute density methane
1536	DO l=1,nlayer
1537	zrho_ch4(ig,l)=zq(ig,l,igcm_ch4_gas)*rho(ig,l)
1538	ENDDO
1539	ENDDO
1540	ENDIF
1541	endif
1542
1543	if (carbox) then
1544	IF (fast) then
1545	DO ig=1,ngrid
1546	vmr_co(ig)=zq(ig,1,igcm_co_gas)*
1547	& mmol(igcm_n2)/mmol(igcm_co_gas)*100.
1548	ENDDO
1549	ELSE
1550	DO ig=1,ngrid
1551	! compute vmr CO
1552	vmr_co(ig)=qcol(ig,igcm_co_gas)*
1553	& g/ps(ig)mmol(igcm_n2)/mmol(igcm_co_gas)100.
1554	! compute density CO
1555	DO l=1,nlayer
1556	zrho_co(ig,l)=zq(ig,l,igcm_co_gas)*rho(ig,l)
1557	ENDDO
1558	ENDDO
1559	ENDIF
1560	endif
1561
1562	zrho_haze(:,:)=0.
1563	zdqrho_photprec(:,:)=0.
1564	IF (haze.and.aerohaze) then
1565	DO ig=1,ngrid
1566	DO l=1,nlayer
1567	zrho_haze(ig,l)=zq(ig,l,igcm_haze)*rho(ig,l)
1568	zdqrho_photprec(ig,l)=zdqphot_prec(ig,l)*rho(ig,l)
1569	ENDDO
1570	ENDDO
1571	ENDIF
1572
1573	IF (fasthaze) then
1574	DO ig=1,ngrid
1575	qcol(ig,igcm_haze)=zq(ig,1,igcm_haze)*pplev(ig,1)/g
1576	qcol(ig,igcm_prec_haze)=zq(ig,1,igcm_prec_haze)*pplev(ig,1)/g
1577	ENDDO
1578	ENDIF
1579
1580	c Info about Ls, declin...
1581	IF (fast) THEN
1582	write(,),'Ls=',zls180./pi,' dec=',declin180./pi
1583	write(,),'zday=',zday,' ps=',globave
1584	IF (lastcall) then
1585	write(,),'lastcall'
1586	ENDIF
1587	ELSE
1588	write(,),'Ls=',zls180./pi,'decli=',declin180./pi,'zday=',zday
1589	ENDIF
1590
1591	lecttsoil=0 ! default value for lecttsoil
1592	call getin("lecttsoil",lecttsoil)
1593	IF (lastcall.and.(ngrid.EQ.1).and.(lecttsoil.eq.1)) THEN
1594	! save tsoil temperature profile for 1D profile
1595	OPEN(13,file='proftsoil.out',form='formatted')
1596	DO i=1,nsoilmx
1597	write(13,*) tsoil(1,i)
1598	ENDDO
1599	CLOSE(13)
1600	ENDIF
1601
1602
1603	IF (ngrid.NE.1) THEN
1604	! -------------------------------------------------------------------
1605	! Writing NetCDF file "RESTARTFI" at the end of the run
1606	! -------------------------------------------------------------------
1607	! Note: 'restartfi' is stored just before dynamics are stored
1608	! in 'restart'. Between now and the writting of 'restart',
1609	! there will have been the itau=itau+1 instruction and
1610	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1611	! thus we store for time=time+dtvr
1612
1613	IF(lastcall) THEN
1614
1615	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1616	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1617
1618	if (paleo) then
1619	! time range for tendencies of ice flux qsurfyear
1620	zdt_tot=year_day ! Last year of simulation
1621
1622	masslost(:)=0.
1623	massacc(:)=0.
1624
1625	DO ig=1,ngrid
1626	! update new reservoir of ice on the surface
1627	DO iq=1,nq
1628	! kg/m2 to be sublimed or condensed during paleoyears
1629	qsurfyear(ig,iq)=qsurfyear(ig,iq)*
1630	& paleoyears365.25/(zdt_totdaysec/86400.)
1631
1632	! special case if we sublime the entire reservoir
1633	IF (-qsurfyear(ig,iq).gt.qsurf(ig,iq)) THEN
1634	masslost(iq)=masslost(iq)+realarea(ig)*
1635	& (-qsurfyear(ig,iq)-qsurf(ig,iq))
1636	qsurfyear(ig,iq)=-qsurf(ig,iq)
1637	ENDIF
1638
1639	IF (qsurfyear(ig,iq).gt.0.) THEN
1640	massacc(iq)=massacc(iq)+realarea(ig)*qsurfyear(ig,iq)
1641	ENDIF
1642
1643	ENDDO
1644	ENDDO
1645
1646	DO ig=1,ngrid
1647	DO iq=1,nq
1648	qsurfpal(ig,iq)=qsurf(ig,iq)+qsurfyear(ig,iq)
1649	IF (qsurfyear(ig,iq).gt.0.) THEN
1650	qsurfpal(ig,iq)=qsurfpal(ig,iq)-
1651	& qsurfyear(ig,iq)*masslost(iq)/massacc(iq)
1652	ENDIF
1653	ENDDO
1654	ENDDO
1655	! Finally ensure conservation of qsurf
1656	DO iq=1,nq
1657	call globalaverage2d(ngrid,qsurf(:,iq),globaveice(iq))
1658	call globalaverage2d(ngrid,qsurfpal(:,iq),
1659	& globavenewice(iq))
1660	IF (globavenewice(iq).gt.0.) THEN
1661	qsurfpal(:,iq)=qsurfpal(:,iq)*
1662	& globaveice(iq)/globavenewice(iq)
1663	ENDIF
1664	ENDDO
1665
1666	! update new geopotential depending on the ice reservoir
1667	phisfipal(:)=phisfinew(:)+qsurfpal(:,igcm_n2)*g/1000.
1668	!phisfipal(ig)=phisfi(ig)
1669
1670
1671	if (kbo.or.triton) then ! case of Triton : we do not change the orbital parameters
1672
1673	pdaypal=pday ! no increment of pdaypal to keep same evolution of the subsolar point
1674	eccpal=1.-periheli/((periheli+aphelie)/2.) !no change of ecc
1675	peri_daypal=peri_day ! no change
1676	oblipal=obliquit ! no change
1677	tpalnew=tpal
1678	adjustnew=adjust
1679
1680	else ! Pluto
1681	! update new pday and tpal (Myr) to be set in startfi controle
1682	pdaypal=int(day_ini+paleoyears*365.25/6.3872)
1683	tpalnew=tpal+paleoyears*1.e-6 ! Myrs
1684
1685	! update new N2 ice adjustment (not tested yet on Pluto)
1686	adjustnew=adjust
1687
1688	! update milankovitch parameters : obliquity,Lsp,ecc
1689	call calcmilank(tpalnew,oblipal,peri_daypal,eccpal)
1690	!peri_daypal=peri_day
1691	!eccpal=0.009
1692
1693	endif
1694
1695	write(,) "Paleo peri=",peri_daypal," tpal=",tpalnew
1696	write(,) "Paleo eccpal=",eccpal," tpal=",tpalnew
1697
1698	! create restartfi
1699	call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq,
1700	. ptimestep,pdaypal,
1701	. ztime_fin,tsurf,tsoil,emis,q2,qsurfpal,
1702	. area,albedodat,tidat,zmea,zstd,zsig,
1703	. zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal,
1704	. peri_daypal)
1705	else
1706
1707	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
1708	. ptimestep,pday,
1709	. ztime_fin,tsurf,tsoil,emis,q2,qsurf,
1710	. area,albedodat,tidat,zmea,zstd,zsig,
1711	. zgam,zthe)
1712
1713	endif
1714
1715	ENDIF
1716
1717	! -----------------------------------------------------------------
1718	! Saving statistics :
1719	! -----------------------------------------------------------------
1720	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
1721	! which can later be used to make the statistic files of the run:
1722	! "stats") only possible in 3D runs !
1723
1724
1725	IF (callstats) THEN
1726	write (,) "stats have been turned off in the code.
1727	& You can manage your own output in physiq.F "
1728	stop
1729
1730	! call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1731	! call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1732	! call wstats(ngrid,"fluxsurf_lw",
1733	! . "Thermal IR radiative flux to surface","W.m-2",2,
1734	! . fluxsurf_lw)
1735	! call wstats(ngrid,"fluxsurf_sw",
1736	! . "Solar radiative flux to surface","W.m-2",2,
1737	! . fluxsurf_sw_tot)
1738	! call wstats(ngrid,"fluxtop_lw",
1739	! . "Thermal IR radiative flux to space","W.m-2",2,
1740	! . fluxtop_lw)
1741	! call wstats(ngrid,"fluxtop_sw",
1742	! . "Solar radiative flux to space","W.m-2",2,
1743	! . fluxtop_sw_tot)
1744
1745	! call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1746	! call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1747	! call wstats(ngrid,"v","Meridional (North-South) wind",
1748	! . "m.s-1",3,zv)
1749
1750	IF(lastcall) THEN
1751	write (,) "Writing stats..."
1752	call mkstats(ierr)
1753	ENDIF
1754	ENDIF !if callstats
1755
1756
1757	! ---------------------------------------------------------------------
1758	! 3D OUTPUTS
1759	! ----------------------------------------------------------------------
1760	! output in netcdf file "DIAGFI", containing any variable for diagnostic
1761	! (output with period "ecritphy", set in "run.def")
1762	! ----------------------------------------------------------------------
1763
1764	! if(MOD(countG3D,saveG3D).eq.0)then
1765	! print*,'countG3D',countG3D
1766
1767	!! BASIC
1768
1769	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
1770	& tsurf)
1771	call WRITEDIAGFI(ngrid,"rice_ch4","ch4 ice mass mean radius",
1772	& "m",3,rice_ch4)
1773	call WRITEDIAGFI(ngrid,"area","area","m",2,area)
1774	call WRITEDIAGFI(ngrid,"mu0","cos sza","m",2,mu0)
1775	call WRITEDIAGFI(ngrid,"fract","fractional day","",2,fract)
1776	call WRITEDIAGFI(ngrid,"declin","declin","deg",0,declin*180./pi)
1777	call WRITEDIAGFI(ngrid,"ls","ls","deg",0,zls*180./pi)
1778	call WRITEDIAGFI(ngrid,"dist_sol","dist_sol","AU",0,dist_sol)
1779	call WRITEDIAGFI(ngrid,"phisfi","phisfi"," ",2,phisfi)
1780	! call WRITEDIAGFI(ngrid,"tsub","subsurface temperature","K",
1781	! & 1,tsub)
1782	! call WRITEDIAGFI(ngrid,"phitop","phitop"," ",2,phitop)
1783
1784	! ENERGY - Total energy balance diagnostics
1785
1786	call WRITEDIAGFI(ngrid,"albedo","albedo","sans dim",2,
1787	& albedo)
1788	call WRITEDIAGFI(ngrid,"emissivite","emissivite","sans dim",2,
1789	& emis)
1790	call WRITEDIAGFI(ngrid,"fluxtop_dn","fluxtop_dn","sans dim",2,
1791	& fluxtop_dn)
1792	call WRITEDIAGFI(ngrid,"ISR","incoming stellar rad.","W m-2",
1793	& 2,fluxtop_dn)
1794	call WRITEDIAGFI(ngrid,"ASR","absorbed stellar rad.","W m-2",
1795	& 2,fluxabs_sw)
1796	call WRITEDIAGFI(ngrid,"OLR","outgoing longwave rad.","W m-2",
1797	& 2,fluxtop_lw)
1798
1799	IF (fast) then
1800	! pressure reprocessed in nogcm.F
1801	call WRITEDIAGFI(ngrid,"globave","surf press","Pa",0,globave)
1802	call WRITEDIAGFI(ngrid,"ps","surf press","Pa",2,ps)
1803	call WRITEDIAGFI(ngrid,"fluxrad","fluxrad",
1804	& "W m-2",2,fluxrad)
1805	call WRITEDIAGFI(ngrid,"fluxgrd","fluxgrd",
1806	& "W m-2",2,fluxgrd)
1807	call WRITEDIAGFI(ngrid,"capcal","capcal",
1808	& "W.s m-2 K-1",2,capcal)
1809	! call WRITEDIAGFI(ngrid,"dplanck","dplanck",
1810	! & "W.s m-2 K-1",2,dplanck)
1811	call WRITEDIAGFI(ngrid,"tsoil","tsoil","K",3,tsoil)
1812
1813	! If we want to see tidat evolution with time:
1814	call WRITEDIAGFI(ngrid,"tidat","tidat","SI",3,tidat_out)
1815
1816
1817	IF (fasthaze) then
1818	call WRITEDIAGFI(ngrid,"gradflux","gradflux",
1819	& "Ph m-2 s-1",2,gradflux)
1820	call WRITEDIAGFI(ngrid,"fluxbot","flux bottom",
1821	& "Ph m-2 s-1",2,fluxbot)
1822	call WRITEDIAGFI(ngrid,"fluxbotsol","flux bottom sol",
1823	& "Ph m-2 s-1",2,fluxlym_sol_bot)
1824	call WRITEDIAGFI(ngrid,"fluxbotipm","flux bottom ipm",
1825	& "Ph m-2 s-1",2,fluxlym_ipm_bot)
1826	call WRITEDIAGFI(ngrid,"fluxlymipm","flux incid ipm",
1827	& "Ph m-2 s-1",2,flym_ipm)
1828	call WRITEDIAGFI(ngrid,"fluxlymsol","flux incid sol",
1829	& "Ph m-2 s-1",2,flym_sol)
1830	call WRITEDIAGFI(ngrid,"tend_prodhaze","prod haze",
1831	& "kg m-2 s-1",2,zdqsprodhaze)
1832	call WRITEDIAGFI(ngrid,"tend_lostch4","tend ch4",
1833	& "kg kg-1 s-1",2,zdqprodhaze(1,igcm_ch4_gas))
1834	call WRITEDIAGFI(ngrid,"tend_prodhaze","tend prod haze",
1835	& "kg kg-1 s-1",2,zdqprodhaze(1,igcm_haze))
1836	ENDIF
1837	IF (paleo) then
1838	call WRITEDIAGFI(ngrid,"qsurfn2_year","qsurfn2_year",
1839	& "kg m-2.plutoyear-1",2,qsurfyear(:,1))
1840	call WRITEDIAGFI(ngrid,"phisfipal","phisfipal",
1841	& " ",2,phisfipal)
1842	call WRITEDIAGFI(ngrid,"zdqflow","zdqflow",
1843	& " ",2,zdqflow(1,igcm_n2))
1844	ENDIF
1845	ELSE
1846	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
1847	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
1848	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1849	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
1850	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
1851	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
1852	call WRITEDIAGFI(ngrid,"p","Pression","Pa",3,pplay)
1853	call WRITEDIAGFI(ngrid,"fluxrad","fluxrad",
1854	& "W m-2",2,fluxrad)
1855	call WRITEDIAGFI(ngrid,"fluxgrd","fluxgrd",
1856	& "W m-2",2,fluxgrd)
1857	ENDIF
1858
1859	! Usually not used :
1860
1861	! call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
1862	call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
1863	! call WRITEDIAGFI(ngrid,"DYN","dynamical heat input","W m-2",
1864	! & 2,fluxdyn)
1865
1866	! TENDANCIES
1867
1868	call WRITEDIAGFI(ngrid,"dps","dps","K",2,pdpsrf)
1869	call WRITEDIAGFI(ngrid,"zdtsw","SW heating","K s-1",3,zdtsw)
1870	call WRITEDIAGFI(ngrid,"zdtlw","LW heating","K s-1",3,zdtlw)
1871	call WRITEDIAGFI(ngrid,"zdtc","tendancy T cond N2","K",3,zdtc)
1872	call WRITEDIAGFI(ngrid,"zdtch4cloud","tendancy T ch4cloud",
1873	& "K",3,zdtch4cloud)
1874	call WRITEDIAGFI(ngrid,"zdtcocloud","tendancy T cocloud",
1875	& "K",3,zdtcocloud)
1876	call WRITEDIAGFI(ngrid,"zdtsurfc","zdtsurfc","K",2,zdtsurfc)
1877	call WRITEDIAGFI(ngrid,"zdtdif","zdtdif","K",3,zdtdif)
1878	call WRITEDIAGFI(ngrid,"zdtconduc","tendancy T conduc",
1879	& "K",3,zdtconduc)
1880	call WRITEDIAGFI(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
1881	call WRITEDIAGFI(ngrid,"zdtrad","rad heating","T s-1",3,dtrad)
1882	call WRITEDIAGFI(ngrid,"zdtadj","conv adj","T s-1",3,zdtadj)
1883	! call WRITEDIAGFI(ngrid,"zdqcloud","ch4 cloud","T s-1",
1884	! & 3,zdqch4cloud(1,1,igcm_ch4_ice))
1885	! call WRITEDIAGFI(ngrid,"zdqcloudgas","ch4 cloud","T s-1",
1886	! & 3,zdqch4cloud(1,1,igcm_ch4_gas))
1887	! call WRITEDIAGFI(ngrid,"zdqcn2_c","zdq condn2","",3,zdqc(:,:,3))
1888	! call WRITEDIAGFI(ngrid,"zdqdif_c","zdqdif","",3,zdqdif(:,:,3))
1889	! call WRITEDIAGFI(ngrid,"zdqadj_c","zdqadj","",3,zdqadj(:,:,3))
1890	!call WRITEDIAGFI(ngrid,"zdqsed_h","zdqsed","",3,zdqsed(:,:,7))
1891	!call WRITEDIAGFI(ngrid,"zdqssed","zdqssed","",2,zdqssed)
1892	call WRITEDIAGFI(ngrid,"zq1temp_ch4"," "," ",2,zq1temp_ch4)
1893	call WRITEDIAGFI(ngrid,"qsat_ch4"," "," ",2,qsat_ch4)
1894	call WRITEDIAGFI(ngrid,"qsat_ch4_l1"," "," ",2,qsat_ch4_l1)
1895	call WRITEDIAGFI(ngrid,"senshf1","senshf1"," ",2,sensiblehf1)
1896	call WRITEDIAGFI(ngrid,"senshf2","senshf2"," ",2,sensiblehf2)
1897
1898
1899	! OTHERS
1900
1901	! call WRITEDIAGFI(ngrid,"dWtotdifv","dWtotdifv","kg m-2",1,dWtot)
1902	! call WRITEDIAGFI(ngrid,"dWtotsdifv","dWtotsdifv","kgm-2",1,dWtots)
1903	! call WRITEDIAGFI(ngrid,"nconsMX","nconsMX","kgm-2s-1",1,nconsMAX)
1904
1905
1906	! TRACERS
1907
1908	if (tracer) then
1909
1910	!!! afficher la tendance sur le flux de la glace
1911	call WRITEDIAGFI(ngridmx,'n2_iceflux',
1912	& 'n2_iceflux',
1913	& "kg m^-2 s^-1",2,flusurf(1,igcm_n2) )
1914
1915	do iq=1,nq
1916	call WRITEDIAGFI(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
1917	call WRITEDIAGFI(ngrid,trim(noms(iq))//'_Tcondn2',
1918	& trim(noms(iq))//'_Tcondn2',
1919	& 'kg/kg',3,zdqc(1,1,iq))
1920
1921	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_surf',
1922	& trim(noms(iq))//'_surf',
1923	& 'kg m^-2',2,qsurf(1,iq) )
1924
1925	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_col',
1926	& trim(noms(iq))//'_col',
1927	& 'kg m^-2',2,qcol(1,iq) )
1928
1929	! call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_aero',
1930	! & trim(noms(iq))//'_aero',
1931	! & 'kg/kg',3,aerosol(1,1,iq))
1932
1933
1934	enddo
1935
1936	call WRITEDIAGFI(ngridmx,'haze_reff',
1937	& 'haze_reff',
1938	& 'm',3,reffrad(1,1,1))
1939
1940	if (methane) then
1941
1942	call WRITEDIAGFI(ngridmx,'ch4_iceflux',
1943	& 'ch4_iceflux',
1944	& "kg m^-2 s^-1",2,flusurf(1,igcm_ch4_ice) )
1945
1946	call WRITEDIAGFI(ngrid,"vmr_ch4","vmr_ch4","%",
1947	& 2,vmr_ch4)
1948	if (.not.fast) THEN
1949	call WRITEDIAGFI(ngrid,"zrho_ch4","zrho_ch4","kg.m-3",
1950	& 3,zrho_ch4(:,:))
1951	endif
1952
1953	! if (fast) THEN
1954	! call WRITEDIAGFI(ngrid,"zq1_ch4","zq ch4","kg m-2",
1955	! & 2,zq(1,1,igcm_ch4_gas))
1956	! endif
1957
1958	! Tendancies
1959	call WRITEDIAGFI(ngrid,"zdqch4cloud","ch4 cloud","T s-1",
1960	& 3,zdqch4cloud(1,1,igcm_ch4_gas))
1961	call WRITEDIAGFI(ngrid,"zdqcn2_ch4","zdq condn2 ch4","",
1962	& 3,zdqc(:,:,igcm_ch4_gas))
1963	call WRITEDIAGFI(ngrid,"zdqdif_ch4","zdqdif ch4","",
1964	& 3,zdqdif(:,:,igcm_ch4_gas))
1965	call WRITEDIAGFI(ngrid,"zdqsdif_ch4_ice","zdqsdif ch4","",
1966	& 2,zdqsdif(:,igcm_ch4_ice))
1967	call WRITEDIAGFI(ngrid,"zdqadj_ch4","zdqadj ch4","",
1968	& 3,zdqadj(:,:,igcm_ch4_gas))
1969	call WRITEDIAGFI(ngrid,"zdqsed_ch4","zdqsed ch4","",
1970	& 3,zdqsed(:,:,igcm_ch4_gas))
1971	call WRITEDIAGFI(ngrid,"zdqssed_ch4","zdqssed ch4","",
1972	& 2,zdqssed(:,igcm_ch4_gas))
1973	call WRITEDIAGFI(ngrid,"zdtch4cloud","ch4 cloud","T s-1",
1974	& 3,zdtch4cloud)
1975
1976	endif
1977
1978	if (carbox) then
1979
1980	call WRITEDIAGFI(ngridmx,'co_iceflux',
1981	& 'co_iceflux',
1982	& "kg m^-2 s^-1",2,flusurf(1,igcm_co_ice) )
1983
1984	call WRITEDIAGFI(ngrid,"vmr_co","vmr_co","%",
1985	& 2,vmr_co)
1986
1987	if (.not.fast) THEN
1988	call WRITEDIAGFI(ngrid,"zrho_co","zrho_co","kg.m-3",
1989	& 3,zrho_co(:,:))
1990	endif
1991
1992	! if (fast) THEN
1993	! call WRITEDIAGFI(ngrid,"zq1_co","zq co","kg m-2",
1994	! & 2,zq(1,1,igcm_co_gas))
1995	! endif
1996
1997	! Tendancies
1998	! call WRITEDIAGFI(ngrid,"zdqcocloud","co cloud","T s-1",
1999	! & 3,zdqcocloud(1,1,igcm_co_gas))
2000	! call WRITEDIAGFI(ngrid,"zdqcn2_co","zdq condn2 co","",
2001	! & 3,zdqc(:,:,igcm_co_gas))
2002	! call WRITEDIAGFI(ngrid,"zdqdif_co","zdqdif co","",
2003	! & 3,zdqdif(:,:,igcm_co_gas))
2004	! call WRITEDIAGFI(ngrid,"zdqsdif_co_ice","zdqsdif co","",
2005	! & 2,zdqsdif(:,igcm_co_ice))
2006	! call WRITEDIAGFI(ngrid,"zdqadj_co","zdqadj co","",
2007	! & 3,zdqadj(:,:,igcm_co_gas))
2008	! call WRITEDIAGFI(ngrid,"zdqsed_co","zdqsed co","",
2009	! & 3,zdqsed(:,:,igcm_co_gas))
2010	! call WRITEDIAGFI(ngrid,"zdqssed_co","zdqssed co","",
2011	! & 2,zdqssed(:,igcm_co_gas))
2012	! call WRITEDIAGFI(ngrid,"zdtcocloud","co cloud","T s-1",
2013	! & 3,zdtcocloud)
2014
2015	endif
2016
2017	if (haze) then
2018
2019	! call WRITEDIAGFI(ngrid,"zrho_haze","zrho_haze","kg.m-3",
2020	! & 3,zrho_haze(:,:))
2021	call WRITEDIAGFI(ngrid,"zdqrho_photprec","zdqrho_photprec",
2022	& "kg.m-3.s-1",3,zdqrho_photprec(:,:))
2023	call WRITEDIAGFI(ngrid,"zdqphot_prec","zdqphot_prec","",
2024	& 3,zdqphot_prec(:,:))
2025	call WRITEDIAGFI(ngrid,"zdqhaze_ch4","zdqhaze_ch4","",
2026	& 3,zdqhaze(:,:,igcm_ch4_gas))
2027	call WRITEDIAGFI(ngrid,"zdqhaze_prec","zdqhaze_prec","",
2028	& 3,zdqhaze(:,:,igcm_prec_haze))
2029	if (igcm_haze.ne.0) then
2030	call WRITEDIAGFI(ngrid,"zdqhaze_haze","zdqhaze_haze","",
2031	& 3,zdqhaze(:,:,igcm_haze))
2032	call WRITEDIAGFI(ngrid,"zdqssed_haze","zdqssed haze",
2033	& "kg/m2/s",2,zdqssed(:,igcm_haze))
2034	endif
2035	call WRITEDIAGFI(ngrid,"zdqphot_ch4","zdqphot_ch4","",
2036	& 3,zdqphot_ch4(:,:))
2037	call WRITEDIAGFI(ngrid,"zdqconv_prec","zdqconv_prec","",
2038	& 3,zdqconv_prec(:,:))
2039	! call WRITEDIAGFI(ngrid,"zdqhaze_col","zdqhaze col","kg/m2/s",
2040	! & 2,zdqhaze_col(:))
2041	endif
2042
2043	if (aerohaze) then
2044	call WRITEDIAGFI(ngridmx,"tau_col",
2045	& "Total aerosol optical depth","opacity",2,tau_col)
2046	endif
2047
2048	! call WRITEDIAGFI(ngridmx,"tau_col",
2049	! & "Total aerosol optical depth","[]",2,tau_col)
2050	endif
2051
2052	! countG3D=1
2053	! else
2054	! countG3D=countG3D+1
2055	! endif ! if time to save
2056
2057	! ----------------------------------------------------------------------
2058	! 1D OUTPUTS
2059	! ----------------------------------------------------------------------
2060	ELSE ! if(ngrid.eq.1)
2061
2062	if(countG1D.eq.saveG1D)then
2063
2064	! BASIC 1D ONLY
2065
2066	call WRITEDIAGFI(ngrid,"ISR","incoming stellar rad.","W m-2",
2067	& 0,fluxtop_dn)
2068	call WRITEDIAGFI(ngrid,"ASR","absorbed stellar rad.","W m-2",
2069	& 0,fluxabs_sw)
2070	call WRITEDIAGFI(ngrid,"OLR","outgoing longwave rad.","W m-2",
2071	& 0,fluxtop_lw)
2072	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
2073	& tsurf)
2074	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",0,ps)
2075	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
2076	call WRITEDIAGFI(ngrid,"p","Pression","Pa",3,pplay)
2077
2078	call WRITEDIAGFI(ngrid,"fluxsurf_sw","sw surface flux",
2079	& "K",0,fluxsurf_sw)
2080	call WRITEDIAGFI(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
2081	call WRITEDIAGFI(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
2082	call WRITEDIAGFI(ngrid,"zdtdif","zdtdif","K",3,zdtdif)
2083	call WRITEDIAGFI(ngrid,"zdtconduc","tendancy T conduc",
2084	& "K",3,zdtconduc)
2085
2086	call WRITEDIAGFI(ngrid,"fluxsurf_sw","sw surface flux",
2087	& "K",0,fluxsurf_sw)
2088	call WRITEDIAGFI(ngrid,"declin","lat subsolaire",
2089	& "deg",0,declin*180./pi)
2090	call WRITEDIAGFI(ngrid,"ls","ls","deg",0,zls*180./pi)
2091	call WRITEDIAGFI(ngrid,"dist_sol","dist_sol","AU",0,dist_sol)
2092	call WRITEDIAGFI(ngrid,"mu0","solar zenith angle",
2093	& "deg",0,mu0)
2094	call WRITEDIAGFI(ngrid,"albedo","albedo",
2095	& "",0,albedo)
2096	call WRITEDIAGFI(ngrid,"emis","emis",
2097	& "",0,emis)
2098	call WRITEDIAGFI(ngrid,"tsoil","tsoil","K",1,tsoil(1,:))
2099	call WRITEDIAGFI(ngrid,"tidat","tidat","SI",1,tidat_out(1,:))
2100
2101	! OUTPUT OF TENDANCIES
2102	! call WRITEDIAGFI(ngrid,"zdqcloud","ch4 cloud","T s-1",
2103	! & 3,zdqcloud(1,1,igcm_ch4_gas))
2104	! call WRITEDIAGFI(ngrid,"zdqcn2","zdq condn2","",1,zdqc(:,:,3))
2105	! call WRITEDIAGFI(ngrid,"zdqdif","zdqdif","",1,zdqdif(:,:,3))
2106	! call WRITEDIAGFI(ngrid,"zdqadj","zdqadj","",1,zdqadj(:,:,3))
2107	! call WRITEDIAGFI(ngrid,"zdqsed","zdqsed","",1,zdqsed(:,:,3))
2108	! call WRITEDIAGFI(ngrid,"zdqc","zdqc","",1,zdqc(:,:,3))
2109	! call WRITEDIAGFI(ngrid,"zdqhaze","zdqhaze","",1,
2110	! & zdqhaze(:,:,3))
2111
2112	! 1D methane
2113	if (methane) then
2114
2115	call WRITEDIAGFI(ngridmx,'ch4_iceflux',
2116	& 'ch4_iceflux',
2117	& "kg m^-2 s^-1",0,flusurf(1,igcm_ch4_ice) )
2118
2119	call WRITEDIAGFI(ngrid,"vmr_ch4","vmr_ch4","%",
2120	& 0,vmr_ch4(1))
2121	call WRITEDIAGFI(ngrid,"zrho_ch4","zrho_ch4","kg.m-3",
2122	& 1,zrho_ch4(1,:))
2123
2124	! Tendancies
2125	call WRITEDIAGFI(ngrid,"zdqch4cloud","ch4 cloud","T s-1",
2126	& 1,zdqch4cloud(1,1,igcm_ch4_gas))
2127	call WRITEDIAGFI(ngrid,"zdqcn2_ch4","zdq condn2 ch4","",
2128	& 1,zdqc(1,:,igcm_ch4_gas))
2129	call WRITEDIAGFI(ngrid,"zdqdif_ch4","zdqdif ch4","",
2130	& 1,zdqdif(1,:,igcm_ch4_gas))
2131	call WRITEDIAGFI(ngrid,"zdqsdif_ch4_ice","zdqsdif ch4","",
2132	& 1,zdqsdif(:,igcm_ch4_ice))
2133	call WRITEDIAGFI(ngrid,"zdqadj_ch4","zdqadj ch4","",
2134	& 1,zdqadj(1,:,igcm_ch4_gas))
2135	call WRITEDIAGFI(ngrid,"zdqsed_ch4","zdqsed ch4","",
2136	& 1,zdqsed(1,:,igcm_ch4_gas))
2137	call WRITEDIAGFI(ngrid,"zdqssed_ch4","zdqssed ch4","",
2138	& 0,zdqssed(1,igcm_ch4_gas))
2139
2140	endif
2141
2142	! 1D Haze
2143	if (haze) then
2144
2145	! call WRITEDIAGFI(ngrid,"zrho_haze","zrho_haze","kg.m-3",
2146	! & 1,zrho_haze(:,:))
2147	call WRITEDIAGFI(ngrid,"zdqrho_photprec","zdqrho_photprec",
2148	& "kg.m-3.s-1",3,zdqrho_photprec)
2149	call WRITEDIAGFI(ngrid,"zdqphot_prec","zdqphot_prec","",
2150	& 3,zdqphot_prec)
2151	call WRITEDIAGFI(ngrid,"zdqhaze_ch4","zdqhaze_ch4","",
2152	& 1,zdqhaze(1,:,igcm_ch4_gas))
2153	call WRITEDIAGFI(ngrid,"zdqhaze_prec","zdqhaze_prec","",
2154	& 1,zdqhaze(1,:,igcm_prec_haze))
2155	call WRITEDIAGFI(ngrid,"zdqhaze_haze","zdqhaze_haze","",
2156	& 1,zdqhaze(1,:,igcm_haze))
2157	call WRITEDIAGFI(ngrid,"zdqphot_ch4","zdqphot_ch4","",
2158	& 3,zdqphot_ch4)
2159	! call WRITEDIAGFI(ngrid,"zdqconv_prec","zdqconv_prec","",
2160	! & 1,zdqconv_prec(1,:))
2161	call WRITEDIAGFI(ngrid,"zdqssed_haze","zdqssed haze","kg/m2/s",
2162	& 0,zdqssed(1,igcm_haze))
2163	! call WRITEDIAGFI(ngrid,"zdqhaze_col","zdqhaze col","kg/m2/s",
2164	! & 0,zdqhaze_col(:))
2165	if (haze_radproffix) then
2166	call WRITEDIAGFI(ngridmx,'haze_reffrad','haze_reffrad','m',
2167	& 3,reffrad(1,1,1))
2168	endif
2169	endif
2170
2171	if (aerohaze) then
2172	call WRITEDIAGFI(ngridmx,"tau_col",
2173	& "Total aerosol optical depth","[]",0,tau_col(1))
2174	call WRITEDIAGFI(ngridmx,"tau_aero",
2175	& "aerosol optical depth","[]",3,aerosol(1,1,1))
2176	endif
2177
2178	! TRACERS 1D ONLY
2179	if (tracer) then
2180	do iq=1,nq
2181	call WRITEDIAGFI(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
2182	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_surf',
2183	& trim(noms(iq))//'_surf',
2184	& 'kg m^-2',2,qsurf(1,iq) )
2185	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_col',
2186	& trim(noms(iq))//'_col',
2187	& 'kg m^-2',2,qcol(1,iq) )
2188	enddo
2189	endif
2190
2191	countG1D=1
2192	else
2193	countG1D=countG1D+1
2194	endif ! if time to save
2195
2196	END IF ! if(ngrid.ne.1)
2197
2198	! countG3D=countG3D+1
2199	icount=icount+1
2200
2201	RETURN
2202	END

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

Download in other formats:

Original Format