Context Navigation

physiq_mod.F90 @ 1722

Last change on this file since 1722 was 1722, checked in by jvatant, 7 years ago
Adapt various modifs of LMDZ.GENERIC to LMDZ.TITAN from r1690-1694-1699-1709-1715 --JVO
Property svn:executable set to ``*
File size: 62.7 KB

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

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

Download in other formats:

Original Format