Context Navigation

physiq_mod.F90 @ 2239

Last change on this file since 2239 was 2138, checked in by jvatant, 6 years ago

Correct outptuts of dtaui/v with correct ponderations
of weights, but now it is attenuation exp(-tau)
--JVO

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

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

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

Download in other formats:

Original Format