Context Navigation

physiq_mod.F90 @ 3497

Last change on this file since 3497 was 3497, checked in by debatzbr, 2 weeks ago

Add AC6H6 condensation in the microphysics

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

File size: 85.5 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, gzlat, grav, BWNV, WAVEI, WAVEV
18	use comchem_h, only: nkim, cnames, nlaykim_up, ykim_up, ykim_tot, botCH4
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 datafile_mod, only: datadir, corrkdir, banddir, haze_opt_file, nudging_file
23	use geometry_mod, only: latitude, latitude_deg, longitude, cell_area
24	USE comgeomfi_h, only: totarea, totarea_planet
25	USE tracer_h
26	use time_phylmdz_mod, only: ecritphy, iphysiq, nday
27	use phyetat0_mod, only: phyetat0
28	use phyredem, only: physdem0, physdem1
29	use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
30	use mod_phys_lmdz_para, only : is_master
31	use planete_mod, only: apoastr, periastr, year_day, peri_day, &
32	obliquit, nres, z0
33	use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
34	use time_phylmdz_mod, only: daysec
35	#ifndef MESOSCALE
36	use logic_mod, only: moyzon_ch
37	use moyzon_mod, only: zphibar, zphisbar, zplevbar, zplaybar, &
38	zzlevbar, zzlaybar, ztfibar, zqfibar
39	#endif
40	use callkeys_mod
41	use phys_state_var_mod
42	use turb_mod, only : q2,sensibFlux,turb_resolved
43	#ifndef MESOSCALE
44	use vertical_layers_mod, only: presnivs, pseudoalt
45	use mod_phys_lmdz_omp_data, ONLY: is_omp_master
46	#else
47	use comm_wrf, only : comm_HR_SW, comm_HR_LW, &
48	comm_FLUXTOP_DN,comm_FLUXABS_SW,&
49	comm_FLUXTOP_LW,comm_FLUXSURF_SW,&
50	comm_FLUXSURF_LW,comm_FLXGRD
51	#endif
52	#ifdef CPP_XIOS
53	use xios_output_mod, only: initialize_xios_output, &
54	update_xios_timestep, &
55	send_xios_field
56	use wxios, only: wxios_context_init, xios_context_finalize
57	#endif
58	use muphy_diag
59	implicit none
60
61
62	!==================================================================
63	!
64	! Purpose
65	! -------
66	! Central subroutine for all the physics parameterisations in the
67	! universal model. Originally adapted from the Mars LMDZ model.
68	!
69	! The model can be run without or with tracer transport
70	! depending on the value of "tracer" in file "callphys.def".
71	!
72	!
73	! It includes:
74	!
75	! I. Initialization :
76	! I.1 Firstcall initializations.
77	! I.2 Initialization for every call to physiq.
78	!
79	! II. Compute radiative transfer tendencies (longwave and shortwave) :
80	! II.a Option 1 : Call correlated-k radiative transfer scheme.
81	! II.b Option 2 : Call Newtonian cooling scheme.
82	! II.c Option 3 : Atmosphere has no radiative effect.
83	!
84	! III. Vertical diffusion (turbulent mixing) :
85	!
86	! IV. Dry Convective adjustment :
87	!
88	! V. Tracers
89	! V.1. Microphysics
90	! V.2. Chemistry
91	! V.3. Updates (pressure variations, surface budget).
92	! V.4. Surface Tracer Update.
93	!
94	! VI. Surface and sub-surface soil temperature.
95	!
96	! VII. Perform diagnostics and write output files.
97	!
98	!
99	! arguments
100	! ---------
101	!
102	! INPUT
103	! -----
104	!
105	! ngrid Size of the horizontal grid.
106	! nlayer Number of vertical layers.
107	! nq Number of advected fields.
108	! nametrac Name of corresponding advected fields.
109	!
110	! firstcall True at the first call.
111	! lastcall True at the last call.
112	!
113	! pday Number of days counted from the North. Spring equinoxe.
114	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
115	! ptimestep timestep (s).
116	!
117	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa).
118	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
119	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2.s-2).
120	!
121	! pu(ngrid,nlayer) u, zonal component of the wind (ms-1).
122	! pv(ngrid,nlayer) v, meridional component of the wind (ms-1).
123	!
124	! pt(ngrid,nlayer) Temperature (K).
125	!
126	! pq(ngrid,nlayer,nq) Advected fields.
127	!
128	! pudyn(ngrid,nlayer) \
129	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
130	! ptdyn(ngrid,nlayer) / corresponding variables.
131	! pqdyn(ngrid,nlayer,nq) /
132	! flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary
133	!
134	! OUTPUT
135	! ------
136	!
137	! pdu(ngrid,nlayer) \
138	! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding
139	! pdt(ngrid,nlayer) / variables due to physical processes.
140	! pdq(ngrid,nlayer) /
141	! pdpsrf(ngrid) /
142	!
143	!
144	! Authors
145	! -------
146	! Frederic Hourdin 15/10/93
147	! Francois Forget 1994
148	! Christophe Hourdin 02/1997
149	! Subroutine completely rewritten by F. Forget (01/2000)
150	! Water ice clouds: Franck Montmessin (update 06/2003)
151	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
152	! New correlated-k radiative scheme: R. Wordsworth (2009)
153	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
154	! Improved water cycle: R. Wordsworth / B. Charnay (2010)
155	! To F90: R. Wordsworth (2010)
156	! New turbulent diffusion scheme: J. Leconte (2012)
157	! Loops converted to F90 matrix format: J. Leconte (2012)
158	! No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
159	! Purge of the code : M. Turbet (2015)
160	! Fork for Titan : J. Vatant d'Ollone (2017)
161	! + clean of all too-generic (ocean, water, co2 ...) routines
162	! + Titan's chemistry
163	! Microphysical moment model: J.Burgalat / B. de Batz de Trenquelléon (2022-2023)
164	! Optics for haze and clouds: B. de Batz de Trenquelléon (2023)
165	!============================================================================================
166
167	! ---------------
168	! 0. DECLARATIONS
169	! ---------------
170
171	include "netcdf.inc"
172
173	! Arguments :
174	! -----------
175
176	! INPUTS:
177	! -------
178
179	integer,intent(in) :: ngrid ! Number of atmospheric columns.
180	integer,intent(in) :: nlayer ! Number of atmospheric layers.
181	integer,intent(in) :: nq ! Number of tracers.
182	character*30,intent(in) :: nametrac(nq) ! Names of the tracers taken from dynamics.
183
184	logical,intent(in) :: firstcall ! Signals first call to physics.
185	logical,intent(in) :: lastcall ! Signals last call to physics.
186
187	real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0.
188	real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
189	real,intent(in) :: ptimestep ! Physics timestep (s).
190	real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
191	real,intent(in) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
192	real,intent(in) :: pphi(ngrid,nlayer) ! Geopotential at mid-layer (m2s-2).
193	real,intent(in) :: pu(ngrid,nlayer) ! Zonal wind component (m/s).
194	real,intent(in) :: pv(ngrid,nlayer) ! Meridional wind component (m/s).
195	real,intent(in) :: pt(ngrid,nlayer) ! Temperature (K).
196	real,intent(in) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
197	real,intent(in) :: flxw(ngrid,nlayer) ! Vertical mass flux (ks/s) at lower boundary of layer
198
199	! OUTPUTS:
200	! --------
201
202	! Physical tendencies :
203
204	real,intent(out) :: pdu(ngrid,nlayer) ! Zonal wind tendencies (m/s/s).
205	real,intent(out) :: pdv(ngrid,nlayer) ! Meridional wind tendencies (m/s/s).
206	real,intent(out) :: pdt(ngrid,nlayer) ! Temperature tendencies (K/s).
207	real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
208	real,intent(out) :: pdpsrf(ngrid) ! Surface pressure tendency (Pa/s).
209
210	! Local saved variables:
211	! ----------------------
212	integer,save :: day_ini ! Initial date of the run (sol since Ls=0).
213	integer,save :: icount ! Counter of calls to physiq during the run.
214	!$OMP THREADPRIVATE(day_ini,icount)
215
216	! Local variables :
217	! -----------------
218
219	real zh(ngrid,nlayer) ! Potential temperature (K).
220	real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
221
222	integer l,ig,ierr,iq,nw,isoil,ilat,lat_idx,i,j
223
224	! FOR DIAGNOSTIC :
225
226	real zls ! Solar longitude (radians).
227	real zlss ! Sub solar point longitude (radians).
228	real zday ! Date (time since Ls=0, calculated in sols).
229	real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers (ref : local surf).
230	real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries (ref : local surf).
231	real zzlay_eff(ngrid,nlayer) ! Effective altitude at the middle of the atmospheric layers (ref : geoid ).
232	real zzlev_eff(ngrid,nlayer+1) ! Effective altitude at the atmospheric layer boundaries ( ref : geoid ).
233
234	! TENDENCIES due to various processes :
235
236	! For Surface Temperature : (K/s)
237	real zdtsurf(ngrid) ! Cumulated tendencies.
238	real zdtsurfmr(ngrid) ! Mass_redistribution routine.
239	real zdtsurfevap(ngrid) ! Evaporation.
240	real zdtsdif(ngrid) ! Turbdiff/vdifc routines.
241
242	! For Atmospheric Temperatures : (K/s)
243	real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
244	real zdtmr(ngrid,nlayer) ! Mass_redistribution routine.
245	real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
246	real zdtlc(ngrid,nlayer) ! Condensation heating rate.
247
248	! For Surface Tracers : (kg/m2/s)
249	real dqsurf(ngrid,nq) ! Cumulated tendencies.
250	real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines.
251	real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine.
252
253	! For Tracers : (kg/kg_of_air/s)
254	real zdqadj(ngrid,nlayer,nq) ! Convadj routine.
255	real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines.
256	real zdqevap(ngrid,nlayer) ! Turbdiff routine.
257	real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine.
258
259	real zdqchi(ngrid,nlayer,nq) ! Chemical tendency (chemistry routine).
260
261	real zdqmufi(ngrid,nlayer,nq) ! Microphysical tendency.
262
263	real zdqfibar(ngrid,nlayer,nq) ! For 2D chemistry
264	real zdqmufibar(ngrid,nlayer,nq) ! For 2D chemistry
265
266	! For Winds : (m/s/s)
267	real zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! Convadj routine.
268	real zdumr(ngrid,nlayer),zdvmr(ngrid,nlayer) ! Mass_redistribution routine.
269	real zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
270	real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines.
271	real zdhadj(ngrid,nlayer) ! Convadj routine.
272	real zdundg(ngrid,nlayer) ! Nudging for zonal wind.
273
274	! For Pressure and Mass :
275	real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
276	real zdmassmr_col(ngrid) ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
277	real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s).
278
279
280	! Local variables for LOCAL CALCULATIONS:
281	! ---------------------------------------
282	real zflubid(ngrid)
283	real zplanck(ngrid),zpopsk(ngrid,nlayer)
284	real ztim1,ztim2,ztim3, z1,z2
285	real ztime_fin
286	real zdh(ngrid,nlayer)
287	real gmplanet
288	real taux(ngrid),tauy(ngrid)
289	real factlat
290	real zundg(nlayer)
291
292
293	! local variables for DIAGNOSTICS : (diagfi & stat)
294	! -------------------------------------------------
295	real ps(ngrid) ! Surface Pressure.
296	real zt(ngrid,nlayer) ! Atmospheric Temperature.
297	real zu(ngrid,nlayer),zv(ngrid,nlayer) ! Zonal and Meridional Winds.
298	real zq(ngrid,nlayer,nq) ! Atmospheric Tracers.
299	real zdtadj(ngrid,nlayer) ! Convadj Diagnostic.
300	real zdtdyn(ngrid,nlayer) ! Dynamical Heating (K/s).
301	real zdudyn(ngrid,nlayer) ! Dynamical Zonal Wind tendency (m.s-2).
302
303	real zhorizwind(ngrid,nlayer) ! Horizontal Wind (sqrt(uu+vv))
304
305	real vmr(ngrid,nlayer) ! volume mixing ratio
306	real time_phys
307
308	real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
309
310	! To test energy conservation (RW+JL)
311	real mass(ngrid,nlayer),massarea(ngrid,nlayer)
312	real dEtot, dEtots, AtmToSurf_TurbFlux
313	real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
314	!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
315	real dEzRadsw(ngrid,nlayer),dEzRadlw(ngrid,nlayer),dEzdiff(ngrid,nlayer)
316	real dEdiffs(ngrid),dEdiff(ngrid)
317
318	! JL12 : conservation test for mean flow kinetic energy has been disabled temporarily
319	real dItot, dItot_tmp, dVtot, dVtot_tmp
320	real dWtot, dWtot_tmp, dWtots, dWtots_tmp
321
322	! For Clear Sky Case.
323	real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxsurfabs_sw1(ngrid) ! For SW/LW flux.
324	real fluxtop_lw1(ngrid), fluxabs_sw1(ngrid) ! For SW/LW flux.
325	real albedo_equivalent1(ngrid) ! For Equivalent albedo calculation.
326	real tf, ntf
327
328	! Miscellaneous :
329	character(len=10) :: tmp1
330	character(len=10) :: tmp2
331	character*2 :: str2
332
333	#ifndef MESOSCALE
334
335	! Local variables for Titan chemistry and microphysics
336	! ----------------------------------------------------
337
338	real,save :: ctimestep ! Chemistry timestep (s)
339	!$OMP THREADPRIVATE(ctimestep)
340
341	! Chemical tracers in molar fraction [mol/mol]
342	real, dimension(ngrid,nlayer,nkim) :: ychim
343	real, dimension(ngrid,nlayer,nkim) :: ychimbar ! For 2D chemistry
344
345	! Molar fraction tendencies (chemistry, condensation and evaporation) for tracers [mol/mol/s]
346	real, dimension(ngrid,nlayer,nq) :: dyccond ! Condensation rate. NB : for all tracers, as we want to use indx on it.
347	real, dimension(ngrid,nlayer,size(ices_indx)) :: dmuficond ! Condensation rate from microphysics [kg/kg/s].
348	real, dimension(ngrid,nlayer,nq) :: dyccondbar ! For 2D chemistry
349	real, dimension(ngrid) :: dycevapCH4 ! Surface "pseudo-evaporation" rate for CH4.
350
351	! Saturation profiles [mol/mol]
352	real, dimension(ngrid,nlayer,nkim) :: ysat
353
354	! Temporary fraction of CH4 [mol/mol]
355	real, dimension(ngrid,nlayer) :: tpq_CH4
356
357	real :: i2e(ngrid,nlayer) ! int 2 ext factor (X.kg-1 -> X.m-3)
358
359	#ifdef USE_QTEST
360	real,save,dimension(:,:,:), allocatable :: tpq ! Tracers for decoupled microphysical tests ( temporary in 01/18 )
361	!$OMP THREADPRIVATE(tpq)
362	real,dimension(ngrid,nlayer,nq) :: dtpq ! (temporary in 01/18)
363	#endif
364
365	logical file_ok
366
367	!-----------------------------------------------------------------------------
368	! Interface to calmufi
369	! --> needed in order to pass assumed-shape arrays. Otherwise we must put calmufi in a module
370	! (to have an explicit interface generated by the compiler).
371	! Or one can put calmufi in MMP_GCM module (in muphytitan).
372	INTERFACE
373	SUBROUTINE calmufi(dt, plev, zlev, play, zlay, g3d, temp, pq, zdqfi, zdq)
374	REAL(kind=8), INTENT(IN) :: dt !! Physics timestep (s).
375	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: plev !! Pressure levels (Pa).
376	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlev !! Altitude levels (m).
377	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: play !! Pressure layers (Pa).
378	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlay !! Altitude at the center of each layer (m).
379	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: g3d !! Latitude-Altitude depending gravitational acceleration (m.s-2).
380	REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp !! Temperature at the center of each layer (K).
381	REAL(kind=8), DIMENSION(:,:,:), INTENT(IN) :: pq !! Tracers ($X.kg^{-1}}$).
382	REAL(kind=8), DIMENSION(:,:,:), INTENT(IN) :: zdqfi !! Tendency from former processes for tracers ($X.kg^{-1}}$).
383	REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq !! Microphysical tendency for tracers ($X.kg^{-1}}$).
384	END SUBROUTINE calmufi
385	END INTERFACE
386
387	#endif
388
389	!==================================================================================================
390
391	! -----------------
392	! I. INITIALISATION
393	! -----------------
394
395	! --------------------------------
396	! I.1 First Call Initialisation.
397	! --------------------------------
398	if (firstcall) then
399
400	#ifdef USE_QTEST
401	allocate(tpq(ngrid,nlayer,nq))
402	tpq(:,:,:) = pq(:,:,:)
403	#endif
404	! Initialisation of nmicro as well as tracers names, indexes ...
405	if (ngrid.ne.1) then ! Already done in rcm1d
406	call initracer2(nq,nametrac) ! WARNING JB (27/03/2018): should be wrapped in an OMP SINGLE statement (see module notes)
407	endif
408
409	! Allocate saved arrays (except for 1D model, where this has already been done)
410	#ifndef MESOSCALE
411	if (ngrid>1) call phys_state_var_init(nq)
412	#endif
413
414	! Variables set to 0
415	! ~~~~~~~~~~~~~~~~~~
416	dtrad(:,:) = 0.D0
417	zdtlc(:,:) = 0.D0
418	fluxrad(:) = 0.D0
419	zdtsw(:,:) = 0.D0
420	zdtlw(:,:) = 0.D0
421	zpopthi(:,:,:,:) = 0.D0
422	zpopthv(:,:,:,:) = 0.D0
423	zpoptti(:,:,:,:) = 0.D0
424	zpopttv(:,:,:,:) = 0.D0
425
426	! Initialize setup for correlated-k radiative transfer
427	! JVO 17 : Was in callcorrk firstcall, but we need spectral intervals for microphysics.
428	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
429
430	if (corrk) then
431
432	call system('rm -f surf_vals_long.out')
433
434	write( tmp1, '(i3)' ) L_NSPECTI
435	write( tmp2, '(i3)' ) L_NSPECTV
436	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
437	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
438
439	call setspi ! Basic infrared properties.
440	call setspv ! Basic visible properties.
441	call sugas_corrk ! Set up gaseous absorption properties.
442
443	OLR_nu(:,:) = 0.D0
444	OSR_nu(:,:) = 0.D0
445
446	int_dtaui(:,:,:) = 0.D0
447	int_dtauv(:,:,:) = 0.D0
448
449	#ifndef MESOSCALE
450	IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
451	haze_opt_file=trim(datadir)//'/optical_tables/HAZE_OPTIC_'//trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))//'.DAT'
452	inquire(file=trim(haze_opt_file),exist=file_ok)
453	if(.not.file_ok) then
454	write(,) 'The file ',TRIM(haze_opt_file),' with the haze optical properties'
455	write(,) 'was not found by optci.F90 ! Check in ', TRIM(datadir)
456	write(,) 'that you have the one corresponding to the given spectral resolution !!'
457	write(,) 'Meanwhile I abort ...'
458	call abort
459	endif
460	ENDIF
461	#endif
462
463	endif
464
465	#ifndef MESOSCALE
466	! Initialize names and timestep for chemistry
467	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
468
469	if (callchim) then
470
471	if (moyzon_ch .and. ngrid.eq.1) then
472	print *, "moyzon_ch=",moyzon_ch," and ngrid=1"
473	print *, "Please desactivate zonal mean for 1D !"
474	stop
475	endif
476
477	! Chemistry timestep
478	ctimestep = ptimestep*REAL(ichim)
479
480	endif
481
482	! Initialize microphysics.
483	! ~~~~~~~~~~~~~~~~~~~~~~~~
484
485	IF (callmufi) THEN
486	! WARNING JB (27/03/2018): inimufi should be wrapped in an OMP SINGLE statement.
487	call inimufi(ptimestep)
488
489	! initialize microphysics diagnostics arrays.
490	call ini_diag_arrays(ngrid,nlayer,nice)
491
492	ENDIF
493	#endif
494
495	#ifdef CPP_XIOS
496	! Initialize XIOS context
497	write(,) "physiq: call wxios_context_init"
498	CALL wxios_context_init
499	#endif
500
501	! Read 'startfi.nc' file.
502	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
503	#ifndef MESOSCALE
504	call phyetat0(startphy_file, &
505	ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, &
506	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,tankCH4)
507	#else
508	emis(:)=0.0
509	q2(:,:)=0.0
510	qsurf(:,:)=0.0
511	tankCH4(:)=0.0
512	day_ini = pday
513	#endif
514
515	#ifndef MESOSCALE
516	if (.not.startphy_file) then
517	! additionnal "academic" initialization of physics
518	if (is_master) write(,) "Physiq: initializing tsurf(:) to pt(:,1) !!"
519	tsurf(:)=pt(:,1)
520	if (is_master) write(,) "Physiq: initializing tsoil(:) to pt(:,1) !!"
521	do isoil=1,nsoilmx
522	tsoil(:,isoil)=tsurf(:)
523	enddo
524	if (is_master) write(,) "Physiq: initializing day_ini to pdat !"
525	day_ini=pday
526	endif
527	#endif
528
529	if (pday.ne.day_ini) then
530	write(,) "ERROR in physiq.F90:"
531	write(,) "bad synchronization between physics and dynamics"
532	write(,) "dynamics day: ",pday
533	write(,) "physics day: ",day_ini
534	stop
535	endif
536	write (,) 'In physiq day_ini =', day_ini
537
538	! Initialize albedo calculation.
539	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
540	albedo(:,:)=0.D0
541	albedo_bareground(:)=0.D0
542	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground)
543
544	! Initialize orbital calculation.
545	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
546	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
547
548
549	if(tlocked)then
550	print*,'Planet is tidally locked at resonance n=',nres
551	print*,'Make sure you have the right rotation rate!!!'
552	endif
553
554
555	! Initialize soil.
556	! ~~~~~~~~~~~~~~~~
557	if (callsoil) then
558
559	call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
560	ptimestep,tsurf,tsoil,capcal,fluxgrd)
561
562	else ! else of 'callsoil'.
563
564	print*,'WARNING! Thermal conduction in the soil turned off'
565	capcal(:)=1.e6
566	fluxgrd(:)=intheat
567	print*,'Flux from ground = ',intheat,' W m^-2'
568
569	endif ! end of 'callsoil'.
570
571	icount=1
572
573
574	! Initialize surface history variable.
575	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
576	qsurf_hist(:,:)=qsurf(:,:)
577
578	! Initialize variable for dynamical heating and zonal wind tendency diagnostic
579	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
580	ztprevious(:,:)=pt(:,:)
581	zuprevious(:,:)=pu(:,:)
582
583
584	if(meanOLR)then
585	call system('rm -f rad_bal.out') ! to record global radiative balance.
586	call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
587	call system('rm -f h2o_bal.out') ! to record global hydrological balance.
588	endif
589
590	! Read NewNudging.dat and initialize the nudging for pu
591	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
592	if (nudging_u) then
593	nudging_file=trim(datadir)//'/NewNudging.dat'
594	inquire(file=trim(nudging_file),exist=file_ok)
595	if(.not.file_ok) then
596	write(,) 'ERROR : The file ',TRIM(nudging_file),' was not found by physiq_mod.F90 ! Check in ', TRIM(datadir)
597	write(,) 'Meanwhile I abort ...'
598	call abort
599	endif
600
601	open(88,file=nudging_file,form='formatted')
602	do ilat = 1, 49
603	read(88,*) u_ref(ilat,:)
604	enddo
605	close(88)
606	endif
607
608	#ifndef MESOSCALE
609	if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
610	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
611	ptimestep,pday+nday,time_phys,cell_area, &
612	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
613	endif
614	#endif
615
616	! XIOS outputs
617	#ifdef CPP_XIOS
618
619	write(,) "physiq: call initialize_xios_output"
620	call initialize_xios_output(pday,ptime,ptimestep,daysec, &
621	presnivs,pseudoalt,wavei,wavev)
622	#endif
623	write(,) "physiq: end of firstcall"
624	endif ! end of 'firstcall'
625
626	! ------------------------------------------------------
627	! I.2 Initializations done at every physical timestep:
628	! ------------------------------------------------------
629
630	#ifdef CPP_XIOS
631	! update XIOS time/calendar
632	call update_xios_timestep
633	#endif
634
635	! Initialize various variables
636	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
637
638	pdt(:,:) = 0.D0
639	zdtsurf(:) = 0.D0
640	zdtsurfevap(:) = 0.D0
641	pdq(:,:,:) = 0.D0
642	dqsurf(:,:) = 0.D0
643	pdu(:,:) = 0.D0
644	pdv(:,:) = 0.D0
645	pdpsrf(:) = 0.D0
646	zflubid(:) = 0.D0
647	taux(:) = 0.D0
648	tauy(:) = 0.D0
649
650	zday=pday+ptime ! Compute time, in sols (and fraction thereof).
651
652	! Compute Stellar Longitude (Ls), and orbital parameters.
653	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
654	if (season) then
655	call stellarlong(zday,zls)
656	else
657	call stellarlong(noseason_day,zls)
658	end if
659
660	call orbite(zls,dist_star,declin,right_ascen)
661
662	if (tlocked) then
663	zlss=Mod(-(2.pi(zday/year_day)nres - right_ascen),2.pi)
664	elseif (diurnal) then
665	zlss=-2.pi(zday-.5)
666	else if(diurnal .eqv. .false.) then
667	zlss=9999.
668	endif
669
670
671	! Compute variations of g with latitude (oblate case).
672	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
673	if (oblate .eqv. .false.) then
674	gzlat(:,:) = g
675	else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
676	print*,'I need values for flatten, J2, Rmean and MassPlanet to compute gzlat (else set oblate=.false.)'
677	call abort
678	else
679	gmplanet = MassPlanetgrav1e24
680	do ig=1,ngrid
681	gzlat(ig,:)= gmplanet/(Rmean*2) (1.D0 + 0.75 J2 - 2.0flatten/3. + (2.flatten - 15./4. J2) * cos(2. * (pi/2. - latitude(ig))))
682	end do
683	endif
684
685	! Compute altitudes with the geopotential coming from the dynamics.
686	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
687
688	if (eff_gz .eqv. .false.) then
689
690	do l=1,nlayer
691	zzlay(:,l) = pphi(:,l) / gzlat(:,l) ! Reference = local surface
692	enddo
693
694	else ! In this case we consider variations of g with altitude
695
696	do l=1,nlayer
697	zzlay(:,l) = gradrad / ( g*rad - ( pphi(:,l) + phisfi(:) )) - rad
698	gzlat(:,l) = gradrad / ( rad + zzlay(:,l) )**2
699	end do
700
701	endif ! if eff_gz
702
703	zzlev(:,1)=0.
704	zzlev(:,nlayer+1)=1.e7 ! Dummy top of last layer above 10000 km...
705	! JVO 19 : This altitude is indeed dummy for the GCM and fits ptop=0
706	! but for upper chemistry that's a pb -> we anyway redefine it just after ..
707
708	do l=2,nlayer
709	do ig=1,ngrid
710	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
711	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
712	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
713	enddo
714	enddo
715
716	! Effective altitudes ( eg needed for chemistry ) with correct g, and with reference to the geoid
717	! JVO 19 : We shall always have correct altitudes in chemistry no matter what's in physics
718	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
719	#ifndef MESOSCALE
720	if (moyzon_ch) then ! Zonal averages
721
722	zzlaybar(1,:)=gradrad/(g*rad-(zphibar(1,:)+zphisbar(1)))-rad ! reference = geoid
723	zzlevbar(1,1)=zphisbar(1)/g
724	DO l=2,nlayer
725	z1=(zplaybar(1,l-1)+zplevbar(1,l))/(zplaybar(1,l-1)-zplevbar(1,l))
726	z2=(zplevbar(1,l) +zplaybar(1,l))/(zplevbar(1,l) -zplaybar(1,l))
727	zzlevbar(1,l)=(z1zzlaybar(1,l-1)+z2zzlaybar(1,l))/(z1+z2)
728	ENDDO
729	zzlevbar(1,nlayer+1)=zzlaybar(1,nlayer)+(zzlaybar(1,nlayer)-zzlevbar(1,nlayer))
730
731	DO ig=2,ngrid
732	if (latitude(ig).ne.latitude(ig-1)) then
733	DO l=1,nlayer
734	zzlaybar(ig,l)=gradrad/(g*rad-(zphibar(ig,l)+zphisbar(ig)))-rad
735	ENDDO
736	zzlevbar(ig,1)=zphisbar(ig)/g
737	DO l=2,nlayer
738	z1=(zplaybar(ig,l-1)+zplevbar(ig,l))/ (zplaybar(ig,l-1)-zplevbar(ig,l))
739	z2=(zplevbar(ig,l) +zplaybar(ig,l))/(zplevbar(ig,l) -zplaybar(ig,l))
740	zzlevbar(ig,l)=(z1zzlaybar(ig,l-1)+z2zzlaybar(ig,l))/(z1+z2)
741	ENDDO
742	zzlevbar(ig,nlayer+1)=zzlaybar(ig,nlayer)+(zzlaybar(ig,nlayer)-zzlevbar(ig,nlayer))
743	else
744	zzlaybar(ig,:)=zzlaybar(ig-1,:)
745	zzlevbar(ig,:)=zzlevbar(ig-1,:)
746	endif
747	ENDDO
748
749	else ! if not moyzon
750	#endif
751
752	DO ig=1,ngrid
753	DO l=1,nlayer
754	zzlay_eff(ig,l)=gradrad/(g*rad-(pphi(ig,l)+phisfi(ig)))-rad ! reference = geoid
755	ENDDO
756	zzlev_eff(ig,1)=phisfi(ig)/g
757	DO l=2,nlayer
758	z1=(pplay(ig,l-1)+pplev(ig,l))/ (pplay(ig,l-1)-pplev(ig,l))
759	z2=(pplev(ig,l) +pplay(ig,l))/(pplev(ig,l) -pplay(ig,l))
760	zzlev_eff(ig,l)=(z1zzlay_eff(ig,l-1)+z2zzlay_eff(ig,l))/(z1+z2)
761	ENDDO
762	zzlev_eff(ig,nlayer+1)=zzlay_eff(ig,nlayer)+(zzlay_eff(ig,nlayer)-zzlev_eff(ig,nlayer))
763	ENDDO
764
765	#ifndef MESOSCALE
766	endif ! moyzon
767	#endif
768
769	! -------------------------------------------------------------------------------------
770	! Compute potential temperature
771	! Note : Potential temperature calculation may not be the same in physiq and dynamic...
772	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
773	do l=1,nlayer
774	do ig=1,ngrid
775	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
776	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
777	mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/gzlat(ig,l)
778	massarea(ig,l)=mass(ig,l)*cell_area(ig)
779	enddo
780	enddo
781
782	! Compute vertical velocity (m/s) from vertical mass flux
783	! w = F / (rhoarea) and rho = P/(rT)
784	! But first linearly interpolate mass flux to mid-layers
785	do l=1,nlayer-1
786	pw(:,l)=0.5*(flxw(:,l)+flxw(:,l+1))
787	enddo
788	pw(:,nlayer)=0.5*flxw(:,nlayer) ! since flxw(nlayer+1)=0
789	do l=1,nlayer
790	pw(:,l)=(pw(:,l)rpt(:,l)) / (pplay(:,l)*cell_area(:))
791	enddo
792
793	!---------------------------------
794	! II. Compute radiative tendencies
795	!---------------------------------
796
797	if (callrad) then
798	if( mod(icount-1,iradia).eq.0.or.lastcall) then
799
800	! Compute local stellar zenith angles
801	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
802	if (tlocked) then
803	! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb
804	ztim1=SIN(declin)
805	ztim2=COS(declin)*COS(zlss)
806	ztim3=COS(declin)*SIN(zlss)
807
808	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
809	ztim1,ztim2,ztim3,mu0,fract, flatten)
810
811	elseif (diurnal) then
812	ztim1=SIN(declin)
813	ztim2=COS(declin)COS(2.pi*(zday-.5))
814	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
815
816	call stelang(ngrid,sinlon,coslon,sinlat,coslat, &
817	ztim1,ztim2,ztim3,mu0,fract, flatten)
818
819	else if(diurnal .eqv. .false.) then
820
821	call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad,flatten)
822	! WARNING: this function appears not to work in 1D
823
824	endif
825
826	! Eclipse incoming sunlight (e.g. Saturn ring shadowing).
827	if(rings_shadow) then
828	call call_rings(ngrid, ptime, pday, diurnal)
829	endif
830
831
832	if (corrk) then
833
834	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
835	! II.a Call correlated-k radiative transfer scheme
836	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
837
838	call call_profilgases(nlayer)
839
840	! standard callcorrk
841	call callcorrk(ngrid,nlayer,pq,nq,qsurf,zday, &
842	albedo,albedo_equivalent,emis,mu0,pplev,pplay,zzlev,&
843	pt,tsurf,fract,dist_star, &
844	zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, &
845	fluxsurfabs_sw,fluxtop_lw, &
846	fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, &
847	int_dtaui,int_dtauv,zpopthi,zpopthv,zpoptti,zpopttv,&
848	lastcall)
849
850	! Radiative flux from the sky absorbed by the surface (W.m-2).
851	GSR=0.0
852	fluxrad_sky(:)=emis(:)*fluxsurf_lw(:)+fluxsurfabs_sw(:)
853
854	!if(noradsurf)then ! no lower surface; SW flux just disappears
855	! GSR = SUM(fluxsurf_sw(:)*cell_area(:))/totarea
856	! fluxrad_sky(:)=emis(:)*fluxsurf_lw(:)
857	! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
858	!endif
859
860	! Net atmospheric radiative heating rate (K.s-1)
861	dtrad(:,:)=zdtsw(:,:)+zdtlw(:,:)
862
863	elseif(newtonian)then
864
865	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
866	! II.b Call Newtonian cooling scheme
867	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
868	call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
869
870	zdtsurf(:) = +(pt(:,1)-tsurf(:))/ptimestep
871	! e.g. surface becomes proxy for 1st atmospheric layer ?
872
873	else
874
875	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
876	! II.c Atmosphere has no radiative effect
877	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
878	fluxtop_dn(:) = fract(:)mu0(:)Fat1AU/dist_star**2
879	if(ngrid.eq.1)then ! / by 4 globally in 1D case...
880	fluxtop_dn(1) = fract(1)Fat1AU/dist_star*2/2.0
881	endif
882	fluxsurf_sw(:) = fluxtop_dn(:)
883	print*,'------------WARNING---WARNING------------' ! by MT2015.
884	print*,'You are in corrk=false mode, '
885	print*,'and the surface albedo is taken equal to the first visible spectral value'
886
887	fluxsurfabs_sw(:) = fluxtop_dn(:)*(1.-albedo(:,1))
888	fluxrad_sky(:) = fluxsurfabs_sw(:)
889	fluxtop_lw(:) = emis(:)sigmatsurf(:)**4
890
891	dtrad(:,:)=0.D0 ! no atmospheric radiative heating
892
893	endif ! end of corrk
894
895	endif ! of if(mod(icount-1,iradia).eq.0)
896
897
898	! Transformation of the radiative tendencies
899	! ------------------------------------------
900	zplanck(:)=tsurf(:)*tsurf(:)
901	zplanck(:)=emis(:)sigmazplanck(:)*zplanck(:)
902	fluxrad(:)=fluxrad_sky(:)-zplanck(:)
903	pdt(:,:)=pdt(:,:)+dtrad(:,:)
904
905	! Test of energy conservation
906	!----------------------------
907	if(enertest)then
908	call planetwide_sumval(cppmassarea(:,:)zdtsw(:,:)/totarea_planet,dEtotSW)
909	call planetwide_sumval(cppmassarea(:,:)zdtlw(:,:)/totarea_planet,dEtotLW)
910	!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
911	call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
912	call planetwide_sumval((fluxsurf_lw(:)emis(:)-zplanck(:))cell_area(:)/totarea_planet,dEtotsLW)
913	dEzRadsw(:,:)=cppmass(:,:)zdtsw(:,:)
914	dEzRadlw(:,:)=cppmass(:,:)zdtlw(:,:)
915	if (is_master) then
916	print*,'---------------------------------------------------------------'
917	print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2'
918	print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2'
919	print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2'
920	print*,'In corrk SW surface heating =',dEtotsSW,' W m-2'
921	print*,'In corrk LW surface heating =',dEtotsLW,' W m-2'
922	print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2'
923	endif
924	endif ! end of 'enertest'
925
926	endif ! of if (callrad)
927
928
929
930	! --------------------------------------------
931	! III. Vertical diffusion (turbulent mixing) :
932	! --------------------------------------------
933	if (calldifv) then
934
935	zflubid(:)=fluxrad(:)+fluxgrd(:)
936
937	! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
938	if (UseTurbDiff) then
939
940	call turbdiff(ngrid,nlayer,nq, &
941	ptimestep,capcal,lwrite, &
942	pplay,pplev,zzlay,zzlev,z0, &
943	pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, &
944	pdt,pdq,zflubid, &
945	zdudif,zdvdif,zdtdif,zdtsdif, &
946	sensibFlux,q2,zdqdif,zdqsdif, &
947	taux,tauy,lastcall)
948
949	else
950
951	zdh(:,:)=pdt(:,:)/zpopsk(:,:)
952
953	call vdifc(ngrid,nlayer,nq,zpopsk, &
954	ptimestep,capcal,lwrite, &
955	pplay,pplev,zzlay,zzlev,z0, &
956	pu,pv,zh,pq,tsurf,emis,qsurf, &
957	zdh,pdq,zflubid, &
958	zdudif,zdvdif,zdhdif,zdtsdif, &
959	sensibFlux,q2,zdqdif,zdqsdif, &
960	taux,tauy,lastcall)
961
962	zdtdif(:,:)=zdhdif(:,:)*zpopsk(:,:) ! for diagnostic only
963	zdqevap(:,:)=0.
964
965	end if !end of 'UseTurbDiff'
966
967	if (.not. turb_resolved) then
968	pdv(:,:)=pdv(:,:)+zdvdif(:,:)
969	pdu(:,:)=pdu(:,:)+zdudif(:,:)
970	pdt(:,:)=pdt(:,:)+zdtdif(:,:)
971	endif
972
973	zdtsurf(:)=zdtsurf(:)+zdtsdif(:)
974
975	if (tracer) then
976	pdq(:,:,:)=pdq(:,:,:)+ zdqdif(:,:,:)
977	dqsurf(:,:)=dqsurf(:,:) + zdqsdif(:,:)
978	end if ! of if (tracer)
979
980
981	! test energy conservation
982	!-------------------------
983	if(enertest)then
984
985	dEzdiff(:,:)=cppmass(:,:)zdtdif(:,:)
986	do ig = 1, ngrid
987	dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
988	dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
989	enddo
990
991	call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
992	dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
993	call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
994	call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
995
996	if (is_master) then
997
998	if (UseTurbDiff) then
999	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1000	print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2'
1001	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1002	else
1003	print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1004	print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2'
1005	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1006	end if
1007	endif ! end of 'is_master'
1008
1009	! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
1010	endif ! end of 'enertest'
1011
1012	else ! calldifv
1013
1014	if(.not.newtonian)then
1015
1016	zdtsurf(:) = zdtsurf(:) + (fluxrad(:) + fluxgrd(:))/capcal(:)
1017
1018	endif
1019
1020	endif ! end of 'calldifv'
1021
1022
1023	!----------------------------------
1024	! IV. Dry convective adjustment :
1025	!----------------------------------
1026
1027	if(calladj) then
1028
1029	zdh(:,:) = pdt(:,:)/zpopsk(:,:)
1030	zduadj(:,:)=0.D0
1031	zdvadj(:,:)=0.D0
1032	zdhadj(:,:)=0.D0
1033
1034
1035	call convadj(ngrid,nlayer,nq,ptimestep, &
1036	pplay,pplev,zpopsk, &
1037	pu,pv,zh,pq, &
1038	pdu,pdv,zdh,pdq, &
1039	zduadj,zdvadj,zdhadj, &
1040	zdqadj)
1041
1042	pdu(:,:) = pdu(:,:) + zduadj(:,:)
1043	pdv(:,:) = pdv(:,:) + zdvadj(:,:)
1044	pdt(:,:) = pdt(:,:) + zdhadj(:,:)*zpopsk(:,:)
1045	zdtadj(:,:) = zdhadj(:,:)*zpopsk(:,:) ! for diagnostic only
1046
1047	if(tracer) then
1048	pdq(:,:,:) = pdq(:,:,:) + zdqadj(:,:,:)
1049	end if
1050
1051	! Test energy conservation
1052	if(enertest)then
1053	call planetwide_sumval(cppmassarea(:,:)zdtadj(:,:)/totarea_planet,dEtot)
1054	if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2'
1055	endif
1056
1057
1058	endif ! end of 'calladj'
1059
1060
1061	!---------------------------------------------
1062	! V. Specific parameterizations for tracers
1063	!---------------------------------------------
1064
1065	if (tracer) then
1066
1067	#ifndef MESOSCALE
1068	! -------------------
1069	! V.1 Microphysics
1070	! -------------------
1071
1072	! We must call microphysics before chemistry, for condensation !
1073	if (callmufi) then
1074
1075	zzlev(:,nlayer+1)=zzlay(:,nlayer)+(zzlay(:,nlayer)-zzlev(:,nlayer)) ! JVO 19 : We assume zzlev isn't reused later on (could be done cleaner)
1076
1077	#ifdef USE_QTEST
1078	dtpq(:,:,:) = 0.D0 ! we want tpq to go only through mufi
1079	call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,tpq,dtpq,zdqmufi)
1080	tpq(:,:,:) = tpq(:,:,:) + zdqmufi(:,:,:)ptimestep ! only manipulation of tpq->ptimestep here
1081
1082	#else
1083	call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,zdqmufi)
1084	pdq(:,:,:) = pdq(:,:,:) + zdqmufi(:,:,:)
1085
1086	! Sanity check (way safer to be done here rather than within YAMMS)
1087	! Important : the sanity check intentionally include the former processes tendency !
1088	! NB : Despite this sanity check there might be still some unphysical values going through :
1089	! - Negatives, but harmless as it will be only for the output files
1090	! just remove them in post-proc.
1091	! - Weird unphysical ratio of m0 and m3, ok for now, but take care of them if
1092	! you want to compute optics from radii.
1093	WHERE ((pq(:,:,1)+pdq(:,:,1)ptimestep < 0.D0) .OR. (pq(:,:,2)+pdq(:,:,2)ptimestep < 0.D0))
1094	pdq(:,:,1) = (epsilon(1.0)-1.D0)*pq(:,:,1)/ptimestep
1095	pdq(:,:,2) = (epsilon(1.0)-1.D0)*pq(:,:,2)/ptimestep
1096	ENDWHERE
1097	WHERE ((pq(:,:,3)+pdq(:,:,3)ptimestep < 0.D0) .OR. (pq(:,:,4)+pdq(:,:,4)ptimestep < 0.D0))
1098	pdq(:,:,3) = (epsilon(1.0)-1.D0)*pq(:,:,3)/ptimestep
1099	pdq(:,:,4) = (epsilon(1.0)-1.D0)*pq(:,:,4)/ptimestep
1100	ENDWHERE
1101	IF (callclouds) THEN
1102	WHERE ((pq(:,:,5)+pdq(:,:,5)ptimestep < 0.D0) .OR. (pq(:,:,6)+pdq(:,:,6)ptimestep < 0.D0))
1103	pdq(:,:,5) = (epsilon(1.0)-1.D0)*pq(:,:,5)/ptimestep
1104	pdq(:,:,6) = (epsilon(1.0)-1.D0)*pq(:,:,6)/ptimestep
1105	ENDWHERE
1106	DO iq = 1, size(ices_indx)
1107	! For ices :
1108	WHERE (pq(:,:,ices_indx(iq))+pdq(:,:,ices_indx(iq))*ptimestep < 0.D0)
1109	pdq(:,:,ices_indx(iq)) = (epsilon(1.0)-1.D0)*pq(:,:,ices_indx(iq))/ptimestep
1110	ENDWHERE
1111	! For gases :
1112	WHERE (pq(:,:,gazs_indx(iq))+pdq(:,:,gazs_indx(iq))*ptimestep < 0.D0)
1113	pdq(:,:,gazs_indx(iq)) = (epsilon(1.0)-1.D0)*pq(:,:,gazs_indx(iq))/ptimestep
1114	ENDWHERE
1115	ENDDO
1116	ENDIF
1117	! In case there is no clouds, in the troposphere the moments are fixed to evacuate all aerosols
1118	IF (.NOT. callclouds) THEN
1119	WHERE (pplay(:,:) .gt. 1000.)
1120	pdq(:,:,1) = 0.
1121	pdq(:,:,2) = 0.
1122	ENDWHERE
1123	WHERE (pplay(:,:) .gt. 5000.)
1124	pdq(:,:,3) = 0.
1125	pdq(:,:,4) = 0.
1126	ENDWHERE
1127	ENDIF
1128	#endif
1129
1130	! Microphysics condensation for 2D fields to sent non-saturated fields to photochem
1131	if (callclouds .and. moyzon_ch .and. mod(icount-1,ichim).eq.0) then
1132	zdqfibar(:,:,:) = 0.D0 ! We work in zonal average -> forget processes other than condensation
1133	call calmufi(ptimestep,zplevbar,zzlevbar,zplaybar,zzlaybar, &
1134	gzlat,ztfibar,zqfibar,zdqfibar,zdqmufibar)
1135	! TODO : Add a sanity check here !
1136	endif
1137
1138	! Condensation heating rate :
1139	if (callclouds) then
1140	! Default value -> no condensation [kg/kg_air/s] :
1141	dmuficond(:,:,:) = 0.D0
1142	do iq = 1, size(ices_indx)
1143	dmuficond(:,:,iq) = zdqmufi(:,:,gazs_indx(iq))
1144	enddo
1145	call cond_muphy(ngrid,nlayer,pt,dmuficond,zdtlc)
1146	!pdt(:,:) = pdt(:,:) + zdtlc(:,:)
1147	endif
1148	endif ! callmufi
1149
1150	! -----------------
1151	! V.2. Chemistry
1152	! -----------------
1153	! NB : Must be call last ( brings fields back to an equilibrium )
1154	if (callchim) then
1155
1156	! o. Convert updated tracers to molar fraction
1157	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1158	do iq = 1,nkim
1159	ychim(:,:,iq) = (pq(:,:,iq+nmicro) + pdq(:,:,iq+nmicro)*ptimestep) / rat_mmol(iq+nmicro)
1160	enddo
1161
1162	! JVO 05/18 : We update zonally averaged fields with condensation
1163	! as it is compulsory to have correct photochem production. But for other
1164	! processes ( convadj ... ) we miss them in any case as we work in zonally/diurnal
1165	! mean -> no fine diurnal/short time evolution, only seasonal evolution only.
1166	if (moyzon_ch .and. mod(icount-1,ichim).eq. 0) then
1167	do iq = 1,nkim
1168	ychimbar(:,:,iq) = zqfibar(:,:,iq+nmicro) / rat_mmol(iq+nmicro)
1169	if ( callclouds ) then
1170	ychimbar(:,:,iq) = ychimbar(:,:,iq) + ( zdqmufibar(:,:,iq+nmicro)*ptimestep / rat_mmol(iq+nmicro) )
1171	endif
1172	enddo
1173	endif
1174
1175	! i. Condensation of the 3D tracers after the transport
1176	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1177
1178	call calc_ysat(ngrid,nlayer,pplay/100.0,pt,ysat) ! Compute saturation profiles for every grid point (!!p in mbar!!)
1179
1180	dyccond(:,:,:) = 0.D0 ! Default value -> no condensation
1181
1182	do iq=1,nkim
1183	where (ychim(:,:,iq).gt.ysat(:,:,iq))
1184	dyccond(:,:,iq+nmicro) = (-ychim(:,:,iq)+ysat(:,:,iq)) / ptimestep
1185	endwhere
1186	enddo
1187
1188	if (callclouds) then
1189	do iq = 1, size(ices_indx)
1190	dyccond(:,:,gazs_indx(iq)) = 0.D0 ! Condensation have been calculated in the cloud microphysics
1191	enddo
1192	endif
1193
1194	do iq=1,nkim
1195	ychim(:,:,iq) = ychim(:,:,iq) + dyccond(:,:,iq+nmicro)*ptimestep ! update molar ychim for following calchim
1196	pdq(:,:,iq+nmicro) = pdq(:,:,iq+nmicro) + dyccond(:,:,iq+nmicro)*rat_mmol(iq+nmicro) ! convert tendencies to mass mixing ratio
1197	enddo
1198
1199
1200	! ii. 2D zonally averaged fields need to condense and evap before photochem
1201	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1202	if (moyzon_ch .and. mod(icount-1,ichim).eq. 0) then
1203
1204	call calc_ysat(ngrid,nlayer,zplaybar/100.0,ztfibar,ysat) ! Compute saturation profiles for every grid point for the zon-ave fields
1205
1206	dyccondbar(:,:,:) = 0.D0 ! Default value -> no condensation
1207
1208	do iq = 1,nkim
1209	where ( ychimbar(:,:,iq).gt.ysat(:,:,iq) )
1210	dyccondbar(:,:,iq+nmicro) = ( -ychimbar(:,:,iq)+ysat(:,:,iq) ) / ptimestep
1211	endwhere
1212	enddo
1213
1214	if (callclouds) then
1215	do iq = 1, size(ices_indx)
1216	dyccondbar(:,:,gazs_indx(iq)) = 0.D0 ! Condensation have been calculated in the cloud microphysics
1217	enddo
1218	endif
1219
1220	do iq=1,nkim
1221	ychimbar(:,:,iq) = ychimbar(:,:,iq) + dyccondbar(:,:,iq+nmicro)*ptimestep
1222	enddo
1223
1224	endif ! if ( moyzon_ch .and. mod(icount-1,ichim).eq.0 )
1225
1226	! iii. Photochemistry ( must be call after condensation (and evap of 2D) )
1227	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1228	if( mod(icount-1,ichim).eq.0. ) then
1229	print *, "We enter in the photochemistry ..."
1230
1231	if (moyzon_ch) then ! 2D zonally averaged chemistry
1232	! Here we send zonal average fields ( corrected with cond ) from dynamics to chem. module
1233	call calchim(ngrid,ychimbar,declin,ctimestep,ztfibar,zphibar,zphisbar, &
1234	zplaybar,zplevbar,zzlaybar,zzlevbar,dycchi)
1235
1236	else ! 3D chemistry (or 1D run)
1237	call calchim(ngrid,ychim,declin,ctimestep,pt,pphi,phisfi, &
1238	pplay,pplev,zzlay_eff,zzlev_eff,dycchi)
1239	endif ! if moyzon
1240	endif ! if (mod(icount-1,ichim).eq.0)
1241
1242	! iv. Surface pseudo-evaporation
1243	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1244	! Infinite tank of CH4
1245	if (resCH4_inf) then
1246	do ig=1,ngrid
1247	if ((ychim(ig,1,7)+dycchi(ig,1,7)*ptimestep) .lt. botCH4) then ! + dycchi because ychim not yet updated
1248	dycevapCH4(ig) = ( -ychim(ig,1,7)+botCH4 ) / ptimestep - dycchi(ig,1,7)
1249	else
1250	dycevapCH4(ig) = 0.D0
1251	endif
1252	enddo
1253
1254	else
1255	! Fill lakes with precipitation :
1256	!if (REAL(latitude_deg(ig)) .ge. 70.0) then
1257	! tankCH4(ig) = 200.0 ! [m]
1258	!else if (REAL(latitude_deg(ig)) .le. -70.0) then
1259	! tankCH4(ig) = 50.0 ! [m]
1260	!else
1261	! tankCH4(ig) = 0.0 ! [m]
1262	!endif
1263	tankCH4(:) = tankCH4(:) + (mmd_ice_prec(:,1) / 422. * ptimestep) ! [m]
1264
1265	! Evaporation of lakes :
1266	if (moyzon_ch) then
1267	tpq_CH4(:,:) = ychimbar(:,:,7) + dycchi(:,:,7)*ptimestep ! + dycchi because ychim not yet updated [mol/mol]
1268	else
1269	tpq_CH4(:,:) = ychim(:,:,7) + dycchi(:,:,7)*ptimestep ! + dycchi because ychim not yet updated [mol/mol]
1270	endif
1271	call evapCH4(ngrid,nlayer,ptimestep,pplev,zzlay,zzlev,&
1272	pu,pv,tsurf,tpq_CH4,tankCH4,dycevapCH4,zdtsurfevap)
1273	zdtsurf(:) = zdtsurf(:) + zdtsurfevap(:)
1274	endif
1275
1276
1277	pdq(:,1,7+nmicro) = pdq(:,1,7+nmicro) + dycevapCH4(:)*rat_mmol(7+nmicro) ! convert tendencies to mass mixing ratio
1278
1279	! v. Updates and positivity check
1280	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1281	zdqchi(:,:,:) = 0.D0 ! -> dycchi is saved but for the nmicro tracers we must update to 0 at each step
1282
1283	do iq=1,nkim
1284	zdqchi(:,:,iq+nmicro) = dycchi(:,:,iq)*rat_mmol(iq+nmicro) ! convert tendencies to mass mixing ratio
1285
1286	where( (pq(:,:,iq+nmicro) + ( pdq(:,:,iq+nmicro)+zdqchi(:,:,iq+nmicro) )*ptimestep ) .LT. 0.) & ! When using zonal means we set the same tendency
1287	zdqchi(:,:,iq+nmicro) = 1.D-30 - pdq(:,:,iq+nmicro) - pq(:,:,iq+nmicro)/ptimestep ! everywhere in longitude -> could lead to negs !
1288	enddo
1289
1290	pdq(:,:,:) = pdq(:,:,:) + zdqchi(:,:,:)
1291
1292	endif ! end of 'callchim'
1293
1294	! END MESOSCALE
1295	#endif
1296
1297	! ---------------
1298	! V.3 Updates
1299	! ---------------
1300
1301	! Updating Atmospheric Mass and Tracers budgets.
1302	if(mass_redistrib) then
1303
1304	zdmassmr(:,:) = mass(:,:) * zdqevap(:,:)
1305
1306	do ig = 1, ngrid
1307	zdmassmr_col(ig)=SUM(zdmassmr(ig,:))
1308	enddo
1309
1310	call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
1311	call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
1312	call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
1313
1314	call mass_redistribution(ngrid,nlayer,nq,ptimestep, &
1315	capcal,pplay,pplev,pt,tsurf,pq,qsurf, &
1316	pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, &
1317	zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
1318
1319	pdq(:,:,:) = pdq(:,:,:) + zdqmr(:,:,:)
1320	dqsurf(:,:) = dqsurf(:,:) + zdqsurfmr(:,:)
1321	pdt(:,:) = pdt(:,:) + zdtmr(:,:)
1322	pdu(:,:) = pdu(:,:) + zdumr(:,:)
1323	pdv(:,:) = pdv(:,:) + zdvmr(:,:)
1324	pdpsrf(:) = pdpsrf(:) + zdpsrfmr(:)
1325	zdtsurf(:) = zdtsurf(:) + zdtsurfmr(:)
1326	endif
1327
1328	! -----------------------------
1329	! V.4. Surface Tracer Update
1330	! -----------------------------
1331
1332	qsurf(:,:) = qsurf(:,:) + ptimestep*dqsurf(:,:)
1333
1334	! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
1335	! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
1336	qsurf_hist(:,:) = qsurf(:,:)
1337
1338	endif! end of if 'tracer'
1339
1340
1341	!------------------------------------------------
1342	! VI. Surface and sub-surface soil temperature
1343	!------------------------------------------------
1344
1345
1346	! Increment surface temperature
1347
1348	tsurf(:)=tsurf(:)+ptimestep*zdtsurf(:)
1349
1350	! Compute soil temperatures and subsurface heat flux.
1351	if (callsoil) then
1352	call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, &
1353	ptimestep,tsurf,tsoil,capcal,fluxgrd)
1354	endif
1355
1356
1357	! Test energy conservation
1358	if(enertest)then
1359	call planetwide_sumval(cell_area(:)capcal(:)zdtsurf(:)/totarea_planet,dEtots)
1360	if (is_master) print*,'Surface energy change =',dEtots,' W m-2'
1361	endif
1362
1363
1364	!---------------------------------------------------
1365	! VII. Perform diagnostics and write output files
1366	!---------------------------------------------------
1367
1368	! Nudging of zonal wind !
1369	! ~~~~~~~~~~~~~~~~~~~~~~~
1370	if (nudging_u) then
1371	zdundg(:,:) = 0.D0
1372	j=1
1373	! boucle sur les points de grille :
1374	do i = 1, ngrid
1375	! interpolation linéaire des données dans le fichier lu (sur 49 latitudes)
1376	do while ((u_ref(j,1).ge.latitude_deg(i)).and.(j.lt.49))
1377	j=j+1
1378	enddo
1379	factlat = (latitude_deg(i)-u_ref(j,1))/(u_ref(j-1,1)-u_ref(j,1))
1380	! Nudging of the first 23 layers only !!
1381	! IF CHANGE IN VERTICAL RESOLUTION IN THE FIRST 23 LEVELS, IT DOES NOT WORK !!!
1382	zundg(1:23) = factlatu_ref(j-1,2:24)+(1-factlat)u_ref(j,2:24)
1383	zdundg(i,1:23) = (zundg(1:23) - (pu(i,1:23)+pdu(i,1:23)*ptimestep)) / nudging_dt
1384	enddo
1385
1386	pdu(:,:) = pdu(:,:) + zdundg(:,:)
1387	endif
1388
1389	! Note : For output only: the actual model integration is performed in the dynamics.
1390
1391	! Temperature, zonal and meridional winds.
1392	zt(:,:) = pt(:,:) + pdt(:,:)*ptimestep
1393	zu(:,:) = pu(:,:) + pdu(:,:)*ptimestep
1394	zv(:,:) = pv(:,:) + pdv(:,:)*ptimestep
1395
1396	! Diagnostic.
1397	zdtdyn(:,:) = (pt(:,:)-ztprevious(:,:)) / ptimestep
1398	ztprevious(:,:) = zt(:,:)
1399
1400	zdudyn(:,:) = (pu(:,:)-zuprevious(:,:)) / ptimestep
1401	zuprevious(:,:) = zu(:,:)
1402
1403	if(firstcall)then
1404	zdtdyn(:,:)=0.D0
1405	zdudyn(:,:)=0.D0
1406	endif
1407
1408	! Horizotal wind
1409	zhorizwind(:,:) = sqrt( zu(:,:)zu(:,:) + zv(:,:)zv(:,:) )
1410
1411	! Dynamical heating diagnostic.
1412	do ig=1,ngrid
1413	fluxdyn(ig)= SUM(zdtdyn(ig,:) mass(ig,:))cpp
1414	enddo
1415
1416	! [Forcage de la photochimie pour les nuages]
1417	if (callclouds) then
1418	do ig = 1, ngrid
1419	do iq = 1, size(ices_indx)
1420	! C2H2 :
1421	!-------
1422	if(trim(nameOfTracer(gazs_indx(iq))) .eq. "C2H2") then
1423	pdq(ig,nlayer-3:,gazs_indx(iq)) = (4.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
1424	endif
1425	! C2H6 :
1426	!-------
1427	if(trim(nameOfTracer(gazs_indx(iq))) .eq. "C2H6") then
1428	pdq(ig,nlayer-3:,gazs_indx(iq)) = (8.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
1429	endif
1430	! HCN :
1431	!------
1432	if(trim(nameOfTracer(gazs_indx(iq))) .eq. "HCN") then
1433	pdq(ig,nlayer-3:,gazs_indx(iq)) = (2.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
1434	endif
1435	! AC6H6 :
1436	!--------
1437	if(trim(nameOfTracer(gazs_indx(iq))) .eq. "AC6H6") then
1438	pdq(ig,nlayer-3:,gazs_indx(iq)) = (2.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
1439	endif
1440	enddo
1441	enddo
1442	endif
1443
1444	zq(:,:,:) = pq(:,:,:) + pdq(:,:,:)*ptimestep
1445
1446	! Surface pressure.
1447	ps(:) = pplev(:,1) + pdpsrf(:)*ptimestep
1448
1449
1450	! Surface and soil temperature information
1451	call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
1452	call planetwide_minval(tsurf(:),Ts2)
1453	call planetwide_maxval(tsurf(:),Ts3)
1454	if(callsoil)then
1455	TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer
1456	if (is_master) then
1457	print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]'
1458	print*,Ts1,Ts2,Ts3,TsS
1459	end if
1460	else
1461	if (is_master) then
1462	print*,' ave[Tsurf] min[Tsurf] max[Tsurf]'
1463	print*,Ts1,Ts2,Ts3
1464	endif
1465	end if
1466
1467
1468	! Check the energy balance of the simulation during the run
1469	if(corrk)then
1470
1471	call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
1472	call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
1473	call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
1474	call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
1475	call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
1476	do ig=1,ngrid
1477	if(fluxtop_dn(ig).lt.0.0)then
1478	print*,'fluxtop_dn has gone crazy'
1479	print*,'fluxtop_dn=',fluxtop_dn(ig)
1480	print*,'temp= ',pt(ig,:)
1481	print*,'pplay= ',pplay(ig,:)
1482	call abort
1483	endif
1484	end do
1485
1486	if(ngrid.eq.1)then
1487	DYN=0.0
1488	endif
1489
1490	if (is_master) then
1491	print*,' ISR ASR OLR GND DYN [W m^-2]'
1492	print*, ISR,ASR,OLR,GND,DYN
1493	endif
1494
1495	if(enertest .and. is_master)then
1496	print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
1497	print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
1498	print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
1499	endif
1500
1501	if(meanOLR .and. is_master)then
1502	if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
1503	! to record global radiative balance
1504	open(92,file="rad_bal.out",form='formatted',position='append')
1505	write(92,*) zday,ISR,ASR,OLR
1506	close(92)
1507	open(93,file="tem_bal.out",form='formatted',position='append')
1508	if(callsoil)then
1509	write(93,*) zday,Ts1,Ts2,Ts3,TsS
1510	else
1511	write(93,*) zday,Ts1,Ts2,Ts3
1512	endif
1513	close(93)
1514	endif
1515	endif
1516
1517	endif ! end of 'corrk'
1518
1519
1520	! Diagnostic to test radiative-convective timescales in code.
1521	if(testradtimes)then
1522	open(38,file="tau_phys.out",form='formatted',position='append')
1523	ig=1
1524	do l=1,nlayer
1525	write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1526	enddo
1527	close(38)
1528	print*,'As testradtimes enabled,'
1529	print*,'exiting physics on first call'
1530	call abort
1531	endif
1532
1533	if (is_master) print,'--> Ls =',zls180./pi
1534
1535
1536	!----------------------------------------------------------------------
1537	! Writing NetCDF file "RESTARTFI" at the end of the run
1538	!----------------------------------------------------------------------
1539
1540	! Note: 'restartfi' is stored just before dynamics are stored
1541	! in 'restart'. Between now and the writting of 'restart',
1542	! there will have been the itau=itau+1 instruction and
1543	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1544	! thus we store for time=time+dtvr
1545
1546
1547
1548	if(lastcall) then
1549	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1550
1551	#ifndef MESOSCALE
1552	if (ngrid.ne.1) then
1553	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1554
1555	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
1556	ptimestep,ztime_fin, &
1557	tsurf,tsoil,emis,q2,qsurf_hist,tankCH4)
1558	endif
1559	#endif
1560
1561	endif ! end of 'lastcall'
1562
1563	#ifndef MESOSCALE
1564	!-----------------------------------------------------------------------------------------------------
1565	! OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
1566	!
1567	! Note 1 : output with period "ecritphy", set in "run.def"
1568	!
1569	! Note 2 : writediagfi can also be called from any other subroutine for any variable,
1570	! but its preferable to keep all the calls in one place ...
1571	!-----------------------------------------------------------------------------------------------------
1572
1573
1574	call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
1575	call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
1576	call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
1577	call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
1578	call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1579	call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
1580	call writediagfi(ngrid,"temp","temperature","K",3,zt)
1581	call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
1582	call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
1583	call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
1584	call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
1585	call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
1586
1587	! Subsurface temperatures
1588	!call writediagsoil(ngrid,"tempsoil","temperature soil","K",3,tsoil)
1589
1590	! Total energy balance diagnostics
1591	if(callrad.and.(.not.newtonian))then
1592	call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
1593	call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
1594	call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
1595	call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
1596
1597	call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
1598	call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
1599	call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
1600	call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
1601	call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
1602	call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
1603
1604	call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
1605	call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
1606	endif ! end of 'callrad'
1607
1608	if(enertest) then
1609	if (calldifv) then
1610	call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
1611	call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
1612	call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
1613	call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
1614	call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
1615	endif
1616
1617	if (corrk) then
1618	call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
1619	call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
1620	endif
1621	endif ! end of 'enertest'
1622
1623	! Diagnostics of optical thickness
1624	! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
1625	if (diagdtau) then
1626	do nw=1,L_NSPECTV
1627	write(str2,'(i2.2)') nw
1628	call writediagfi(ngrid,'dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',1,int_dtauv(:,nlayer:1:-1,nw))
1629	enddo
1630	do nw=1,L_NSPECTI
1631	write(str2,'(i2.2)') nw
1632	call writediagfi(ngrid,'dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',1,int_dtaui(:,nlayer:1:-1,nw))
1633	enddo
1634	endif
1635
1636	! Temporary inclusions for winds diagnostics.
1637	call writediagfi(ngrid,"zdudif","Turbdiff tend. zon. wind","m s-2",3,zdudif)
1638	call writediagfi(ngrid,"zdudyn","Dyn. tend. zon. wind","m s-2",3,zdudyn)
1639
1640	! Temporary inclusions for heating diagnostics.
1641	call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
1642	call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
1643	call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
1644	call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
1645
1646	! For Debugging.
1647	call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
1648
1649	! Output tracers.
1650	if (tracer) then
1651
1652	if (callmufi) then
1653	! Microphysical tracers are expressed in unit/m3.
1654	! convert X.kg-1 --> X.m-3 (whereas for optics was -> X.m-2)
1655	i2e(:,:) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / gzlat(:,1:nlayer) /(zzlev(:,2:nlayer+1)-zzlev(:,1:nlayer))
1656
1657	#ifdef USE_QTEST
1658	! Microphysical tracers passed through dyn+phys(except mufi)
1659	call writediagfi(ngrid,"mu_m0as_dp","Dynphys only spherical mode 0th order moment",'m-3',3,zq(:,:,micro_indx(1))*i2e)
1660	call writediagfi(ngrid,"mu_m3as_dp","Dynphys only spherical mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(2))*i2e)
1661	call writediagfi(ngrid,"mu_m0af_dp","Dynphys only fractal mode 0th order moment",'m-3',3,zq(:,:,micro_indx(3))*i2e)
1662	call writediagfi(ngrid,"mu_m3af_dp","Dynphys only fractal mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(4))*i2e)
1663	! Microphysical tracers passed through mufi only
1664	call writediagfi(ngrid,"mu_m0as_mo","Mufi only spherical mode 0th order moment",'m-3',3,tpq(:,:,micro_indx(1))*i2e)
1665	call writediagfi(ngrid,"mu_m3as_mo","Mufi only spherical mode 3rd order moment",'m3/m3',3,tpq(:,:,micro_indx(2))*i2e)
1666	call writediagfi(ngrid,"mu_m0af_mo","Mufi only fractal mode 0th order moment",'m-3',3,tpq(:,:,micro_indx(3))*i2e)
1667	call writediagfi(ngrid,"mu_m3af_mo","Mufi only fractal mode 3rd order moment",'m3/m3',3,tpq(:,:,micro_indx(4))*i2e)
1668	#else
1669	call writediagfi(ngrid,"mu_m0as","Spherical mode 0th order moment",'m-3',3,zq(:,:,micro_indx(1))*i2e)
1670	call writediagfi(ngrid,"mu_m3as","Spherical mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(2))*i2e)
1671	call writediagfi(ngrid,"mu_m0af","Fractal mode 0th order moment",'m-3',3,zq(:,:,micro_indx(3))*i2e)
1672	call writediagfi(ngrid,"mu_m3af","Fractal mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(4))*i2e)
1673	#endif
1674
1675	! Microphysical diagnostics
1676	call writediagfi(ngrid,"mmd_aer_prec","Total aerosols precipitations",'m',2,mmd_aer_prec)
1677	call writediagfi(ngrid,"mmd_aer_s_flux","Spherical aerosols sedimentation flux",'kg.m-2.s-1',3,mmd_aer_s_flux)
1678	call writediagfi(ngrid,"mmd_aer_f_flux","Fractal aerosols sedimentation flux",'kg.m-2.s-1',3,mmd_aer_f_flux)
1679	call writediagfi(ngrid,"mmd_rc_sph","Spherical mode caracteristic radius",'m',3,mmd_rc_sph)
1680	call writediagfi(ngrid,"mmd_rc_fra","Fractal mode caracteristic radius",'m',3,mmd_rc_fra)
1681	endif ! end of 'callmufi'
1682
1683	! Chemical tracers
1684	if (callchim) then
1685	do iq=1,nkim
1686	call writediagfi(ngrid,cnames(iq),cnames(iq),'mol/mol',3,zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro))
1687	call writediagfi(ngrid,'dqcond_'//cnames(iq),'dqcond_'//cnames(iq),'mol/mol/s',3,dyccond(:,:,iq+nmicro))
1688	enddo
1689	call writediagfi(ngrid,"evapCH4","Surface CH4 pseudo-evaporation rate",'mol/mol/s',2,dycevapCH4)
1690	endif
1691
1692	endif ! end of 'tracer'
1693
1694	#ifdef CPP_XIOS
1695	!-----------------------------------------------------------------------------------------------------
1696	! XIOS outputs
1697	!-----------------------------------------------------------------------------------------------------
1698	! Send fields to XIOS: (NB these fields must also be defined as
1699	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
1700
1701	!--------------------------------------------------------
1702	! General diagnostics :
1703	!--------------------------------------------------------
1704	CALL send_xios_field("ls",zls*180./pi)
1705	CALL send_xios_field("lss",zlss*180./pi)
1706	CALL send_xios_field("RA",right_ascen*180./pi)
1707	CALL send_xios_field("Declin",declin*180./pi)
1708
1709	! Atmosphere (3D) :
1710	CALL send_xios_field("temp",zt)
1711	CALL send_xios_field("teta",zh)
1712	CALL send_xios_field("p",pplay)
1713	CALL send_xios_field("u",zu)
1714	CALL send_xios_field("v",zv)
1715	CALL send_xios_field("w",pw)
1716
1717	CALL send_xios_field("area",cell_area)
1718	CALL send_xios_field("pphi",pphi)
1719
1720	! Surface (2D) :
1721	CALL send_xios_field("ps",ps)
1722	CALL send_xios_field("tsurf",tsurf)
1723	CALL send_xios_field("pphis",phisfi)
1724
1725	! Total energy balance diagnostics (2D) :
1726	IF (callrad.and.(.not.newtonian)) THEN
1727	CALL send_xios_field("ISR_TOA",fluxtop_dn)
1728	CALL send_xios_field("OLR_TOA",fluxtop_lw)
1729	ENDIF
1730
1731	!--------------------------------------------------------
1732	! Winds trends :
1733	!--------------------------------------------------------
1734	! Atmosphere (3D) :
1735	! du_tot = zdudyn + pdu
1736	CALL send_xios_field("dudyn",zdudyn)
1737	CALL send_xios_field("pdu",pdu)
1738
1739	! pdu = zdudif + zduadj + zdundg
1740	CALL send_xios_field("dudif",zdudif)
1741	CALL send_xios_field("duadj",zduadj)
1742	IF (nudging_u) THEN
1743	CALL send_xios_field("dundg",zdundg)
1744	ENDIF
1745
1746	! zhorizwind = sqrt(uu + vv)
1747	CALL send_xios_field("horizwind",zhorizwind)
1748
1749
1750	!--------------------------------------------------------
1751	! Heating trends :
1752	!--------------------------------------------------------
1753	! Atmosphere (3D) :
1754	! dt_tot = dtdyn + pdt
1755	CALL send_xios_field("dtdyn",zdtdyn)
1756	CALL send_xios_field("pdt",pdt)
1757
1758	! pdt = dtrad + zdtdif + dtadj + zdtlc
1759	CALL send_xios_field("dtrad",dtrad)
1760	CALL send_xios_field("dtdif",zdtdif)
1761	CALL send_xios_field("dtadj",zdtadj(:,:))
1762	IF (callclouds) THEN
1763	CALL send_xios_field("dtlc",zdtlc)
1764	ENDIF
1765
1766	! dtrad = zdtsw + zdtlw
1767	CALL send_xios_field("dtsw",zdtsw)
1768	CALL send_xios_field("dtlw",zdtlw)
1769
1770	! Surface (2D) :
1771	IF(enertest) THEN
1772	IF (calldifv) THEN
1773	CALL send_xios_field("sensibFlux",sensibFlux)
1774	ENDIF
1775	ENDIF
1776	CALL send_xios_field("fluxsurf_lw",fluxsurf_lw(:))
1777	CALL send_xios_field("fluxsurfabs_sw",fluxsurfabs_sw(:))
1778	CALL send_xios_field("emis",emis(:))
1779
1780	! dtsurf = dtsdif + dtsurfevap
1781	CALL send_xios_field("dtsurf",zdtsurf(:))
1782	CALL send_xios_field("dtsdif",zdtsdif(:))
1783	CALL send_xios_field("dtsurfevap",zdtsurfevap(:))
1784
1785
1786	!--------------------------------------------------------
1787	! Optical diagnostics :
1788	!--------------------------------------------------------
1789	! Haze opacity :
1790	CALL send_xios_field('ttauhv_08',zpopthv(:,:,8,2)) ! 08 --> 1.983 um
1791	CALL send_xios_field('ttauhv_15',zpopthv(:,:,15,2)) ! 15 --> 1.000 um
1792	CALL send_xios_field('ttauhv_20',zpopthv(:,:,20,2)) ! 20 --> 0.671 um
1793	CALL send_xios_field('ttauhv_23',zpopthv(:,:,23,2)) ! 23 --> 0.346 um
1794	CALL send_xios_field('ttauhi_02',zpopthi(:,:,2,2)) ! 02 --> 175.3 um
1795	CALL send_xios_field('ttauhi_17',zpopthi(:,:,17,2)) ! 17 --> 11.00 um
1796	CALL send_xios_field('ttauhi_23',zpopthi(:,:,23,2)) ! 23 --> 4.849 um
1797	! Haze extinction :
1798	CALL send_xios_field('kkhv_08',zpopthv(:,:,8,3))
1799	CALL send_xios_field('kkhv_15',zpopthv(:,:,15,3))
1800	CALL send_xios_field('kkhv_20',zpopthv(:,:,20,3))
1801	CALL send_xios_field('kkhv_23',zpopthv(:,:,23,3))
1802	CALL send_xios_field('kkhi_02',zpopthi(:,:,2,3))
1803	CALL send_xios_field('kkhi_17',zpopthi(:,:,17,3))
1804	CALL send_xios_field('kkhi_23',zpopthi(:,:,23,3))
1805	! Haze single scattering albedo :
1806	CALL send_xios_field('wwhv_08',zpopthv(:,:,8,4))
1807	CALL send_xios_field('wwhv_15',zpopthv(:,:,15,4))
1808	CALL send_xios_field('wwhv_20',zpopthv(:,:,20,4))
1809	CALL send_xios_field('wwhv_23',zpopthv(:,:,23,4))
1810	CALL send_xios_field('wwhi_02',zpopthi(:,:,2,4))
1811	CALL send_xios_field('wwhi_17',zpopthi(:,:,17,4))
1812	CALL send_xios_field('wwhi_23',zpopthi(:,:,23,4))
1813	! Haze asymmetry parameter :
1814	CALL send_xios_field('gghv_08',zpopthv(:,:,8,5))
1815	CALL send_xios_field('gghv_15',zpopthv(:,:,15,5))
1816	CALL send_xios_field('gghv_20',zpopthv(:,:,20,5))
1817	CALL send_xios_field('gghv_23',zpopthv(:,:,23,5))
1818	CALL send_xios_field('gghi_02',zpopthi(:,:,2,5))
1819	CALL send_xios_field('gghi_17',zpopthi(:,:,17,5))
1820	CALL send_xios_field('gghi_23',zpopthi(:,:,23,5))
1821
1822	! Diagnostics for haze and clouds :
1823	IF (callclouds) THEN
1824	! Opacity :
1825	CALL send_xios_field('ttauv_08',zpopttv(:,:,8,2)) ! 08 --> 1.983 um
1826	CALL send_xios_field('ttauv_15',zpopttv(:,:,15,2)) ! 15 --> 1.000 um
1827	CALL send_xios_field('ttauv_20',zpopttv(:,:,20,2)) ! 20 --> 0.671 um
1828	CALL send_xios_field('ttauv_23',zpopttv(:,:,23,2)) ! 23 --> 0.346 um
1829	CALL send_xios_field('ttaui_02',zpoptti(:,:,2,2)) ! 02 --> 175.3 um
1830	CALL send_xios_field('ttaui_17',zpoptti(:,:,17,2)) ! 17 --> 11.00 um
1831	CALL send_xios_field('ttaui_23',zpoptti(:,:,23,2)) ! 23 --> 4.849 um
1832	! Extinction :
1833	CALL send_xios_field('kkv_08',zpopttv(:,:,8,3))
1834	CALL send_xios_field('kkv_15',zpopttv(:,:,15,3))
1835	CALL send_xios_field('kkv_20',zpopttv(:,:,20,3))
1836	CALL send_xios_field('kkv_23',zpopttv(:,:,23,3))
1837	CALL send_xios_field('kki_02',zpoptti(:,:,2,3))
1838	CALL send_xios_field('kki_17',zpoptti(:,:,17,3))
1839	CALL send_xios_field('kki_23',zpoptti(:,:,23,3))
1840	! Single scattering albedo :
1841	CALL send_xios_field('wwv_08',zpopttv(:,:,8,4))
1842	CALL send_xios_field('wwv_15',zpopttv(:,:,15,4))
1843	CALL send_xios_field('wwv_20',zpopttv(:,:,20,4))
1844	CALL send_xios_field('wwv_23',zpopttv(:,:,23,4))
1845	CALL send_xios_field('wwi_02',zpoptti(:,:,2,4))
1846	CALL send_xios_field('wwi_17',zpoptti(:,:,17,4))
1847	CALL send_xios_field('wwi_23',zpoptti(:,:,23,4))
1848	! Asymmetry parameter :
1849	CALL send_xios_field('ggv_08',zpopttv(:,:,8,5))
1850	CALL send_xios_field('ggv_15',zpopttv(:,:,15,5))
1851	CALL send_xios_field('ggv_20',zpopttv(:,:,20,5))
1852	CALL send_xios_field('ggv_23',zpopttv(:,:,23,5))
1853	CALL send_xios_field('ggi_02',zpoptti(:,:,2,5))
1854	CALL send_xios_field('ggi_17',zpoptti(:,:,17,5))
1855	CALL send_xios_field('ggi_23',zpoptti(:,:,23,5))
1856	! DRAYAER, TAUGAS, DCONT :
1857	CALL send_xios_field('drayaerv_20',zpopttv(:,:,20,6)) ! 20 --> 0.671um
1858	CALL send_xios_field('taugasv_20',zpopttv(:,:,20,7))
1859	CALL send_xios_field('dcontv_20',zpopttv(:,:,20,8))
1860	CALL send_xios_field('drayaeri_17',zpoptti(:,:,17,6)) ! 17 --> 11.00um
1861	CALL send_xios_field('taugasi_17',zpoptti(:,:,17,7))
1862	CALL send_xios_field('dconti_17',zpoptti(:,:,17,8))
1863	ENDIF
1864
1865	! Diagnostics for haze and clouds (4D) :
1866	CALL send_xios_field('dtauhi',zpopthi(:,:,:,1))
1867	CALL send_xios_field('tauhi',zpopthi(:,:,:,2))
1868	CALL send_xios_field('khi',zpopthi(:,:,:,3))
1869	CALL send_xios_field('whi',zpopthi(:,:,:,4))
1870	CALL send_xios_field('ghi',zpopthi(:,:,:,5))
1871	CALL send_xios_field('dtauhv',zpopthv(:,:,:,1))
1872	CALL send_xios_field('tauhv',zpopthv(:,:,:,2))
1873	CALL send_xios_field('khv',zpopthv(:,:,:,3))
1874	CALL send_xios_field('whv',zpopthv(:,:,:,4))
1875	CALL send_xios_field('ghv',zpopthv(:,:,:,5))
1876	IF (callclouds) THEN
1877	CALL send_xios_field('dtaui',zpoptti(:,:,:,1))
1878	CALL send_xios_field('taui',zpoptti(:,:,:,2))
1879	CALL send_xios_field('ki',zpoptti(:,:,:,3))
1880	CALL send_xios_field('wi',zpoptti(:,:,:,4))
1881	CALL send_xios_field('gi',zpoptti(:,:,:,5))
1882	CALL send_xios_field('dtauv',zpopttv(:,:,:,1))
1883	CALL send_xios_field('tauv',zpopttv(:,:,:,2))
1884	CALL send_xios_field('kv',zpopttv(:,:,:,3))
1885	CALL send_xios_field('wv',zpopttv(:,:,:,4))
1886	CALL send_xios_field('gv',zpopttv(:,:,:,5))
1887	ENDIF
1888
1889	!--------------------------------------------------------
1890	! Microphysical tracers :
1891	!--------------------------------------------------------
1892	IF (callmufi) THEN
1893	! Atmosphere (3D) :
1894	! Moments M0 and M3 :
1895	CALL send_xios_field("mu_m0as",zq(:,:,micro_indx(1))*i2e)
1896	CALL send_xios_field("mu_m3as",zq(:,:,micro_indx(2))*i2e)
1897	CALL send_xios_field("mu_m0af",zq(:,:,micro_indx(3))*i2e)
1898	CALL send_xios_field("mu_m3af",zq(:,:,micro_indx(4))*i2e)
1899	IF (callclouds) THEN
1900	CALL send_xios_field("mu_m0n",zq(:,:,micro_indx(5))*i2e)
1901	CALL send_xios_field("mu_m3n",zq(:,:,micro_indx(6))*i2e)
1902	DO iq = 1, size(ices_indx)
1903	CALL send_xios_field(TRIM(nameOfTracer(ices_indx(iq))),zq(:,:,ices_indx(iq))*i2e)
1904	ENDDO
1905	ENDIF
1906
1907	! Microphysical diagnostics :
1908	CALL send_xios_field("rc_sph",mmd_rc_sph(:,:))
1909	CALL send_xios_field("rc_fra",mmd_rc_fra(:,:))
1910	CALL send_xios_field("vsed_aers",mmd_aer_s_w(:,:))
1911	CALL send_xios_field("vsed_aerf",mmd_aer_f_w(:,:))
1912	CALL send_xios_field("flux_aers",mmd_aer_s_flux(:,:))
1913	CALL send_xios_field("flux_aerf",mmd_aer_f_flux(:,:))
1914	IF (callclouds) THEN
1915	CALL send_xios_field("rc_cld",mmd_rc_cld(:,:))
1916	CALL send_xios_field("vsed_ccn",mmd_ccn_w(:,:))
1917	CALL send_xios_field("flux_ccn",mmd_ccn_flux(:,:))
1918	DO iq = 1, size(ices_indx)
1919	CALL send_xios_field('flux_i'//TRIM(nameOfTracer(gazs_indx(iq))),mmd_ice_fluxes(:,:,iq))
1920	CALL send_xios_field(TRIM(nameOfTracer(gazs_indx(iq)))//'_sat',mmd_gazs_sat(:,:,iq))
1921	ENDDO
1922	ENDIF
1923
1924	! Surface (2D) :
1925	CALL send_xios_field("aer_prec",mmd_aer_prec(:))
1926	IF (callclouds) THEN
1927	CALL send_xios_field("ccn_prec",mmd_ccn_prec(:))
1928	DO iq = 1, size(ices_indx)
1929	CALL send_xios_field('i'//TRIM(nameOfTracer(gazs_indx(iq)))//'_prec',mmd_ice_prec(:,iq))
1930	ENDDO
1931	ENDIF
1932	ENDIF ! of 'if callmufi'
1933
1934	!--------------------------------------------------------
1935	! Chemical tracers :
1936	!--------------------------------------------------------
1937	IF (callchim) THEN
1938	! Surface (2D) :
1939	CALL send_xios_field("evapCH4",dycevapCH4(:)) ! Pseudo-evaporation flux (mol/mol/s)
1940	CALL send_xios_field("tankCH4",tankCH4(:)) ! CH4 tank at the surface (m)
1941
1942	! Atmosphere (3D) :
1943	! Chemical species :
1944	DO iq = 1, nkim
1945	! If no cloud : gzs_indx uninitialized
1946	CALL send_xios_field(trim(cnames(iq)),zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro)) ! kg/kg -> mol/mol
1947	ENDDO
1948
1949	! Condensation tendencies from microphysics (mol/mol/s) :
1950	IF (callclouds) THEN
1951	DO iq = 1, size(ices_indx)
1952	CALL send_xios_field('dmuficond_'//trim(nameOfTracer(gazs_indx(iq))),dmuficond(:,:,iq)/rat_mmol(gazs_indx(iq))) ! kg/kg/s -> mol/mol/s
1953	ENDDO
1954	ENDIF
1955
1956	! Condensation tendencies (mol/mol/s) :
1957	CALL send_xios_field("dqcond_CH4",dyccond(:,:,7+nmicro))
1958	CALL send_xios_field("dqcond_C2H2",dyccond(:,:,10+nmicro))
1959	CALL send_xios_field("dqcond_C2H4",dyccond(:,:,12+nmicro))
1960	CALL send_xios_field("dqcond_C2H6",dyccond(:,:,14+nmicro))
1961	CALL send_xios_field("dqcond_C3H6",dyccond(:,:,17+nmicro))
1962	CALL send_xios_field("dqcond_C4H4",dyccond(:,:,21+nmicro))
1963	CALL send_xios_field("dqcond_CH3CCH",dyccond(:,:,23+nmicro))
1964	CALL send_xios_field("dqcond_C3H8",dyccond(:,:,24+nmicro))
1965	CALL send_xios_field("dqcond_C4H2",dyccond(:,:,25+nmicro))
1966	CALL send_xios_field("dqcond_C4H6",dyccond(:,:,26+nmicro))
1967	CALL send_xios_field("dqcond_C4H10",dyccond(:,:,27+nmicro))
1968	CALL send_xios_field("dqcond_AC6H6",dyccond(:,:,28+nmicro))
1969	CALL send_xios_field("dqcond_HCN",dyccond(:,:,35+nmicro))
1970	CALL send_xios_field("dqcond_CH3CN",dyccond(:,:,39+nmicro))
1971	CALL send_xios_field("dqcond_HC3N",dyccond(:,:,41+nmicro))
1972	CALL send_xios_field("dqcond_NCCN",dyccond(:,:,42+nmicro))
1973	CALL send_xios_field("dqcond_C4N2",dyccond(:,:,43+nmicro))
1974
1975	! Upper atmosphere chemistry variables (3D) :
1976	DO iq = 1, nkim
1977	CALL send_xios_field(trim(cnames(iq))//"_up",ykim_up(iq,:,:)) ! mol/mol
1978	ENDDO
1979
1980	! Total atmosphere chemistry variables (3D) :
1981	! Append fields in ykim_tot for output on the total vertical grid (0->1300km)
1982	DO iq = 1, nkim
1983	! GCM levels
1984	DO l = 1, nlayer
1985	ykim_tot(iq,:,l) = zq(:,l,iq+nmicro)/rat_mmol(iq+nmicro)
1986	ENDDO
1987	! Upper levels
1988	DO l = 1, nlaykim_up
1989	ykim_tot(iq,:,nlayer+l) = ykim_up(iq,:,l)
1990	ENDDO
1991	CALL send_xios_field(trim(cnames(iq))//"_tot",ykim_tot(iq,:,:)) ! mol/mol
1992	ENDDO
1993	ENDIF ! of 'if callchim'
1994
1995
1996	if (lastcall.and.is_omp_master) then
1997	write(,) "physiq: call xios_context_finalize"
1998	call xios_context_finalize
1999	endif
2000	#endif
2001	#else
2002	!MESOSCALE outputs
2003	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
2004	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
2005	comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid)
2006	comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid)
2007	comm_FLUXTOP_LW(1:ngrid)=fluxtop_lw(1:ngrid)
2008	comm_FLUXSURF_SW(1:ngrid)=fluxsurf_sw(1:ngrid)
2009	comm_FLUXSURF_LW(1:ngrid)=fluxsurf_lw(1:ngrid)
2010	comm_FLXGRD(1:ngrid)=fluxgrd(1:ngrid)
2011	sensibFlux(1:ngrid) = zflubid(1:ngrid) - capcal(1:ngrid)*zdtsdif(1:ngrid)
2012	#endif
2013
2014	icount=icount+1
2015
2016	end subroutine physiq
2017
2018	end module physiq_mod

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

Download in other formats:

Original Format