Context Navigation

physiq_mod.F90 @ 3026

Last change on this file since 3026 was 2742, checked in by aslmd, 2 years ago

small bug fix on the name of in WRF planet_type check callphysi_mod for Titan

MESOSCALE interface: modifications of local time management for idealized studies (LES) of Titan. removed undefined variables pday ptime.

MESOSCALE. Titan. bug fix: added a test to set surface roughness to a default value in case it is not set thus zero. moved the setting of Z0 outside the loop since it is not (yet?) varying with space in the Titan model

MESOSCALE. Titan. adding a comment for future work in turb_mod

added missing implicit none and kps kpe (found by comparing with updata_outputs from Mars). a couple cosmetic spaces to have something closer to Mars for comparisons.

added diagnostics in update_outputs that were sent to comm_wrf module but not invoked in update_outputs to be available for I/O

alternate computation of sensible heat flux in Titan LES. equivalent to Mars. results are strictly similar to what could be obtained within vdifc or turbdiff computations.

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

File size: 69.7 KB

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

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

Download in other formats:

Original Format