Context Navigation

radiation_save.F90 @ 3981

Last change on this file since 3981 was 3908, checked in by idelkadi, 3 years ago

Online implementation of the radiative transfer code ECRAD in the LMDZ model.

Inclusion of the ecrad directory containing the sources of the ECRAD code
- interface routine : radiation_scheme.F90
Adaptation of compilation scripts :
- compilation under CPP key CPP_ECRAD
- compilation with option "-rad ecard" or "-ecard true"
- The "-rad old/rtm/ecran" build option will need to replace the "-rrtm true" and "-ecrad true" options in the future.
Runing LMDZ simulations with ecrad, you need :
- logical key iflag_rrtm = 2 in physiq.def
- namelist_ecrad (DefLists?)
- the directory "data" containing the configuration files is temporarily placed in ../libfphylmd/ecrad/
Compilation and execution are tested in the 1D case. The repository under svn would allow to continue the implementation work: tests, verification of the results, ...

File size: 43.0 KB

Line
1	! radiation_save.F90 - Save data to NetCDF files
2	!
3	! (C) Copyright 2014- ECMWF.
4	!
5	! This software is licensed under the terms of the Apache Licence Version 2.0
6	! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
7	!
8	! In applying this licence, ECMWF does not waive the privileges and immunities
9	! granted to it by virtue of its status as an intergovernmental organisation
10	! nor does it submit to any jurisdiction.
11	!
12	! Author: Robin Hogan
13	! Email: r.j.hogan@ecmwf.int
14	!
15	! Modifications
16	! 2017-04-22 R. Hogan Adapt for new way of describing longwave properties
17	! 2019-01-02 R. Hogan Only save cloud properties if do_clouds==.true.
18
19	module radiation_save
20
21	use parkind1, only : jprb
22
23	implicit none
24
25	! Two available subroutines: save final fluxes and save intermediate
26	! radiative properties
27	public :: save_fluxes, save_radiative_properties, save_inputs
28
29	contains
30
31	!---------------------------------------------------------------------
32	! Save fluxes in "flux" to NetCDF file_name, plus pressure from the
33	! thermodynamics object
34	subroutine save_fluxes(file_name, config, thermodynamics, flux, &
35	& iverbose, is_hdf5_file, experiment_name)
36
37	use yomhook, only : lhook, dr_hook
38
39	use easy_netcdf
40
41	use radiation_io, only : nulout
42	use radiation_config, only : config_type, IGasModelMonochromatic
43	use radiation_thermodynamics, only : thermodynamics_type
44	use radiation_flux, only : flux_type
45
46	character(len=*), intent(in) :: file_name
47	type(config_type), intent(in) :: config
48	type(thermodynamics_type), intent(in) :: thermodynamics
49	type(flux_type), intent(in) :: flux
50	integer, optional, intent(in) :: iverbose
51	logical, optional, intent(in) :: is_hdf5_file
52	character(len=*), optional, intent(in) :: experiment_name
53
54	type(netcdf_file) :: out_file
55	integer :: ncol, n_lev_plus1
56	character(5), parameter :: default_lw_units_str = 'W m-2'
57	character(5) :: lw_units_str
58	integer :: i_local_verbose
59
60	real(jprb) :: hook_handle
61
62	if (lhook) call dr_hook('radiation_save:save_fluxes',0,hook_handle)
63
64	if (present(iverbose)) then
65	i_local_verbose = iverbose
66	else
67	i_local_verbose = config%iverbose
68	end if
69
70	! Work out array dimensions
71	if (config%do_sw) then
72	ncol = size(flux%sw_up,1)
73	n_lev_plus1 = size(flux%sw_up,2)
74	elseif (config%do_lw) then
75	ncol = size(flux%lw_up,1)
76	n_lev_plus1 = size(flux%lw_up,2)
77	else
78	if (i_local_verbose >= 1) then
79	write(nulout,'(a,a,a)') 'Warning: neither longwave nor shortwave computed so ', &
80	& file_name,' not written'
81	end if
82	return
83	end if
84
85	if (config%i_gas_model == IGasModelMonochromatic &
86	.and. config%mono_lw_wavelength > 0.0_jprb) then
87	lw_units_str = 'W m-3'
88	else
89	lw_units_str = default_lw_units_str
90	end if
91
92	! Open the file
93	call out_file%create(trim(file_name), iverbose=i_local_verbose, is_hdf5_file=is_hdf5_file)
94
95	! Variables stored internally with column varying fastest, but in
96	! output file column varies most slowly so need to transpose
97	call out_file%transpose_matrices(.true.)
98
99	! Spectral fluxes in memory are dimensioned (nband,ncol,nlev), but
100	! are reoriented in the output file to be (nband,nlev,ncol), where
101	! the convention here is first dimension varying fastest
102	call out_file%permute_3d_arrays( (/ 1, 3, 2 /) )
103
104	! Define dimensions
105	call out_file%define_dimension("column", ncol)
106	call out_file%define_dimension("half_level", n_lev_plus1)
107
108	if (config%do_save_spectral_flux) then
109	if (config%do_lw) then
110	call out_file%define_dimension("band_lw", config%n_spec_lw)
111	end if
112	if (config%do_sw) then
113	call out_file%define_dimension("band_sw", config%n_spec_sw)
114	end if
115	else if (config%do_surface_sw_spectral_flux) then
116	if (config%do_sw) then
117	call out_file%define_dimension("band_sw", config%n_bands_sw)
118	end if
119	end if
120
121	if (config%do_lw .and. config%do_canopy_fluxes_lw) then
122	call out_file%define_dimension("canopy_band_lw", &
123	& size(flux%lw_dn_surf_canopy, 1))
124	end if
125	if (config%do_sw .and. config%do_canopy_fluxes_sw) then
126	call out_file%define_dimension("canopy_band_sw", &
127	& size(flux%sw_dn_diffuse_surf_canopy, 1))
128	end if
129
130	! Put global attributes
131	call out_file%put_global_attributes( &
132	& title_str="Radiative flux profiles from the ecRad offline radiation model", &
133	& references_str="Hogan, R. J., and A. Bozzo, 2018: A flexible and efficient radiation " &
134	& //"scheme for the ECMWF model. J. Adv. Modeling Earth Sys., 10, 1990–2008", &
135	& source_str="ecRad offline radiation model")
136
137	! Save "experiment" global attribute if present and not empty
138	if (present(experiment_name)) then
139	if (experiment_name /= " ") then
140	call out_file%put_global_attribute("experiment", experiment_name)
141	end if
142	end if
143
144	! Define variables
145	call out_file%define_variable("pressure_hl", &
146	& dim2_name="column", dim1_name="half_level", &
147	& units_str="Pa", long_name="Pressure", &
148	& standard_name="air_pressure")
149
150	if (config%do_lw) then
151	call out_file%define_variable("flux_up_lw", &
152	& dim2_name="column", dim1_name="half_level", &
153	& units_str=lw_units_str, long_name="Upwelling longwave flux", &
154	& standard_name="upwelling_longwave_flux_in_air")
155	call out_file%define_variable("flux_dn_lw", &
156	& dim2_name="column", dim1_name="half_level", &
157	& units_str=lw_units_str, long_name="Downwelling longwave flux", &
158	& standard_name="downwelling_longwave_flux_in_air")
159	if (config%do_clear) then
160	call out_file%define_variable("flux_up_lw_clear", &
161	& dim2_name="column", dim1_name="half_level", &
162	& units_str=lw_units_str, &
163	& long_name="Upwelling clear-sky longwave flux")
164	call out_file%define_variable("flux_dn_lw_clear", &
165	& dim2_name="column", dim1_name="half_level", &
166	& units_str=lw_units_str, &
167	& long_name="Downwelling clear-sky longwave flux")
168	end if
169
170	if (config%do_lw_derivatives) then
171	call out_file%define_variable("lw_derivative", &
172	& dim2_name="column", dim1_name="half_level", &
173	& units_str="1", &
174	& long_name="Derivative of upwelling LW flux w.r.t. surface value")
175	end if
176
177	if (config%do_save_spectral_flux) then
178	call out_file%define_variable("spectral_flux_up_lw", &
179	& dim3_name="column", dim2_name="half_level", &
180	& dim1_name="band_lw", units_str=lw_units_str, &
181	& long_name="Spectral upwelling longwave flux")
182	call out_file%define_variable("spectral_flux_dn_lw", &
183	& dim3_name="column", dim2_name="half_level", &
184	& dim1_name="band_lw", units_str=lw_units_str, &
185	& long_name="Spectral downwelling longwave flux")
186	if (config%do_clear) then
187	call out_file%define_variable("spectral_flux_up_lw_clear", &
188	& dim3_name="column", dim2_name="half_level", &
189	& dim1_name="band_lw", units_str=lw_units_str, &
190	& long_name="Spectral upwelling clear-sky longwave flux")
191	call out_file%define_variable("spectral_flux_dn_lw_clear", &
192	& dim3_name="column", dim2_name="half_level", &
193	& dim1_name="band_lw", units_str=lw_units_str, &
194	& long_name="Spectral downwelling clear-sky longwave flux")
195	end if
196	end if
197
198	if (config%do_canopy_fluxes_lw) then
199	call out_file%define_variable("canopy_flux_dn_lw_surf", &
200	& dim2_name="column", dim1_name="canopy_band_lw", units_str=lw_units_str, &
201	& long_name="Surface downwelling longwave flux in canopy bands")
202	end if
203
204	end if
205
206	if (config%do_sw) then
207	call out_file%define_variable("flux_up_sw", &
208	& dim2_name="column", dim1_name="half_level", &
209	& units_str="W m-2", long_name="Upwelling shortwave flux", &
210	& standard_name="upwelling_shortwave_flux_in_air")
211	call out_file%define_variable("flux_dn_sw", &
212	& dim2_name="column", dim1_name="half_level", &
213	& units_str="W m-2", long_name="Downwelling shortwave flux", &
214	& standard_name="downwelling_shortwave_flux_in_air")
215	if (config%do_sw_direct) then
216	call out_file%define_variable("flux_dn_direct_sw", &
217	& dim2_name="column", dim1_name="half_level", &
218	& units_str="W m-2", &
219	& long_name="Downwelling direct shortwave flux")
220	end if
221	if (config%do_clear) then
222	call out_file%define_variable("flux_up_sw_clear", &
223	& dim2_name="column", dim1_name="half_level", &
224	& units_str="W m-2", &
225	& long_name="Upwelling clear-sky shortwave flux")
226	call out_file%define_variable("flux_dn_sw_clear", &
227	& dim2_name="column", dim1_name="half_level", &
228	& units_str="W m-2", &
229	& long_name="Downwelling clear-sky shortwave flux")
230	if (config%do_sw_direct) then
231	call out_file%define_variable("flux_dn_direct_sw_clear", &
232	& dim2_name="column", dim1_name="half_level", &
233	& units_str="W m-2", &
234	& long_name="Downwelling clear-sky direct shortwave flux")
235	end if
236	end if
237
238	if (config%do_save_spectral_flux) then
239	call out_file%define_variable("spectral_flux_up_sw", &
240	& dim3_name="column", dim2_name="half_level", &
241	& dim1_name="band_sw", units_str="W m-2", &
242	& long_name="Spectral upwelling shortwave flux")
243	call out_file%define_variable("spectral_flux_dn_sw", &
244	& dim3_name="column", dim2_name="half_level", &
245	& dim1_name="band_sw", units_str="W m-2", &
246	& long_name="Spectral downwelling shortwave flux")
247	if (config%do_sw_direct) then
248	call out_file%define_variable("spectral_flux_dn_direct_sw", &
249	& dim3_name="column", dim2_name="half_level", &
250	& dim1_name="band_sw", units_str="W m-2", &
251	& long_name="Spectral downwelling direct shortwave flux")
252	end if
253	if (config%do_clear) then
254	call out_file%define_variable("spectral_flux_up_sw_clear", &
255	& dim3_name="column", dim2_name="half_level", &
256	& dim1_name="band_sw", units_str="W m-2", &
257	& long_name="Spectral upwelling clear-sky shortwave flux")
258	call out_file%define_variable("spectral_flux_dn_sw_clear", &
259	& dim3_name="column", dim2_name="half_level", &
260	& dim1_name="band_sw", units_str="W m-2", &
261	& long_name="Spectral downwelling clear-sky shortwave flux")
262	if (config%do_sw_direct) then
263	call out_file%define_variable("spectral_flux_dn_direct_sw_clear", &
264	& dim3_name="column", dim2_name="half_level", &
265	& dim1_name="band_sw", units_str="W m-2", &
266	& long_name="Spectral downwelling clear-sky direct shortwave flux")
267	end if
268	end if
269	else if (config%do_surface_sw_spectral_flux) then
270	call out_file%define_variable("spectral_flux_dn_sw_surf", &
271	& dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
272	& long_name="Spectral downwelling shortwave flux at surface")
273	call out_file%define_variable("spectral_flux_dn_direct_sw_surf", &
274	& dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
275	& long_name="Spectral downwelling direct shortwave flux at surface")
276	if (config%do_clear) then
277	call out_file%define_variable("spectral_flux_dn_sw_surf_clear", &
278	& dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
279	& long_name="Spectral downwelling clear-sky shortwave flux at surface")
280	call out_file%define_variable("spectral_flux_dn_direct_sw_surf_clear", &
281	& dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
282	& long_name="Spectral downwelling clear-sky direct shortwave flux at surface")
283	end if
284	end if
285
286	if (config%do_canopy_fluxes_sw) then
287	call out_file%define_variable("canopy_flux_dn_diffuse_sw_surf", &
288	& dim2_name="column", dim1_name="canopy_band_sw", units_str="W m-2", &
289	& long_name="Surface downwelling diffuse shortwave flux in canopy bands")
290	call out_file%define_variable("canopy_flux_dn_direct_sw_surf", &
291	& dim2_name="column", dim1_name="canopy_band_sw", units_str="W m-2", &
292	& long_name="Surface downwelling direct shortwave flux in canopy bands")
293	end if
294
295	end if
296
297	if (config%do_lw .and. config%do_clouds) then
298	call out_file%define_variable("cloud_cover_lw", &
299	& dim1_name="column", units_str="1", &
300	& long_name="Total cloud cover diagnosed by longwave solver", &
301	& standard_name="cloud_area_fraction")
302	end if
303	if (config%do_sw .and. config%do_clouds) then
304	call out_file%define_variable("cloud_cover_sw", &
305	& dim1_name="column", units_str="1", &
306	& long_name="Total cloud cover diagnosed by shortwave solver", &
307	& standard_name="cloud_area_fraction")
308	end if
309
310	! Write variables
311
312	call out_file%put("pressure_hl", thermodynamics%pressure_hl)
313
314	if (config%do_lw) then
315	call out_file%put("flux_up_lw", flux%lw_up)
316	call out_file%put("flux_dn_lw", flux%lw_dn)
317	if (config%do_clear) then
318	call out_file%put("flux_up_lw_clear", flux%lw_up_clear)
319	call out_file%put("flux_dn_lw_clear", flux%lw_dn_clear)
320	end if
321
322	if (config%do_lw_derivatives) then
323	call out_file%put("lw_derivative", flux%lw_derivatives)
324	end if
325
326	if (config%do_save_spectral_flux) then
327	call out_file%put("spectral_flux_up_lw", flux%lw_up_band)
328	call out_file%put("spectral_flux_dn_lw", flux%lw_dn_band)
329	if (config%do_clear) then
330	call out_file%put("spectral_flux_up_lw_clear", flux%lw_up_clear_band)
331	call out_file%put("spectral_flux_dn_lw_clear", flux%lw_dn_clear_band)
332	end if
333	end if
334
335	if (config%do_canopy_fluxes_lw) then
336	call out_file%put("canopy_flux_dn_lw_surf", flux%lw_dn_surf_canopy, &
337	& do_transp = .false.)
338	end if
339
340	end if
341
342	if (config%do_sw) then
343	call out_file%put("flux_up_sw", flux%sw_up)
344	call out_file%put("flux_dn_sw", flux%sw_dn)
345	if (config%do_sw_direct) then
346	call out_file%put("flux_dn_direct_sw", flux%sw_dn_direct)
347	end if
348	if (config%do_clear) then
349	call out_file%put("flux_up_sw_clear", flux%sw_up_clear)
350	call out_file%put("flux_dn_sw_clear", flux%sw_dn_clear)
351	if (config%do_sw_direct) then
352	call out_file%put("flux_dn_direct_sw_clear", flux%sw_dn_direct_clear)
353	end if
354	end if
355
356	if (config%do_save_spectral_flux) then
357	call out_file%put("spectral_flux_up_sw", flux%sw_up_band)
358	call out_file%put("spectral_flux_dn_sw", flux%sw_dn_band)
359	if (config%do_sw_direct) then
360	call out_file%put("spectral_flux_dn_direct_sw", &
361	& flux%sw_dn_direct_band)
362	end if
363	if (config%do_clear) then
364	call out_file%put("spectral_flux_up_sw_clear", flux%sw_up_clear_band)
365	call out_file%put("spectral_flux_dn_sw_clear", flux%sw_dn_clear_band)
366	if (config%do_sw_direct) then
367	call out_file%put("spectral_flux_dn_direct_sw_clear", &
368	& flux%sw_dn_direct_clear_band)
369	end if
370	end if
371	else if (config%do_surface_sw_spectral_flux) then
372	call out_file%put("spectral_flux_dn_sw_surf", flux%sw_dn_surf_band, &
373	& do_transp=.false.)
374	call out_file%put("spectral_flux_dn_direct_sw_surf", flux%sw_dn_direct_surf_band, &
375	& do_transp=.false.)
376	if (config%do_clear) then
377	call out_file%put("spectral_flux_dn_sw_surf_clear", flux%sw_dn_surf_clear_band, &
378	& do_transp=.false.)
379	call out_file%put("spectral_flux_dn_direct_sw_surf_clear", &
380	& flux%sw_dn_direct_surf_clear_band, do_transp=.false.)
381	end if
382	end if
383
384	if (config%do_canopy_fluxes_sw) then
385	call out_file%put("canopy_flux_dn_diffuse_sw_surf", flux%sw_dn_diffuse_surf_canopy, &
386	& do_transp = .false.)
387	call out_file%put("canopy_flux_dn_direct_sw_surf", flux%sw_dn_direct_surf_canopy, &
388	& do_transp = .false.)
389	end if
390
391	end if
392
393	if (config%do_lw .and. config%do_clouds) then
394	call out_file%put("cloud_cover_lw", flux%cloud_cover_lw)
395	end if
396	if (config%do_sw .and. config%do_clouds) then
397	call out_file%put("cloud_cover_sw", flux%cloud_cover_sw)
398	end if
399
400	! Close file
401	call out_file%close()
402
403	if (lhook) call dr_hook('radiation_save:save_fluxes',1,hook_handle)
404
405	end subroutine save_fluxes
406
407
408	!---------------------------------------------------------------------
409	! Save intermediate radiative properties, specifically the
410	! scattering and absorption properties at each g-point/band
411	subroutine save_radiative_properties(file_name, nlev, &
412	& istartcol, iendcol, &
413	& config, single_level, thermodynamics, cloud, &
414	& planck_hl, lw_emission, lw_albedo, &
415	& sw_albedo_direct, sw_albedo_diffuse, &
416	& incoming_sw, &
417	& od_lw, ssa_lw, g_lw, &
418	& od_sw, ssa_sw, g_sw, &
419	& od_lw_cloud, ssa_lw_cloud, g_lw_cloud, &
420	& od_sw_cloud, ssa_sw_cloud, g_sw_cloud)
421
422	use radiation_config, only : config_type
423	use radiation_single_level, only : single_level_type
424	use radiation_thermodynamics,only : thermodynamics_type
425	use radiation_cloud, only : cloud_type
426	use easy_netcdf
427
428	character(len=*), intent(in) :: file_name
429	type(config_type), intent(in) :: config
430	type(single_level_type), intent(in) :: single_level
431	type(thermodynamics_type),intent(in) :: thermodynamics
432	type(cloud_type), intent(in) :: cloud
433
434	integer, intent(in) :: nlev, istartcol, iendcol
435
436	! Input variables, as defined in radiation_interface.F90
437
438	! Layer optical depth, single scattering albedo and asymmetry factor of
439	! gases and aerosols at each shortwave g-point
440	real(jprb), intent(in), dimension(config%n_g_sw,nlev,istartcol:iendcol) :: od_sw, ssa_sw, g_sw
441
442	! Layer optical depth, single scattering albedo and asymmetry factor of
443	! hydrometeors in each shortwave band
444	real(jprb), intent(in), dimension(config%n_bands_sw,nlev,istartcol:iendcol) :: &
445	& od_sw_cloud, ssa_sw_cloud, g_sw_cloud
446
447	! Direct and diffuse surface albedo, and the incoming shortwave
448	! flux into a plane perpendicular to the incoming radiation at
449	! top-of-atmosphere in each of the shortwave g-points
450	real(jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol) &
451	& :: sw_albedo_direct, sw_albedo_diffuse, incoming_sw
452
453	! Layer optical depth, single scattering albedo and asymmetry factor of
454	! gases and aerosols at each longwave g-point, where the latter
455	! two variables are only defined if aerosol longwave scattering is
456	! enabled (otherwise both are treated as zero).
457	real(jprb), intent(in), dimension(config%n_g_lw,nlev,istartcol:iendcol) :: od_lw
458	real(jprb), intent(in), dimension(config%n_g_lw_if_scattering,nlev,istartcol:iendcol) :: &
459	& ssa_lw, g_lw
460
461	! Layer optical depth, single scattering albedo and asymmetry factor of
462	! hydrometeors in each longwave band, where the latter two
463	! variables are only defined if hydrometeor longwave scattering is
464	! enabled (otherwise both are treated as zero).
465	real(jprb), intent(in), dimension(config%n_bands_lw,nlev,istartcol:iendcol) :: od_lw_cloud
466	real(jprb), intent(in), dimension(config%n_bands_lw_if_scattering,nlev,istartcol:iendcol) :: &
467	& ssa_lw_cloud, g_lw_cloud
468
469	! The Planck function (emitted flux from a black body) at half
470	! levels and at the surface at each longwave g-point
471	real(jprb), intent(in), dimension(config%n_g_lw,nlev+1,istartcol:iendcol) :: planck_hl
472
473	! Emission (Planck*emissivity) and albedo (1-emissivity) at the
474	! surface at each longwave g-point
475	real(jprb), dimension(config%n_g_lw, istartcol:iendcol) :: lw_emission, lw_albedo
476
477	integer :: n_col_local ! Number of columns from istartcol to iendcol
478
479	! Object for output NetCDF file
480	type(netcdf_file) :: out_file
481
482	n_col_local = iendcol + 1 - istartcol
483
484	! Alas the NetCDF library is not thread-safe for writing, so we
485	! must write radiative-property files serially
486
487	!$OMP CRITICAL
488
489	! Open the file
490	call out_file%create(trim(file_name), iverbose=config%iverbose)
491
492	! Configure matrix and 3D-array orientation
493	call out_file%transpose_matrices(.true.)
494
495	! Sometimes the Planck function values are very large or small
496	! even if the fluxes are within a manageable range
497	call out_file%double_precision(.true.)
498
499	! Define dimensions
500	! call out_file%define_dimension("column", n_col_local)
501	call out_file%define_dimension("column", 0) ! "Unlimited" dimension
502	call out_file%define_dimension("level", nlev)
503	call out_file%define_dimension("half_level", nlev+1)
504	if (config%do_clouds) then
505	call out_file%define_dimension("level_interface", nlev-1)
506	end if
507
508	if (config%do_lw) then
509	call out_file%define_dimension("gpoint_lw", config%n_g_lw)
510	if (config%do_clouds) then
511	call out_file%define_dimension("band_lw", config%n_bands_lw)
512	end if
513	end if
514	if (config%do_sw) then
515	call out_file%define_dimension("gpoint_sw", config%n_g_sw)
516	if (config%do_clouds) then
517	call out_file%define_dimension("band_sw", config%n_bands_sw)
518	end if
519	end if
520
521	! Put global attributes
522	call out_file%put_global_attributes( &
523	& title_str="Spectral radiative properties from the ecRad offline radiation model", &
524	& references_str="Hogan, R. J., and A. Bozzo, 2018: A flexible and efficient radiation " &
525	& //"scheme for the ECMWF model. J. Adv. Modeling Earth Sys., 10, 1990–2008", &
526	& source_str="ecRad offline radiation model")
527
528	! Define variables
529	call out_file%define_variable("pressure_hl", &
530	& dim2_name="column", dim1_name="half_level", &
531	& units_str="Pa", long_name="Pressure on half-levels")
532
533	if (allocated(thermodynamics%h2o_sat_liq) .and. config%use_aerosols) then
534	call out_file%define_variable("q_sat_liquid", &
535	& dim2_name="column", dim1_name="level", &
536	& units_str="kg kg-1", long_name="Specific humidity at liquid saturation")
537	end if
538
539	if (config%do_sw) then
540	call out_file%define_variable("cos_solar_zenith_angle", &
541	& dim1_name="column", units_str="1", &
542	& long_name="Cosine of the solar zenith angle")
543	end if
544
545	if (config%do_clouds) then
546	call out_file%define_variable("cloud_fraction", &
547	& dim2_name="column", dim1_name="level", &
548	& units_str="1", long_name="Cloud fraction")
549	call out_file%define_variable("overlap_param", &
550	& dim2_name="column", dim1_name="level_interface", &
551	& units_str="1", long_name="Cloud overlap parameter")
552	end if
553
554	if (config%do_lw) then
555	call out_file%define_variable("planck_hl", &
556	& dim3_name="column", dim2_name="half_level", dim1_name="gpoint_lw", &
557	& units_str="W m-2", long_name="Planck function on half-levels")
558	call out_file%define_variable("lw_emission", &
559	& dim2_name="column", dim1_name="gpoint_lw", &
560	& units_str="W m-2", long_name="Longwave surface emission")
561	call out_file%define_variable("lw_emissivity", &
562	& dim2_name="column", dim1_name="gpoint_lw", &
563	& units_str="1", long_name="Surface longwave emissivity")
564
565	call out_file%define_variable("od_lw", &
566	& dim3_name="column", dim2_name="level", dim1_name="gpoint_lw", &
567	& units_str="1", long_name="Clear-sky longwave optical depth")
568	if (config%do_lw_aerosol_scattering) then
569	call out_file%define_variable("ssa_lw", &
570	& dim3_name="column", dim2_name="level", dim1_name="gpoint_lw", &
571	& units_str="1", long_name="Clear-sky longwave single scattering albedo")
572	call out_file%define_variable("asymmetry_lw", &
573	& dim3_name="column", dim2_name="level", dim1_name="gpoint_lw", &
574	& units_str="1", long_name="Clear-sky longwave asymmetry factor")
575	end if
576
577	if (config%do_clouds) then
578	call out_file%define_variable("od_lw_cloud", &
579	& dim3_name="column", dim2_name="level", dim1_name="band_lw", &
580	& units_str="1", long_name="In-cloud longwave optical depth")
581	if (config%do_lw_cloud_scattering) then
582	call out_file%define_variable("ssa_lw_cloud", &
583	& dim3_name="column", dim2_name="level", dim1_name="band_lw", &
584	& units_str="1", long_name="Cloud longwave single scattering albedo")
585	call out_file%define_variable("asymmetry_lw_cloud", &
586	& dim3_name="column", dim2_name="level", dim1_name="band_lw", &
587	& units_str="1", long_name="Cloud longwave asymmetry factor")
588	end if
589	end if ! do_clouds
590	end if ! do_lw
591
592	if (config%do_sw) then
593	call out_file%define_variable("incoming_sw", &
594	& dim2_name="column", dim1_name="gpoint_sw", &
595	& units_str="W m-2", long_name="Incoming shortwave flux at top-of-atmosphere in direction of sun")
596
597	call out_file%define_variable("sw_albedo", &
598	& dim2_name="column", dim1_name="gpoint_sw", &
599	& units_str="1", long_name="Surface shortwave albedo to diffuse radiation")
600	call out_file%define_variable("sw_albedo_direct", &
601	& dim2_name="column", dim1_name="gpoint_sw", &
602	& units_str="1", long_name="Surface shortwave albedo to direct radiation")
603
604	call out_file%define_variable("od_sw", &
605	& dim3_name="column", dim2_name="level", dim1_name="gpoint_sw", &
606	& units_str="1", long_name="Clear-sky shortwave optical depth")
607	call out_file%define_variable("ssa_sw", &
608	& dim3_name="column", dim2_name="level", dim1_name="gpoint_sw", &
609	& units_str="1", long_name="Clear-sky shortwave single scattering albedo")
610	call out_file%define_variable("asymmetry_sw", &
611	& dim3_name="column", dim2_name="level", dim1_name="gpoint_sw", &
612	& units_str="1", long_name="Clear-sky shortwave asymmetry factor")
613
614	if (config%do_clouds) then
615	call out_file%define_variable("od_sw_cloud", &
616	& dim3_name="column", dim2_name="level", dim1_name="band_sw", &
617	& units_str="1", long_name="In-cloud shortwave optical depth")
618	call out_file%define_variable("ssa_sw_cloud", &
619	& dim3_name="column", dim2_name="level", dim1_name="band_sw", &
620	& units_str="1", long_name="Cloud shortwave single scattering albedo")
621	call out_file%define_variable("asymmetry_sw_cloud", &
622	& dim3_name="column", dim2_name="level", dim1_name="band_sw", &
623	& units_str="1", long_name="Cloud shortwave asymmetry factor")
624	end if
625	end if
626
627	if (config%do_clouds) then
628	if (allocated(cloud%fractional_std)) then
629	call out_file%define_variable("fractional_std", &
630	& dim2_name="column", dim1_name="level", units_str="1", &
631	& long_name="Fractional standard deviation of cloud optical depth")
632	end if
633	if (allocated(cloud%inv_cloud_effective_size)) then
634	call out_file%define_variable("inv_cloud_effective_size", &
635	& dim2_name="column", dim1_name="level", units_str="m-1", &
636	& long_name="Inverse of cloud effective horizontal size")
637	end if
638	if (allocated(cloud%inv_inhom_effective_size)) then
639	call out_file%define_variable("inv_inhom_effective_size", &
640	& dim2_name="column", dim1_name="level", units_str="m-1", &
641	& long_name="Inverse of cloud inhomogeneity effective horizontal size")
642	end if
643	end if
644
645	! Write variables
646	call out_file%put("pressure_hl", thermodynamics%pressure_hl(istartcol:iendcol,:))
647
648	if (allocated(thermodynamics%h2o_sat_liq) .and. config%use_aerosols) then
649	call out_file%put("q_sat_liquid", thermodynamics%h2o_sat_liq(istartcol:iendcol,:))
650	end if
651
652	if (config%do_clouds) then
653	call out_file%put("cloud_fraction", cloud%fraction(istartcol:iendcol,:))
654	call out_file%put("overlap_param", cloud%overlap_param(istartcol:iendcol,:))
655	end if
656
657	if (config%do_sw) then
658	call out_file%put("cos_solar_zenith_angle", single_level%cos_sza(istartcol:iendcol))
659	call out_file%put("sw_albedo", sw_albedo_diffuse, do_transp=.false.)
660	call out_file%put("sw_albedo_direct", sw_albedo_direct, do_transp=.false.)
661	end if
662
663	if (config%do_lw) then
664	call out_file%put("lw_emissivity", 1.0_jprb - lw_albedo, do_transp=.false.)
665	call out_file%put("planck_hl", planck_hl)
666	call out_file%put("lw_emission", lw_emission, do_transp=.false.)
667
668	call out_file%put("od_lw", od_lw)
669	if (config%do_lw_aerosol_scattering) then
670	call out_file%put("ssa_lw", ssa_lw)
671	call out_file%put("asymmetry_lw", g_lw)
672	end if
673
674	if (config%do_clouds) then
675	call out_file%put("od_lw_cloud", od_lw_cloud)
676	if (config%do_lw_cloud_scattering) then
677	call out_file%put("ssa_lw_cloud", ssa_lw_cloud)
678	call out_file%put("asymmetry_lw_cloud", g_lw_cloud)
679	end if
680	end if
681	end if
682
683	if (config%do_sw) then
684	call out_file%put("incoming_sw", incoming_sw, do_transp=.false.)
685
686	call out_file%put("od_sw", od_sw)
687	call out_file%put("ssa_sw", ssa_sw)
688	call out_file%put("asymmetry_sw", g_sw)
689
690	if (config%do_clouds) then
691	call out_file%put("od_sw_cloud", od_sw_cloud)
692	call out_file%put("ssa_sw_cloud", ssa_sw_cloud)
693	call out_file%put("asymmetry_sw_cloud", g_sw_cloud)
694	end if
695	end if
696
697	if (config%do_clouds) then
698	if (allocated(cloud%fractional_std)) then
699	call out_file%put("fractional_std", cloud%fractional_std(istartcol:iendcol,:))
700	end if
701	if (allocated(cloud%inv_cloud_effective_size)) then
702	call out_file%put("inv_cloud_effective_size", cloud%inv_cloud_effective_size(istartcol:iendcol,:))
703	end if
704	if (allocated(cloud%inv_inhom_effective_size)) then
705	call out_file%put("inv_inhom_effective_size", cloud%inv_inhom_effective_size(istartcol:iendcol,:))
706	end if
707	end if
708
709	! Close the file
710	call out_file%close()
711
712	!$OMP END CRITICAL
713
714	end subroutine save_radiative_properties
715
716
717	!---------------------------------------------------------------------
718	! Save inputs to the radiation scheme
719	subroutine save_inputs(file_name, config, single_level, thermodynamics, &
720	& gas, cloud, aerosol, lat, lon, iverbose)
721	use yomhook, only : lhook, dr_hook
722
723	use radiation_config, only : config_type
724	use radiation_single_level, only : single_level_type
725	use radiation_thermodynamics, only : thermodynamics_type
726	use radiation_gas
727	use radiation_cloud, only : cloud_type
728	use radiation_aerosol, only : aerosol_type
729	use easy_netcdf
730
731	character(len=*), intent(in) :: file_name
732	type(config_type), intent(in) :: config
733	type(single_level_type), intent(in) :: single_level
734	type(thermodynamics_type), intent(in) :: thermodynamics
735	type(gas_type), intent(inout):: gas
736	type(cloud_type), intent(in) :: cloud
737	type(aerosol_type), optional, intent(in) :: aerosol
738	real(jprb), optional, intent(in) :: lat(:), lon(:)
739	integer, optional, intent(in) :: iverbose
740
741	real(jprb), allocatable :: mixing_ratio(:,:)
742	real(jprb), allocatable :: aerosol_mmr(:,:,:)
743	real(jprb), allocatable :: seed(:)
744	integer :: i_local_verbose
745	integer :: ncol, nlev
746	integer :: jgas
747	character(32) :: var_name, long_name
748
749	! Object for output NetCDF file
750	type(netcdf_file) :: out_file
751
752	logical :: do_aerosol
753
754	real(jprb) :: hook_handle
755
756	if (lhook) call dr_hook('radiation_save:save_inputs',0,hook_handle)
757
758	if (present(iverbose)) then
759	i_local_verbose = iverbose
760	else
761	i_local_verbose = config%iverbose
762	end if
763
764	! Work out array dimensions
765	ncol = size(thermodynamics%pressure_hl,1)
766	nlev = size(thermodynamics%pressure_hl,2)
767	nlev = nlev - 1
768
769	do_aerosol = config%use_aerosols .and. present(aerosol)
770
771	! Open the file
772	call out_file%create(trim(file_name), iverbose=i_local_verbose)
773
774	! Variables stored internally with column varying fastest, but in
775	! output file column varies most slowly so need to transpose
776	call out_file%transpose_matrices(.true.)
777
778	! In the case of aerosols we convert dimensions (ncol,nlev,ntype)
779	! in memory to (nlev,ntype,ncol) in file (in both cases the first
780	! dimension varying fastest).
781	call out_file%permute_3d_arrays( (/ 2, 3, 1 /) ) ! For aerosols
782
783	! Define dimensions
784	call out_file%define_dimension("column", ncol)
785	call out_file%define_dimension("level", nlev)
786	call out_file%define_dimension("half_level", nlev+1)
787	if (allocated(cloud%overlap_param)) then
788	call out_file%define_dimension("level_interface", nlev-1)
789	end if
790	call out_file%define_dimension("sw_albedo_band", &
791	& size(single_level%sw_albedo,2))
792	call out_file%define_dimension("lw_emissivity_band", &
793	& size(single_level%lw_emissivity,2))
794
795	if (do_aerosol) then
796	call out_file%define_dimension("aerosol_type", size(aerosol%mixing_ratio,3))
797	end if
798
799	! Put global attributes
800	call out_file%put_global_attributes( &
801	& title_str="Input profiles to the ecRad offline radiation model", &
802	& references_str="Hogan, R. J., and A. Bozzo, 2018: A flexible and efficient radiation " &
803	& //"scheme for the ECMWF model. J. Adv. Modeling Earth Sys., 10, 1990–2008", &
804	& source_str="ecRad offline radiation model")
805
806	! Define single-level variables
807	call out_file%define_variable("solar_irradiance", &
808	& units_str="W m-2", long_name="Solar irradiance at Earth's orbit")
809	if (present(lat)) then
810	call out_file%define_variable("lat", &
811	& dim1_name="column", units_str="degrees_north", long_name="Latitude")
812	end if
813	if (present(lon)) then
814	call out_file%define_variable("lon", &
815	& dim1_name="column", units_str="degrees_east", long_name="Longitude")
816	end if
817	call out_file%define_variable("skin_temperature", &
818	& dim1_name="column", units_str="K", long_name="Skin_temperature")
819	if (config%do_sw) then
820	call out_file%define_variable("cos_solar_zenith_angle", &
821	& dim1_name="column", units_str="1", &
822	& long_name="Cosine of the solar zenith angle")
823	end if
824
825	if (allocated(single_level%sw_albedo_direct)) then
826	call out_file%define_variable("sw_albedo", &
827	& dim2_name="column", dim1_name="sw_albedo_band", &
828	& units_str="1", long_name="Shortwave surface albedo to diffuse radiation")
829	call out_file%define_variable("sw_albedo_direct", &
830	& dim2_name="column", dim1_name="sw_albedo_band", &
831	& units_str="1", long_name="Shortwave surface albedo to direct radiation")
832	else
833	call out_file%define_variable("sw_albedo", &
834	& dim2_name="column", dim1_name="sw_albedo_band", &
835	& units_str="1", long_name="Shortwave surface albedo")
836
837	end if
838	call out_file%define_variable("lw_emissivity", &
839	& dim2_name="column", dim1_name="lw_emissivity_band", &
840	& units_str="1", long_name="Longwave surface emissivity")
841
842	if (allocated(single_level%iseed)) then
843	call out_file%define_variable("iseed", &
844	& dim1_name="column", units_str="1", is_double=.true., &
845	& long_name="Seed for random-number generator")
846	end if
847
848	! Define thermodynamic variables on half levels
849	call out_file%define_variable("pressure_hl", &
850	& dim2_name="column", dim1_name="half_level", &
851	& units_str="Pa", long_name="Pressure")
852	call out_file%define_variable("temperature_hl", &
853	& dim2_name="column", dim1_name="half_level", &
854	& units_str="K", long_name="Temperature")
855
856	! Define gas mixing ratios on full levels
857	call out_file%define_variable("q", &
858	& dim2_name="column", dim1_name="level", &
859	& units_str="1", long_name="Specific humidity")
860	call out_file%define_variable("o3_mmr", &
861	& dim2_name="column", dim1_name="level", &
862	& units_str="1", long_name="Ozone mass mixing ratio")
863	do jgas = 1,NMaxGases
864	if (gas%is_present(jgas) .and. jgas /= IH2O .and. jgas /= IO3) then
865	write(var_name,'(a,a)') trim(GasLowerCaseName(jgas)), '_vmr'
866	write(long_name,'(a,a)') trim(GasName(jgas)), ' volume mixing ratio'
867	call out_file%define_variable(trim(var_name), &
868	& dim2_name="column", dim1_name="level", &
869	& units_str="1", long_name=trim(long_name))
870	end if
871	end do
872
873	if (config%do_clouds) then
874	! Define cloud variables on full levels
875	call out_file%define_variable("cloud_fraction", &
876	& dim2_name="column", dim1_name="level", &
877	& units_str="1", long_name="Cloud fraction")
878	call out_file%define_variable("q_liquid", &
879	& dim2_name="column", dim1_name="level", &
880	& units_str="1", long_name="Gridbox-mean liquid water mixing ratio")
881	call out_file%define_variable("q_ice", &
882	& dim2_name="column", dim1_name="level", &
883	& units_str="1", long_name="Gridbox-mean ice water mixing ratio")
884	call out_file%define_variable("re_liquid", &
885	& dim2_name="column", dim1_name="level", &
886	& units_str="m", long_name="Ice effective radius")
887	if (allocated(cloud%re_ice)) then
888	call out_file%define_variable("re_ice", &
889	& dim2_name="column", dim1_name="level", &
890	& units_str="m", long_name="Ice effective radius")
891	end if
892	if (allocated(cloud%overlap_param)) then
893	call out_file%define_variable("overlap_param", &
894	& dim2_name="column", dim1_name="level_interface", &
895	& units_str="1", long_name="Cloud overlap parameter")
896	end if
897	if (allocated(cloud%fractional_std)) then
898	call out_file%define_variable("fractional_std", &
899	& dim2_name="column", dim1_name="level", units_str="1", &
900	& long_name="Fractional standard deviation of cloud optical depth")
901	end if
902	if (allocated(cloud%inv_cloud_effective_size)) then
903	call out_file%define_variable("inv_cloud_effective_size", &
904	& dim2_name="column", dim1_name="level", units_str="m-1", &
905	& long_name="Inverse of cloud effective horizontal size")
906	end if
907	if (allocated(cloud%inv_inhom_effective_size)) then
908	call out_file%define_variable("inv_inhom_effective_size", &
909	& dim2_name="column", dim1_name="level", units_str="m-1", &
910	& long_name="Inverse of cloud inhomogeneity effective horizontal size")
911	end if
912	end if ! do_clouds
913
914	! Define aerosol mass mixing ratio
915	if (do_aerosol) then
916	call out_file%define_variable("aerosol_mmr", &
917	& dim3_name="column", dim2_name="aerosol_type", &
918	& dim1_name="level", units_str="kg kg-1", &
919	& long_name="Aerosol mass mixing ratio")
920	end if
921
922	! Write variables
923	call out_file%put("solar_irradiance", single_level%solar_irradiance)
924	if (present(lat)) then
925	call out_file%put("lat", lat)
926	end if
927	if (present(lon)) then
928	call out_file%put("lon", lon)
929	end if
930	call out_file%put("skin_temperature", single_level%skin_temperature)
931	if (config%do_sw) then
932	call out_file%put("cos_solar_zenith_angle", single_level%cos_sza)
933	end if
934	call out_file%put("sw_albedo", single_level%sw_albedo)
935	if (allocated(single_level%sw_albedo_direct)) then
936	call out_file%put("sw_albedo_direct", single_level%sw_albedo_direct)
937	end if
938	call out_file%put("lw_emissivity", single_level%lw_emissivity)
939	if (config%do_clouds .and. allocated(single_level%iseed)) then
940	allocate(seed(ncol))
941	seed = single_level%iseed
942	call out_file%put("iseed", seed)
943	deallocate(seed)
944	end if
945
946	call out_file%put("pressure_hl", thermodynamics%pressure_hl)
947	call out_file%put("temperature_hl", thermodynamics%temperature_hl)
948
949	allocate(mixing_ratio(ncol,nlev))
950	call gas%get(IH2O, IMassMixingRatio, mixing_ratio)
951	call out_file%put("q", mixing_ratio)
952	call gas%get(IO3, IMassMixingRatio, mixing_ratio)
953	call out_file%put("o3_mmr", mixing_ratio)
954	do jgas = 1,NMaxGases
955	if (gas%is_present(jgas) .and. jgas /= IH2O .and. jgas /= IO3) then
956	write(var_name,'(a,a)') trim(GasLowerCaseName(jgas)), '_vmr'
957	call gas%get(jgas, IVolumeMixingRatio, mixing_ratio)
958	call out_file%put(trim(var_name), mixing_ratio)
959	end if
960	end do
961	deallocate(mixing_ratio)
962
963	if (config%do_clouds) then
964	call out_file%put("cloud_fraction", cloud%fraction)
965	call out_file%put("q_liquid", cloud%q_liq)
966	call out_file%put("q_ice", cloud%q_ice)
967	call out_file%put("re_liquid", cloud%re_liq)
968	if (allocated(cloud%re_ice)) then
969	call out_file%put("re_ice", cloud%re_ice)
970	end if
971	if (allocated(cloud%overlap_param)) then
972	call out_file%put("overlap_param", cloud%overlap_param)
973	end if
974	if (allocated(cloud%fractional_std)) then
975	call out_file%put("fractional_std", cloud%fractional_std)
976	end if
977	if (allocated(cloud%inv_cloud_effective_size)) then
978	call out_file%put("inv_cloud_effective_size", cloud%inv_cloud_effective_size)
979	end if
980	if (allocated(cloud%inv_inhom_effective_size)) then
981	call out_file%put("inv_inhom_effective_size", cloud%inv_inhom_effective_size)
982	end if
983	end if
984
985	if (do_aerosol) then
986	allocate(aerosol_mmr(ncol, nlev, size(aerosol%mixing_ratio,3)))
987	aerosol_mmr = 0.0_jprb
988	aerosol_mmr(:,aerosol%istartlev:aerosol%iendlev,:) = aerosol%mixing_ratio
989	call out_file%put("aerosol_mmr", aerosol_mmr)
990	deallocate(aerosol_mmr)
991	end if
992
993	! Close the file
994	call out_file%close()
995
996	if (lhook) call dr_hook('radiation_save:save_inputs',1,hook_handle)
997
998	end subroutine save_inputs
999
1000	end module radiation_save

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

Download in other formats:

Original Format