Context Navigation

physiq_mod.F90 @ 1915

Last change on this file since 1915 was 1915, checked in by jvatant, 7 years ago

Add main current send_xios_field outputs
--JVO

Property svn:executable set to *
Property svn:mergeinfo set to
/trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 merged eligible

File size: 63.2 KB

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

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

Download in other formats:

Original Format