Context Navigation

physiq_mod.F90 @ 1790

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

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

Download in other formats:

Original Format