Context Navigation

physiq_mod.F90 @ 3709

Last change on this file since 3709 was 3682, checked in by lrosset, 8 months ago

Titan microphysics and physics : Outputs the nucleation rates and growth rates for the condensible species. (LR)

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

File size: 86.7 KB

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

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

Download in other formats:

Original Format