Context Navigation

physiq_mod.F90 @ 1808

Last change on this file since 1808 was 1808, checked in by jvatant, 7 years ago
Minor modifs : + correct a bug for ifort compiling in muphy + added some outputs for chemistry JVO
Property svn:executable set to ``*
File size: 62.7 KB

Line
1	module physiq_mod
2
3	implicit none
4
5	contains
6
7	subroutine physiq(ngrid,nlayer,nq, &
8	nametrac, &
9	firstcall,lastcall, &
10	pday,ptime,ptimestep, &
11	pplev,pplay,pphi, &
12	pu,pv,pt,pq, &
13	flxw, &
14	pdu,pdv,pdt,pdq,pdpsrf)
15
16	use radinc_h, only : L_NSPECTI,L_NSPECTV
17	use radcommon_h, only: sigma, glat, grav, BWNV
18	use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe
19	use comdiurn_h, only: coslat, sinlat, coslon, sinlon
20	use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
21	use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
22	use geometry_mod, only: latitude, longitude, cell_area
23	USE comgeomfi_h, only: totarea, totarea_planet
24	USE tracer_h
25	use time_phylmdz_mod, only: ecritphy, iphysiq, nday
26	use phyetat0_mod, only: phyetat0
27	use phyredem, only: physdem0, physdem1
28	use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
29	use mod_phys_lmdz_para, only : is_master
30	use planete_mod, only: apoastr, periastr, year_day, peri_day, &
31	obliquit, nres, z0
32	use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
33	use time_phylmdz_mod, only: daysec
34	use logic_mod, only: moyzon_ch, moyzon_mu
35	use moyzon_mod, only: tmoy, playmoy, zphibar, zphisbar, zplevbar, &
36	zplaybar, zzlevbar, zzlaybar, ztfibar
37	use callkeys_mod
38	use vertical_layers_mod, only: presnivs, pseudoalt
39	#ifdef CPP_XIOS
40	use xios_output_mod, only: initialize_xios_output, &
41	update_xios_timestep, &
42	send_xios_field
43	#endif
44	implicit none
45
46
47	!==================================================================
48	!
49	! Purpose
50	! -------
51	! Central subroutine for all the physics parameterisations in the
52	! universal model. Originally adapted from the Mars LMDZ model.
53	!
54	! The model can be run without or with tracer transport
55	! depending on the value of "tracer" in file "callphys.def".
56	!
57	!
58	! It includes:
59	!
60	! I. Initialization :
61	! I.1 Firstcall initializations.
62	! I.2 Initialization for every call to physiq.
63	!
64	! II. Compute radiative transfer tendencies (longwave and shortwave) :
65	! II.a Option 1 : Call correlated-k radiative transfer scheme.
66	! II.b Option 2 : Call Newtonian cooling scheme.
67	! II.c Option 3 : Atmosphere has no radiative effect.
68	!
69	! III. Vertical diffusion (turbulent mixing) :
70	!
71	! IV. Dry Convective adjustment :
72	!
73	! V. Tracers
74	! V.1. Chemistry
75	! V.2. Microphysics
76	! V.3. Updates (pressure variations, surface budget).
77	! V.4. Surface Tracer Update.
78	!
79	! VI. Surface and sub-surface soil temperature.
80	!
81	! VII. Perform diagnostics and write output files.
82	!
83	!
84	! arguments
85	! ---------
86	!
87	! INPUT
88	! -----
89	!
90	! ngrid Size of the horizontal grid.
91	! nlayer Number of vertical layers.
92	! nq Number of advected fields.
93	! nametrac Name of corresponding advected fields.
94	!
95	! firstcall True at the first call.
96	! lastcall True at the last call.
97	!
98	! pday Number of days counted from the North. Spring equinoxe.
99	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
100	! ptimestep timestep (s).
101	!
102	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa).
103	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
104	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2.s-2).
105	!
106	! pu(ngrid,nlayer) u, zonal component of the wind (ms-1).
107	! pv(ngrid,nlayer) v, meridional component of the wind (ms-1).
108	!
109	! pt(ngrid,nlayer) Temperature (K).
110	!
111	! pq(ngrid,nlayer,nq) Advected fields.
112	!
113	! pudyn(ngrid,nlayer) \
114	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
115	! ptdyn(ngrid,nlayer) / corresponding variables.
116	! pqdyn(ngrid,nlayer,nq) /
117	! flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
118	!
119	! OUTPUT
120	! ------
121	!
122	! pdu(ngrid,nlayer) \
123	! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding
124	! pdt(ngrid,nlayer) / variables due to physical processes.
125	! pdq(ngrid,nlayer) /
126	! pdpsrf(ngrid) /
127	!
128	!
129	! Authors
130	! -------
131	! Frederic Hourdin 15/10/93
132	! Francois Forget 1994
133	! Christophe Hourdin 02/1997
134	! Subroutine completely rewritten by F. Forget (01/2000)
135	! Water ice clouds: Franck Montmessin (update 06/2003)
136	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
137	! New correlated-k radiative scheme: R. Wordsworth (2009)
138	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
139	! Improved water cycle: R. Wordsworth / B. Charnay (2010)
140	! To F90: R. Wordsworth (2010)
141	! New turbulent diffusion scheme: J. Leconte (2012)
142	! Loops converted to F90 matrix format: J. Leconte (2012)
143	! No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
144	! Purge of the code : M. Turbet (2015)
145	! Fork for Titan : J. Vatant d'Ollone (2017)
146	! + clean of all too-generic (ocean, water, co2 ...) routines
147	! + Titan's chemistry
148	!============================================================================================
149
150
151	! 0. Declarations :
152	! ------------------
153
154	include "netcdf.inc"
155
156	! Arguments :
157	! -----------
158
159	! INPUTS:
160	! -------
161
162
163	integer,intent(in) :: ngrid ! Number of atmospheric columns.
164	integer,intent(in) :: nlayer ! Number of atmospheric layers.
165	integer,intent(in) :: nq ! Number of tracers.
166	character*20,intent(in) :: nametrac(nq) ! Names of the tracers taken from dynamics.
167
168	logical,intent(in) :: firstcall ! Signals first call to physics.
169	logical,intent(in) :: lastcall ! Signals last call to physics.
170
171	real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0.
172	real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
173	real,intent(in) :: ptimestep ! Physics timestep (s).
174	real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
175	real,intent(in) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
176	real,intent(in) :: pphi(ngrid,nlayer) ! Geopotential at mid-layer (m2s-2).
177	real,intent(in) :: pu(ngrid,nlayer) ! Zonal wind component (m/s).
178	real,intent(in) :: pv(ngrid,nlayer) ! Meridional wind component (m/s).
179	real,intent(in) :: pt(ngrid,nlayer) ! Temperature (K).
180	real,intent(in) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
181	real,intent(in) :: flxw(ngrid,nlayer) ! Vertical mass flux (ks/s) at lower boundary of layer
182
183	! OUTPUTS:
184	! --------
185
186	! Physical tendencies :
187
188	real,intent(out) :: pdu(ngrid,nlayer) ! Zonal wind tendencies (m/s/s).
189	real,intent(out) :: pdv(ngrid,nlayer) ! Meridional wind tendencies (m/s/s).
190	real,intent(out) :: pdt(ngrid,nlayer) ! Temperature tendencies (K/s).
191	real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
192	real,intent(out) :: pdpsrf(ngrid) ! Surface pressure tendency (Pa/s).
193
194	! Local saved variables:
195	! ----------------------
196
197	integer,save :: day_ini ! Initial date of the run (sol since Ls=0).
198	integer,save :: icount ! Counter of calls to physiq during the run.
199	!$OMP THREADPRIVATE(day_ini,icount)
200
201	real, dimension(:),allocatable,save :: tsurf ! Surface temperature (K).
202	real, dimension(:,:),allocatable,save :: tsoil ! Sub-surface temperatures (K).
203	real, dimension(:,:),allocatable,save :: albedo ! Surface Spectral albedo. By MT2015.
204	real, dimension(:),allocatable,save :: albedo_equivalent ! Spectral Mean albedo.
205
206	!$OMP THREADPRIVATE(tsurf,tsoil,albedo,albedo_equivalent)
207
208	real,dimension(:),allocatable,save :: albedo_bareground ! Bare Ground Albedo. By MT 2015.
209
210	!$OMP THREADPRIVATE(albedo_bareground)
211
212	real,dimension(:),allocatable,save :: emis ! Thermal IR surface emissivity.
213	real,dimension(:,:),allocatable,save :: dtrad ! Net atmospheric radiative heating rate (K.s-1).
214	real,dimension(:),allocatable,save :: fluxrad_sky ! Radiative flux from sky absorbed by surface (W.m-2).
215	real,dimension(:),allocatable,save :: fluxrad ! Net radiative surface flux (W.m-2).
216	real,dimension(:),allocatable,save :: capcal ! Surface heat capacity (J m-2 K-1).
217	real,dimension(:),allocatable,save :: fluxgrd ! Surface conduction flux (W.m-2).
218	real,dimension(:,:),allocatable,save :: qsurf ! Tracer on surface (e.g. kg.m-2).
219	real,dimension(:,:),allocatable,save :: q2 ! Turbulent Kinetic Energy.
220
221	!$OMP THREADPRIVATE(emis,dtrad,fluxrad_sky,fluxrad,capcal,fluxgrd,qsurf,q2)
222
223
224	! Local variables :
225	! -----------------
226
227	real zh(ngrid,nlayer) ! Potential temperature (K).
228	real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
229
230	integer l,ig,ierr,iq,nw,isoil
231
232	! FOR DIAGNOSTIC :
233
234	real,dimension(:),allocatable,save :: fluxsurf_lw ! Incident Long Wave (IR) surface flux (W.m-2).
235	real,dimension(:),allocatable,save :: fluxsurf_sw ! Incident Short Wave (stellar) surface flux (W.m-2).
236	real,dimension(:),allocatable,save :: fluxsurfabs_sw ! Absorbed Short Wave (stellar) flux by the surface (W.m-2).
237	real,dimension(:),allocatable,save :: fluxtop_lw ! Outgoing LW (IR) flux to space (W.m-2).
238	real,dimension(:),allocatable,save :: fluxabs_sw ! Absorbed SW (stellar) flux (W.m-2).
239	real,dimension(:),allocatable,save :: fluxtop_dn ! Incoming SW (stellar) radiation at the top of the atmosphere (W.m-2).
240	real,dimension(:),allocatable,save :: fluxdyn ! Horizontal heat transport by dynamics (W.m-2).
241	real,dimension(:,:),allocatable,save :: OLR_nu ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1)).
242	real,dimension(:,:),allocatable,save :: OSR_nu ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1)).
243	real,dimension(:,:),allocatable,save :: zdtlw ! LW heating tendencies (K/s).
244	real,dimension(:,:),allocatable,save :: zdtsw ! SW heating tendencies (K/s).
245	real,dimension(:),allocatable,save :: sensibFlux ! Turbulent flux given by the atmosphere to the surface (W.m-2).
246
247	!$OMP THREADPRIVATE(fluxsurf_lw,fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,&
248
249	!$OMP zdtlw,zdtsw,sensibFlux)
250
251	real zls ! Solar longitude (radians).
252	real zlss ! Sub solar point longitude (radians).
253	real zday ! Date (time since Ls=0, calculated in sols).
254	real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers.
255	real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries.
256
257	! TENDENCIES due to various processes :
258
259	! For Surface Temperature : (K/s)
260	real zdtsurf(ngrid) ! Cumulated tendencies.
261	real zdtsurfmr(ngrid) ! Mass_redistribution routine.
262	real zdtsdif(ngrid) ! Turbdiff/vdifc routines.
263
264	! For Atmospheric Temperatures : (K/s)
265	real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
266	real zdtmr(ngrid,nlayer) ! Mass_redistribution routine.
267	real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
268
269	! For Surface Tracers : (kg/m2/s)
270	real dqsurf(ngrid,nq) ! Cumulated tendencies.
271	real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines.
272	real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine.
273
274	! For Tracers : (kg/kg_of_air/s)
275	real zdqadj(ngrid,nlayer,nq) ! Convadj routine.
276	real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines.
277	real zdqevap(ngrid,nlayer) ! Turbdiff routine.
278	real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine.
279
280	real zdqchi(ngrid,nlayer,nq) ! Chemical tendency ( chemistry routine ).
281	real zdqcond(ngrid,nlayer,nq) ! Condensation tendency ( chemistry routine ).
282
283	real zdqmufi(ngrid,nlayer,nq) ! Microphysical tendency.
284
285	! For Winds : (m/s/s)
286	real zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! Convadj routine.
287	real zdumr(ngrid,nlayer),zdvmr(ngrid,nlayer) ! Mass_redistribution routine.
288	real zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
289	real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
290	real zdhadj(ngrid,nlayer) ! Convadj routine.
291
292	! For Pressure and Mass :
293	real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
294	real zdmassmr_col(ngrid) ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
295	real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s).
296
297
298
299	! Local variables for LOCAL CALCULATIONS:
300	! ---------------------------------------
301	real zflubid(ngrid)
302	real zplanck(ngrid),zpopsk(ngrid,nlayer)
303	real ztim1,ztim2,ztim3, z1,z2
304	real ztime_fin
305	real zdh(ngrid,nlayer)
306	real gmplanet
307	real taux(ngrid),tauy(ngrid)
308
309
310	! local variables for DIAGNOSTICS : (diagfi & stat)
311	! -------------------------------------------------
312	real ps(ngrid) ! Surface Pressure.
313	real zt(ngrid,nlayer) ! Atmospheric Temperature.
314	real zu(ngrid,nlayer),zv(ngrid,nlayer) ! Zonal and Meridional Winds.
315	real zq(ngrid,nlayer,nq) ! Atmospheric Tracers.
316	real zdtadj(ngrid,nlayer) ! Convadj Diagnostic.
317	real zdtdyn(ngrid,nlayer) ! Dynamical Heating (K/s).
318	real zdudyn(ngrid,nlayer) ! Dynamical Zonal Wind tendency (m.s-2).
319	real,allocatable,dimension(:,:),save :: ztprevious ! Previous loop Atmospheric Temperature (K) ! Useful for Dynamical Heating calculation.
320	real,allocatable,dimension(:,:),save :: zuprevious ! Previous loop Zonal Wind (m.s-1) ! Useful for Zonal Wind tendency calculation.
321	!$OMP THREADPRIVATE(ztprevious,zuprevious)
322
323	real vmr(ngrid,nlayer) ! volume mixing ratio
324	real time_phys
325
326	real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
327
328	! to test energy conservation (RW+JL)
329	real mass(ngrid,nlayer),massarea(ngrid,nlayer)
330	real dEtot, dEtots, AtmToSurf_TurbFlux
331	real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
332	!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
333	real dEzRadsw(ngrid,nlayer),dEzRadlw(ngrid,nlayer),dEzdiff(ngrid,nlayer)
334	real dEdiffs(ngrid),dEdiff(ngrid)
335
336	!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
337
338	real dItot, dItot_tmp, dVtot, dVtot_tmp
339
340	real dWtot, dWtot_tmp, dWtots, dWtots_tmp
341
342
343	! For Clear Sky Case.
344	real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxsurfabs_sw1(ngrid) ! For SW/LW flux.
345	real fluxtop_lw1(ngrid), fluxabs_sw1(ngrid) ! For SW/LW flux.
346	real albedo_equivalent1(ngrid) ! For Equivalent albedo calculation.
347	real tf, ntf
348
349	real,allocatable,dimension(:,:),save :: qsurf_hist
350	!$OMP THREADPRIVATE(qsurf_hist)
351
352	! Local variables for Titan chemistry and microphysics (JVO 2017)
353	! ----------------------------------------------------------------------------
354
355	character*10,dimension(:),allocatable,save :: nomqy
356
357	!$OMP THREADPRIVATE(nomqy)
358
359	real ctimestep ! Chemistry timestep (s)
360
361	! Grandeurs en moyennes zonales ------------------------
362	real temp_eq(nlayer), press_eq(nlayer)
363	real zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
364	real ztemp(ngrid,nlayer)
365
366	real , allocatable, dimension(:,:,:),save :: ychim
367
368	! 2D vmr tendencies ( chemistry and condensation ) for chem. tracers
369	real,dimension(:,:,:),allocatable,save :: dycchi ! Saved since chemistry is not called every step
370	real dyccond(ngrid,nlayer,nq)
371
372	real,dimension(:,:),allocatable,save :: qysat
373
374	!$OMP THREADPRIVATE(ychim,dycchi,qysat)
375
376	real,dimension(:),allocatable,save :: tankCH4 ! Depth of surface methane tank
377	!$OMP THREADPRIVATE(tankCH4)
378
379	!-----------------------------------------------------------------------------
380	! Interface to calmufi
381	! --> needed in order to pass assumed-shape arrays. Otherwise we must put calmufi in a module
382	! (to have an explicit generated by the compiler).
383	! Or one can put calmufi in MMP_GCM module (in muphytitan).
384	INTERFACE
385	SUBROUTINE calmufi(plev, zlev, play, zlay, temp, pq, zdq)
386	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: plev !! Pressure levels (Pa).
387	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlev !! Altitude levels (m).
388	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: play !! Pressure layers (Pa).
389	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlay !! Altitude at the center of each layer (m).
390	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp !! Temperature at the center of each layer (K).
391	REAL(kind=8), DIMENSION(:,:,:), INTENT(IN) :: pq !! Tracers ($kg.kg^{-1}}$).
392	REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq !! Tendency for tracers ($kg.kg^{-1}}$).
393	END SUBROUTINE calmufi
394	END INTERFACE
395
396	!==================================================================================================
397
398	! -----------------
399	! I. INITIALISATION
400	! -----------------
401
402	! --------------------------------
403	! I.1 First Call Initialisation.
404	! --------------------------------
405	if (firstcall) then
406	! initracer:
407	! initialize nmicro,
408	call initracer2(nq,nametrac)
409
410	! ----------------------------------------------------------------------------
411
412	! Allocate saved arrays.
413	ALLOCATE(tsurf(ngrid))
414	ALLOCATE(tsoil(ngrid,nsoilmx))
415	ALLOCATE(albedo(ngrid,L_NSPECTV))
416	ALLOCATE(albedo_equivalent(ngrid))
417	ALLOCATE(albedo_bareground(ngrid))
418	ALLOCATE(emis(ngrid))
419	ALLOCATE(dtrad(ngrid,nlayer))
420	ALLOCATE(fluxrad_sky(ngrid))
421	ALLOCATE(fluxrad(ngrid))
422	ALLOCATE(capcal(ngrid))
423	ALLOCATE(fluxgrd(ngrid))
424	ALLOCATE(qsurf(ngrid,nq))
425	ALLOCATE(q2(ngrid,nlayer+1))
426	ALLOCATE(tankCH4(ngrid))
427	ALLOCATE(ztprevious(ngrid,nlayer))
428	ALLOCATE(zuprevious(ngrid,nlayer))
429	ALLOCATE(qsurf_hist(ngrid,nq))
430	ALLOCATE(fluxsurf_lw(ngrid))
431	ALLOCATE(fluxsurf_sw(ngrid))
432	ALLOCATE(fluxsurfabs_sw(ngrid))
433	ALLOCATE(fluxtop_lw(ngrid))
434	ALLOCATE(fluxabs_sw(ngrid))
435	ALLOCATE(fluxtop_dn(ngrid))
436	ALLOCATE(fluxdyn(ngrid))
437	ALLOCATE(OLR_nu(ngrid,L_NSPECTI))
438	ALLOCATE(OSR_nu(ngrid,L_NSPECTV))
439	ALLOCATE(sensibFlux(ngrid))
440	ALLOCATE(zdtlw(ngrid,nlayer))
441	ALLOCATE(zdtsw(ngrid,nlayer))
442	ALLOCATE(dycchi(ngrid,nlayer,nq-nmicro))
443	ALLOCATE(qysat(nlayer,nq-nmicro))
444	ALLOCATE(nomqy(nq-nmicro+1)) ! +1 because of hv
445
446	! This is defined in comsaison_h
447	ALLOCATE(mu0(ngrid))
448	ALLOCATE(fract(ngrid))
449	! This is defined in radcommon_h
450	ALLOCATE(glat(ngrid))
451
452
453	! Variables set to 0
454	! ~~~~~~~~~~~~~~~~~~
455	dtrad(:,:) = 0.0
456	fluxrad(:) = 0.0
457	zdtsw(:,:) = 0.0
458	zdtlw(:,:) = 0.0
459
460	! Initialize names, timestep and saturation profiles for chemistry
461	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
462
463	if ( callchim .and. (nq.gt.nmicro) ) then
464
465	ctimestep = ptimestep*REAL(ichim)
466
467	do iq=nmicro+1,nq
468	nomqy(iq-nmicro) = nametrac(iq)
469	enddo
470
471	nomqy(nq-nmicro+1) = "HV"
472
473	! qysat is taken at the equator ( small variations of t,p)
474	temp_eq(:) = tmoy(:)
475	press_eq(:) = playmoy(:)/100. ! en mbar
476
477	call inicondens(nq-nmicro,press_eq,temp_eq,nomqy,qysat)
478
479	zdqchi(:,:,:) = 0.0
480	zdqcond(:,:,:) = 0.0
481
482	endif
483
484	! Initialize microphysics.
485	! ~~~~~~~~~~~~~~~~~~~~~~~~
486
487	IF ( callmufi ) THEN
488
489	call inimufi(nq,ptimestep)
490
491	ENDIF
492
493
494	! Read 'startfi.nc' file.
495	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
496	call phyetat0(startphy_file,ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
497	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,tankCH4)
498	if (.not.startphy_file) then
499	! additionnal "academic" initialization of physics
500	if (is_master) write(,) "Physiq: initializing tsurf(:) to pt(:,1) !!"
501	tsurf(:)=pt(:,1)
502	if (is_master) write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
503	do isoil=1,nsoilmx
504	tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
505	enddo
506	if (is_master) write(,) "Physiq: initializing day_ini to pdat !"
507	day_ini=pday
508	endif
509
510	if (pday.ne.day_ini) then
511	write(,) "ERROR in physiq.F90:"
512	write(,) "bad synchronization between physics and dynamics"
513	write(,) "dynamics day: ",pday
514	write(,) "physics day: ",day_ini
515	stop
516	endif
517
518	write (,) 'In physiq day_ini =', day_ini
519
520	! Initialize albedo calculation.
521	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
522	albedo(:,:)=0.0
523	albedo_bareground(:)=0.0
524	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground)
525
526	! Initialize orbital calculation.
527	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
528	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
529
530	if(tlocked)then
531	print*,'Planet is tidally locked at resonance n=',nres
532	print*,'Make sure you have the right rotation rate!!!'
533	endif
534
535
536	! Initialize soil.
537	! ~~~~~~~~~~~~~~~~
538	if (callsoil) then
539
540	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
541	ptimestep,tsurf,tsoil,capcal,fluxgrd)
542
543	else ! else of 'callsoil'.
544
545	print*,'WARNING! Thermal conduction in the soil turned off'
546	capcal(:)=1.e6
547	fluxgrd(:)=intheat
548	print*,'Flux from ground = ',intheat,' W m^-2'
549
550	endif ! end of 'callsoil'.
551
552	icount=1
553
554
555	! Initialize surface history variable.
556	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
557	qsurf_hist(:,:)=qsurf(:,:)
558
559	! Initialize variable for dynamical heating and zonal wind tendency diagnostic
560	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
561	ztprevious(:,:)=pt(:,:)
562	zuprevious(:,:)=pu(:,:)
563
564
565	if(meanOLR)then
566	call system('rm -f rad_bal.out') ! to record global radiative balance.
567	call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
568	call system('rm -f h2o_bal.out') ! to record global hydrological balance.
569	endif
570
571	if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
572	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
573	ptimestep,pday+nday,time_phys,cell_area, &
574	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
575	endif
576
577	! Sanity check for microphysics
578	if ( ((.not.moyzon_mu).and.(callmufi)) ) then
579	print *, "moyzon_mu=",moyzon_mu," and callmufi=",callmufi
580	print *, "Please activate zonal mean to run microphysics (for now) !"
581	stop
582	endif
583
584	! Sanity check for chemistry
585	if ( (.not.moyzon_ch) .and. (callchim) ) then
586	print *, "moyzon_ch=",moyzon_ch," and callchim=",callchim
587	print *, "Please activate zonal mean to run chemistry !"
588	stop
589	endif
590
591
592	! XIOS outputs
593	#ifdef CPP_XIOS
594
595	write(,) "physiq: call initialize_xios_output"
596	call initialize_xios_output(pday,ptime,ptimestep,daysec, &
597	presnivs,pseudoalt)
598	#endif
599	endif ! end of 'firstcall'
600
601	! ------------------------------------------------------
602	! I.2 Initializations done at every physical timestep:
603	! ------------------------------------------------------
604
605	#ifdef CPP_XIOS
606	! update XIOS time/calendar
607	call update_xios_timestep
608	#endif
609
610	! Initialize various variables
611	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
612
613	pdt(1:ngrid,1:nlayer) = 0.0
614	zdtsurf(1:ngrid) = 0.0
615	pdq(1:ngrid,1:nlayer,1:nq) = 0.0
616	dqsurf(1:ngrid,1:nq)= 0.0
617	pdu(1:ngrid,1:nlayer) = 0.0
618	pdv(1:ngrid,1:nlayer) = 0.0
619	pdpsrf(1:ngrid) = 0.0
620	zflubid(1:ngrid) = 0.0
621	taux(1:ngrid) = 0.0
622	tauy(1:ngrid) = 0.0
623
624	zday=pday+ptime ! Compute time, in sols (and fraction thereof).
625
626	! Compute Stellar Longitude (Ls), and orbital parameters.
627	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
628	if (season) then
629	call stellarlong(zday,zls)
630	else
631	call stellarlong(float(day_ini),zls)
632	end if
633
634	call orbite(zls,dist_star,declin,right_ascen)
635
636	if (tlocked) then
637	zlss=Mod(-(2.pi(zday/year_day)nres - right_ascen),2.pi)
638	elseif (diurnal) then
639	zlss=-2.pi(zday-.5)
640	else if(diurnal .eqv. .false.) then
641	zlss=9999.
642	endif
643
644
645	! Compute variations of g with latitude (oblate case).
646	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
647	if (oblate .eqv. .false.) then
648	glat(:) = g
649	else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
650	print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
651	call abort
652	else
653	gmplanet = MassPlanetgrav1e24
654	do ig=1,ngrid
655	glat(ig)= gmplanet/(Rmean*2) (1.D0 + 0.75 J2 - 2.0flatten/3. + (2.flatten - 15./4. J2) * cos(2. * (pi/2. - latitude(ig))))
656	end do
657	endif
658
659	! Compute geopotential between layers.
660	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
661	zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
662	do l=1,nlayer
663	zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
664	enddo
665
666	zzlev(1:ngrid,1)=0.
667	zzlev(1:ngrid,nlayer+1)=1.e7 ! Dummy top of last layer above 10000 km...
668
669	do l=2,nlayer
670	do ig=1,ngrid
671	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
672	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
673	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
674	enddo
675	enddo
676
677	! -------------------------------Taken from old Titan --------------------------
678	! zonal averages needed
679	if (moyzon_ch .or. moyzon_mu) then
680
681	zzlaybar(1,:)=(zphibar(1,:)+zphisbar(1))/g
682	! SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
683	! zzlaybar(1,:)=RGRARA/(RG*RA-(zphibar(1,:)+zphisbar(1)))-RA
684	zzlevbar(1,1)=zphisbar(1)/g
685	DO l=2,nlayer
686	z1=(zplaybar(1,l-1)+zplevbar(1,l))/(zplevbar(1,l-1)-zplevbar(1,l))
687	z2=(zplevbar(1,l) +zplaybar(1,l))/(zplevbar(1,l) -zplaybar(1,l))
688	zzlevbar(1,l)=(z1zzlaybar(1,l-1)+z2zzlaybar(1,l))/(z1+z2)
689	ENDDO
690	zzlevbar(1,nlayer+1)=zzlaybar(1,nlayer)+(zzlaybar(1,nlayer)-zzlevbar(1,nlayer))
691
692	DO ig=2,ngrid
693	if (latitude(ig).ne.latitude(ig-1)) then
694	DO l=1,nlayer
695	zzlaybar(ig,l)=(zphibar(ig,l)+zphisbar(ig))/g
696	! SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
697	!zzlaybar(ig,l)=RGRARA/(RG*RA-(zphibar(ig,l)+zphisbar(ig)))-RA
698	ENDDO
699	zzlevbar(ig,1)=zphisbar(ig)/g
700	DO l=2,nlayer
701	z1=(zplaybar(ig,l-1)+zplevbar(ig,l))/ (zplevbar(ig,l-1)-zplevbar(ig,l))
702	z2=(zplevbar(ig,l) +zplaybar(ig,l))/(zplevbar(ig,l) -zplaybar(ig,l))
703	zzlevbar(ig,l)=(z1zzlaybar(ig,l-1)+z2zzlaybar(ig,l))/(z1+z2)
704	ENDDO
705	zzlevbar(ig,nlayer+1)=zzlaybar(ig,nlayer)+(zzlaybar(ig,nlayer)-zzlevbar(ig,nlayer))
706	else
707	zzlaybar(ig,:)=zzlaybar(ig-1,:)
708	zzlevbar(ig,:)=zzlevbar(ig-1,:)
709	endif
710	ENDDO
711
712	endif ! moyzon
713
714	! -------------------------------------------------------------------------------------
715	! Compute potential temperature
716	! Note : Potential temperature calculation may not be the same in physiq and dynamic...
717	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
718	do l=1,nlayer
719	do ig=1,ngrid
720	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
721	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
722	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
723	massarea(ig,l)=mass(ig,l)*cell_area(ig)
724	enddo
725	enddo
726
727	! Compute vertical velocity (m/s) from vertical mass flux
728	! w = F / (rhoarea) and rho = P/(rT)
729	! But first linearly interpolate mass flux to mid-layers
730	do l=1,nlayer-1
731	pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
732	enddo
733	pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
734	do l=1,nlayer
735	pw(1:ngrid,l)=(pw(1:ngrid,l)rpt(1:ngrid,l)) / &
736	(pplay(1:ngrid,l)*cell_area(1:ngrid))
737	enddo
738
739	!---------------------------------
740	! II. Compute radiative tendencies
741	!---------------------------------
742
743	if (callrad) then
744	if( mod(icount-1,iradia).eq.0.or.lastcall) then
745
746	! Compute local stellar zenith angles
747	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
748	if (tlocked) then
749	! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb
750	ztim1=SIN(declin)
751	ztim2=COS(declin)*COS(zlss)
752	ztim3=COS(declin)*SIN(zlss)
753
754	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
755	ztim1,ztim2,ztim3,mu0,fract, flatten)
756
757	elseif (diurnal) then
758	ztim1=SIN(declin)
759	ztim2=COS(declin)COS(2.pi*(zday-.5))
760	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
761
762	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
763	ztim1,ztim2,ztim3,mu0,fract, flatten)
764	else if(diurnal .eqv. .false.) then
765
766	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad,flatten)
767	! WARNING: this function appears not to work in 1D
768
769	endif
770
771	! Eclipse incoming sunlight (e.g. Saturn ring shadowing).
772	if(rings_shadow) then
773	call call_rings(ngrid, ptime, pday, diurnal)
774	endif
775
776
777	if (corrk) then
778
779	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
780	! II.a Call correlated-k radiative transfer scheme
781	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
782
783	call call_profilgases(nlayer)
784
785	! standard callcorrk
786	call callcorrk(ngrid,nlayer,pq,nq,qsurf, &
787	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
788	tsurf,fract,dist_star, &
789	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, &
790	fluxsurfabs_sw,fluxtop_lw, &
791	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, &
792	firstcall,lastcall)
793
794	! Radiative flux from the sky absorbed by the surface (W.m-2).
795	GSR=0.0
796	fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurfabs_sw(1:ngrid)
797
798	!if(noradsurf)then ! no lower surface; SW flux just disappears
799	! GSR = SUM(fluxsurf_sw(1:ngrid)*cell_area(1:ngrid))/totarea
800	! fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
801	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
802	!endif
803
804	! Net atmospheric radiative heating rate (K.s-1)
805	dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
806
807	elseif(newtonian)then
808
809	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
810	! II.b Call Newtonian cooling scheme
811	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
812	call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
813
814	zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep
815	! e.g. surface becomes proxy for 1st atmospheric layer ?
816
817	else
818
819	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
820	! II.c Atmosphere has no radiative effect
821	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
822	fluxtop_dn(1:ngrid) = fract(1:ngrid)mu0(1:ngrid)Fat1AU/dist_star**2
823	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
824	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
825	endif
826	fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
827	print*,'------------WARNING---WARNING------------' ! by MT2015.
828	print*,'You are in corrk=false mode, '
829	print*,'and the surface albedo is taken equal to the first visible spectral value'
830
831	fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
832	fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid)
833	fluxtop_lw(1:ngrid) = emis(1:ngrid)sigmatsurf(1:ngrid)**4
834
835	dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating
836
837	endif ! end of corrk
838
839	endif ! of if(mod(icount-1,iradia).eq.0)
840
841
842	! Transformation of the radiative tendencies
843	! ------------------------------------------
844	zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
845	zplanck(1:ngrid)=emis(1:ngrid)sigmazplanck(1:ngrid)*zplanck(1:ngrid)
846	fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
847	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
848
849	! Test of energy conservation
850	!----------------------------
851	if(enertest)then
852	call planetwide_sumval(cppmassarea(:,:)zdtsw(:,:)/totarea_planet,dEtotSW)
853	call planetwide_sumval(cppmassarea(:,:)zdtlw(:,:)/totarea_planet,dEtotLW)
854	!call planetwide_sumval(fluxsurf_sw(:)(1.-albedo_equivalent(:))cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
855	call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
856	call planetwide_sumval((fluxsurf_lw(:)emis(:)-zplanck(:))cell_area(:)/totarea_planet,dEtotsLW)
857	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
858	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
859	if (is_master) then
860	print*,'---------------------------------------------------------------'
861	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
862	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
863	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
864	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
865	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
866	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
867	endif
868	endif ! end of 'enertest'
869
870	endif ! of if (callrad)
871
872
873
874	! --------------------------------------------
875	! III. Vertical diffusion (turbulent mixing) :
876	! --------------------------------------------
877
878	if (calldifv) then
879
880	zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
881
882	! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
883	if (UseTurbDiff) then
884
885	call turbdiff(ngrid,nlayer,nq, &
886	ptimestep,capcal,lwrite, &
887	pplay,pplev,zzlay,zzlev,z0, &
888	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
889	pdt,pdq,zflubid, &
890	zdudif,zdvdif,zdtdif,zdtsdif, &
891	sensibFlux,q2,zdqdif,zdqsdif, &
892	taux,tauy,lastcall)
893
894	else
895
896	zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
897
898	call vdifc(ngrid,nlayer,nq,zpopsk, &
899	ptimestep,capcal,lwrite, &
900	pplay,pplev,zzlay,zzlev,z0, &
901	pu,pv,zh,pq,tsurf,emis,qsurf, &
902	zdh,pdq,zflubid, &
903	zdudif,zdvdif,zdhdif,zdtsdif, &
904	sensibFlux,q2,zdqdif,zdqsdif, &
905	taux,tauy,lastcall)
906
907	zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
908	zdqevap(1:ngrid,1:nlayer)=0.
909
910	end if !end of 'UseTurbDiff'
911
912
913	pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
914	pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
915	pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
916	zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
917
918	if (tracer) then
919	pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
920	dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
921	end if ! of if (tracer)
922
923
924	! test energy conservation
925	!-------------------------
926	if(enertest)then
927
928	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
929	do ig = 1, ngrid
930	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
931	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
932	enddo
933
934	call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
935	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
936	call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
937	call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
938
939	if (is_master) then
940
941	if (UseTurbDiff) then
942	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
943	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
944	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
945	else
946	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
947	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
948	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
949	end if
950	endif ! end of 'is_master'
951
952	! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
953	endif ! end of 'enertest'
954
955	else ! calldifv
956
957	if(.not.newtonian)then
958
959	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
960
961	endif
962
963	endif ! end of 'calldifv'
964
965
966	!----------------------------------
967	! IV. Dry convective adjustment :
968	!----------------------------------
969
970	if(calladj) then
971
972	zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
973	zduadj(1:ngrid,1:nlayer)=0.0
974	zdvadj(1:ngrid,1:nlayer)=0.0
975	zdhadj(1:ngrid,1:nlayer)=0.0
976
977
978	call convadj(ngrid,nlayer,nq,ptimestep, &
979	pplay,pplev,zpopsk, &
980	pu,pv,zh,pq, &
981	pdu,pdv,zdh,pdq, &
982	zduadj,zdvadj,zdhadj, &
983	zdqadj)
984
985	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
986	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer)
987	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
988	zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
989
990	if(tracer) then
991	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
992	end if
993
994	! Test energy conservation
995	if(enertest)then
996	call planetwide_sumval(cppmassarea(:,:)zdtadj(:,:)/totarea_planet,dEtot)
997	if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
998	endif
999
1000
1001	endif ! end of 'calladj'
1002
1003
1004	!---------------------------------------------
1005	! V. Specific parameterizations for tracers
1006	!---------------------------------------------
1007
1008	if (tracer) then
1009
1010	! -------------------------
1011	! V.1. Chemistry
1012	! -------------------------
1013
1014	if (callchim) then
1015
1016	! Using zonal mean for calchim
1017	zplev(:,:) = zplevbar(:,:)
1018	zplay(:,:) = zplaybar(:,:)
1019	zzlev(:,:) = zzlevbar(:,:)
1020	zzlay(:,:) = zzlaybar(:,:)
1021	ztemp(:,:) = ztfibar(:,:)
1022
1023	allocate(ychim(ngrid,nlayer,nq-nmicro))
1024
1025	if (nq.gt.nmicro) then
1026	do iq = nmicro+1,nq
1027	ychim(:,:,iq-nmicro) = pq(:,:,iq) * rat_mmol(iq) ! convert to molar fraction
1028	enddo
1029	endif
1030
1031	! Condensation tendency after the transport
1032	do iq=nmicro+1,nq
1033	do l=1,nlayer
1034	do ig=1,ngrid
1035	if ( ychim(ig,l,iq-nmicro).gt.qysat(l,iq-nmicro) ) then
1036	dyccond(ig,l,iq)= ( -ychim(ig,l,iq-nmicro)+qysat(l,iq-nmicro) ) / ptimestep
1037	endif
1038	enddo
1039	enddo
1040	enddo
1041
1042	if ( callclouds ) dyccond(:,:,ices_indx) = 0.0 ! Condensation will be calculated in the cloud microphysics
1043
1044	if( mod(icount-1,ichim).eq.0. ) then
1045
1046	print *, "We enter in the chemistry ..."
1047	call calchim(ngrid,nq-nmicro,ychim,nomqy,declin,zls,ctimestep, &
1048	ztemp,zplay,zplev,zzlay,zzlev,dycchi,nlayer+70)
1049
1050	! JVO 2017 : NLEV = nlayer+70, en accord avec le C. Quid si nlay=/ 55 ?
1051
1052	endif
1053
1054	! TEMPORARY COMMENT
1055	! These conversions as well as 2D/1D should be put in phytrac
1056	! ( GCM-friendly tracers and tendencies -> format for photochem and microphys )
1057
1058	if (nq.gt.nmicro) then
1059	! We convert tendencies to mass mixing ratio
1060	do iq=nmicro+1,nq
1061	zdqchi(:,:,iq) = dycchi(:,:,iq-nmicro) / rat_mmol(iq)
1062	zdqcond(:,:,iq) = dyccond(:,:,iq) / rat_mmol(iq)
1063	enddo
1064
1065	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + &
1066	zdqchi(1:ngrid,1:nlayer,1:nq) + zdqcond(1:ngrid,1:nlayer,1:nq)
1067
1068	endif
1069
1070	endif ! end of 'callchim'
1071
1072	! -------------------
1073	! V.2 Microphysics
1074	! -------------------
1075
1076	if (callmufi) then
1077
1078	! Using zonal mean for microphysics
1079	zplev(:,:) = zplevbar(:,:)
1080	zplay(:,:) = zplaybar(:,:)
1081	zzlev(:,:) = zzlevbar(:,:)
1082	zzlay(:,:) = zzlaybar(:,:)
1083	ztemp(:,:) = ztfibar(:,:)
1084
1085	! Inside this routine we will split 2D->1D, intensive->extensive and separate different types of tracers
1086	! Should be put in phytrac
1087
1088	call calmufi(zplev,zzlev,zplay,zzlay,ztemp,pq,zdqmufi)
1089
1090	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmufi(1:ngrid,1:nlayer,1:nq)
1091
1092	endif ! end of 'callmufi'
1093
1094	! ---------------
1095	! V.3 Updates
1096	! ---------------
1097
1098	! Updating Atmospheric Mass and Tracers budgets.
1099	if(mass_redistrib) then
1100
1101	zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * zdqevap(1:ngrid,1:nlayer)
1102
1103	do ig = 1, ngrid
1104	zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
1105	enddo
1106
1107	call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
1108	call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
1109	call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
1110
1111	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
1112	capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
1113	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
1114	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
1115
1116	pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
1117	dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
1118	pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
1119	pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
1120	pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
1121	pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
1122	zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
1123
1124	endif
1125
1126	! -----------------------------
1127	! V.4. Surface Tracer Update
1128	! -----------------------------
1129
1130	qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
1131
1132	! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
1133	! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
1134	qsurf_hist(:,:) = qsurf(:,:)
1135
1136	endif! end of if 'tracer'
1137
1138
1139	!------------------------------------------------
1140	! VI. Surface and sub-surface soil temperature
1141	!------------------------------------------------
1142
1143
1144	! Increment surface temperature
1145
1146	tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
1147
1148	! Compute soil temperatures and subsurface heat flux.
1149	if (callsoil) then
1150	call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, &
1151	ptimestep,tsurf,tsoil,capcal,fluxgrd)
1152	endif
1153
1154
1155	! Test energy conservation
1156	if(enertest)then
1157	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf(:)/totarea_planet,dEtots)
1158	if (is_master) print*,'Surface energy change =',dEtots,' W m-2'
1159	endif
1160
1161
1162	!---------------------------------------------------
1163	! VII. Perform diagnostics and write output files
1164	!---------------------------------------------------
1165
1166	! Note : For output only: the actual model integration is performed in the dynamics.
1167
1168
1169
1170	! Temperature, zonal and meridional winds.
1171	zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
1172	zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
1173	zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
1174
1175	! Diagnostic.
1176	zdtdyn(1:ngrid,1:nlayer) = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep
1177	ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
1178
1179	zdudyn(1:ngrid,1:nlayer) = (pu(1:ngrid,1:nlayer)-zuprevious(1:ngrid,1:nlayer)) / ptimestep
1180	zuprevious(1:ngrid,1:nlayer) = zu(1:ngrid,1:nlayer)
1181
1182	if(firstcall)then
1183	zdtdyn(1:ngrid,1:nlayer)=0.0
1184	zdudyn(1:ngrid,1:nlayer)=0.0
1185	endif
1186
1187	! Dynamical heating diagnostic.
1188	do ig=1,ngrid
1189	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp
1190	enddo
1191
1192	! Tracers.
1193	zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
1194
1195	! Surface pressure.
1196	ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
1197
1198
1199
1200	! Surface and soil temperature information
1201	call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
1202	call planetwide_minval(tsurf(:),Ts2)
1203	call planetwide_maxval(tsurf(:),Ts3)
1204	if(callsoil)then
1205	TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
1206	if (is_master) then
1207	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
1208	print*,Ts1,Ts2,Ts3,TsS
1209	end if
1210	else
1211	if (is_master) then
1212	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
1213	print*,Ts1,Ts2,Ts3
1214	endif
1215	end if
1216
1217
1218	! Check the energy balance of the simulation during the run
1219	if(corrk)then
1220
1221	call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
1222	call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
1223	call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
1224	call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
1225	call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
1226	do ig=1,ngrid
1227	if(fluxtop_dn(ig).lt.0.0)then
1228	print*,'fluxtop_dn has gone crazy'
1229	print*,'fluxtop_dn=',fluxtop_dn(ig)
1230	print*,'temp= ',pt(ig,:)
1231	print*,'pplay= ',pplay(ig,:)
1232	call abort
1233	endif
1234	end do
1235
1236	if(ngrid.eq.1)then
1237	DYN=0.0
1238	endif
1239
1240	if (is_master) then
1241	print*,' ISR ASR OLR GND DYN [W m^-2]'
1242	print*, ISR,ASR,OLR,GND,DYN
1243	endif
1244
1245	if(enertest .and. is_master)then
1246	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
1247	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
1248	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
1249	endif
1250
1251	if(meanOLR .and. is_master)then
1252	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
1253	! to record global radiative balance
1254	open(92,file="rad_bal.out",form='formatted',position='append')
1255	write(92,*) zday,ISR,ASR,OLR
1256	close(92)
1257	open(93,file="tem_bal.out",form='formatted',position='append')
1258	if(callsoil)then
1259	write(93,*) zday,Ts1,Ts2,Ts3,TsS
1260	else
1261	write(93,*) zday,Ts1,Ts2,Ts3
1262	endif
1263	close(93)
1264	endif
1265	endif
1266
1267	endif ! end of 'corrk'
1268
1269
1270	! Diagnostic to test radiative-convective timescales in code.
1271	if(testradtimes)then
1272	open(38,file="tau_phys.out",form='formatted',position='append')
1273	ig=1
1274	do l=1,nlayer
1275	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1276	enddo
1277	close(38)
1278	print*,'As testradtimes enabled,'
1279	print*,'exiting physics on first call'
1280	call abort
1281	endif
1282
1283
1284	if (is_master) print,'--> Ls =',zls180./pi
1285
1286
1287	!----------------------------------------------------------------------
1288	! Writing NetCDF file "RESTARTFI" at the end of the run
1289	!----------------------------------------------------------------------
1290
1291	! Note: 'restartfi' is stored just before dynamics are stored
1292	! in 'restart'. Between now and the writting of 'restart',
1293	! there will have been the itau=itau+1 instruction and
1294	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1295	! thus we store for time=time+dtvr
1296
1297
1298
1299	if(lastcall) then
1300	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1301
1302	if (ngrid.ne.1) then
1303	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1304
1305	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
1306	ptimestep,ztime_fin, &
1307	tsurf,tsoil,emis,q2,qsurf_hist,tankCH4)
1308	endif
1309
1310	endif ! end of 'lastcall'
1311
1312
1313	!-----------------------------------
1314	! Saving statistics :
1315	!-----------------------------------
1316
1317	! Note :("stats" stores and accumulates 8 key variables in file "stats.nc"
1318	! which can later be used to make the statistic files of the run:
1319	! "stats") only possible in 3D runs !!!
1320
1321
1322	if (callstats) then
1323
1324	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1325	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1326	call wstats(ngrid,"fluxsurf_lw", &
1327	"Thermal IR radiative flux to surface","W.m-2",2, &
1328	fluxsurf_lw)
1329	call wstats(ngrid,"fluxtop_lw", &
1330	"Thermal IR radiative flux to space","W.m-2",2, &
1331	fluxtop_lw)
1332
1333	! call wstats(ngrid,"fluxsurf_sw", &
1334	! "Solar radiative flux to surface","W.m-2",2, &
1335	! fluxsurf_sw_tot)
1336	! call wstats(ngrid,"fluxtop_sw", &
1337	! "Solar radiative flux to space","W.m-2",2, &
1338	! fluxtop_sw_tot)
1339
1340
1341	call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
1342	call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
1343	call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
1344	!call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
1345	!call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
1346	call wstats(ngrid,"p","Pressure","Pa",3,pplay)
1347	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1348	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1349	call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
1350	call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
1351	call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
1352
1353	if (tracer) then
1354	do iq=1,nq
1355	call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
1356	call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', &
1357	'kg m^-2',2,qsurf(1,iq) )
1358
1359	! call wstats(ngrid,trim(noms(iq))//'_reff', &
1360	! trim(noms(iq))//'_reff', &
1361	! 'm',3,reffrad(1,1,iq))
1362
1363	end do
1364
1365	endif ! end of 'tracer'
1366
1367	if(lastcall) then
1368	write (,) "Writing stats..."
1369	call mkstats(ierr)
1370	endif
1371
1372	endif ! end of 'callstats'
1373
1374
1375	!-----------------------------------------------------------------------------------------------------
1376	! OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
1377	!
1378	! Note 1 : output with period "ecritphy", set in "run.def"
1379	!
1380	! Note 2 : writediagfi can also be called from any other subroutine for any variable,
1381	! but its preferable to keep all the calls in one place ...
1382	!-----------------------------------------------------------------------------------------------------
1383
1384
1385	call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
1386	call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
1387	call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
1388	call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
1389	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1390	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
1391	call writediagfi(ngrid,"temp","temperature","K",3,zt)
1392	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
1393	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
1394	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
1395	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
1396	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
1397
1398	! Subsurface temperatures
1399	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1400	! call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
1401
1402
1403
1404	! Total energy balance diagnostics
1405	if(callrad.and.(.not.newtonian))then
1406
1407	!call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
1408	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
1409	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
1410	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
1411
1412	! call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
1413	! call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
1414	! call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
1415	! call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
1416	! call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
1417	! call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
1418
1419
1420	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
1421
1422	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
1423
1424	endif ! end of 'callrad'
1425
1426	if(enertest) then
1427
1428	if (calldifv) then
1429
1430	call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
1431	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
1432
1433	! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
1434	! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
1435	! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
1436
1437	endif
1438
1439	if (corrk) then
1440	call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
1441	call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
1442	endif
1443
1444	endif ! end of 'enertest'
1445
1446	! Temporary inclusions for winds diagnostics.
1447	call writediagfi(ngrid,"zdudif","Turbdiff tend. zon. wind","m s-2",3,zdudif)
1448	call writediagfi(ngrid,"zdudyn","Dyn. tend. zon. wind","m s-2",3,zdudyn)
1449
1450	! Temporary inclusions for heating diagnostics.
1451	call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
1452	call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
1453	call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
1454	call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
1455
1456	! For Debugging.
1457	call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
1458
1459
1460	! Output tracers.
1461	if (tracer) then
1462
1463	if (callmufi) then ! For now we assume an given order for tracers !
1464	call writediagfi(ngrid,"mu_m0as","Spherical mode 0th order moment",'kg/kg',3,zq(:,:,1))
1465	call writediagfi(ngrid,"mu_m3as","Spherical mode 3rd order moment",'kg/kg',3,zq(:,:,2))
1466	call writediagfi(ngrid,"mu_m0af","Fractal mode 0th order moment",'kg/kg',3,zq(:,:,3))
1467	call writediagfi(ngrid,"mu_m3af","Fractal mode 3rd order moment",'kg/kg',3,zq(:,:,4))
1468	endif ! end of 'callmufi'
1469
1470	if ( callchim .and. (.not. callmufi) ) then
1471	call writediagfi(ngrid,"C2H2","C2H2",'kg/kg',3,zq(:,:,10))
1472	call writediagfi(ngrid,"C2H4","C2H4",'kg/kg',3,zq(:,:,12))
1473	call writediagfi(ngrid,"C2H6","C2H6",'kg/kg',3,zq(:,:,14))
1474	call writediagfi(ngrid,"C4H2","C4H2",'kg/kg',3,zq(:,:,26))
1475	call writediagfi(ngrid,"HCN","HCN",'kg/kg',3,zq(:,:,36))
1476	call writediagfi(ngrid,"HC3N","HC3N",'kg/kg',3,zq(:,:,42))
1477	else if ( callchim .and. callmufi ) then
1478	call writediagfi(ngrid,"C2H2","C2H2",'kg/kg',3,zq(:,:,6))
1479	call writediagfi(ngrid,"C2H4","C2H4",'kg/kg',3,zq(:,:,8))
1480	call writediagfi(ngrid,"C2H6","C2H6",'kg/kg',3,zq(:,:,10))
1481	call writediagfi(ngrid,"C4H2","C4H2",'kg/kg',3,zq(:,:,22))
1482	call writediagfi(ngrid,"HCN","HCN",'kg/kg',3,zq(:,:,32))
1483	call writediagfi(ngrid,"HC3N","HC3N",'kg/kg',3,zq(:,:,38))
1484	endif
1485
1486	endif ! end of 'tracer'
1487
1488	! XIOS outputs
1489	#ifdef CPP_XIOS
1490	! Send fields to XIOS: (NB these fields must also be defined as
1491	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
1492	CALL send_xios_field("ls",zls)
1493
1494	CALL send_xios_field("ps",ps)
1495	CALL send_xios_field("area",cell_area)
1496
1497	CALL send_xios_field("temperature",zt)
1498	CALL send_xios_field("u",zu)
1499	CALL send_xios_field("v",zv)
1500
1501	#endif
1502
1503	icount=icount+1
1504
1505	end subroutine physiq
1506
1507	end module physiq_mod

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

Download in other formats:

Original Format