Context Navigation

physiq_mod.F90 @ 3340

Last change on this file since 3340 was 3340, checked in by slebonnois, 6 months ago

SL: Lucie's modifications to run without clouds and still get good tropospheric T profile

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

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

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

Download in other formats:

Original Format