Context Navigation

physiq.F_bugT

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

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

Download in other formats:

Original Format