Context Navigation

radiation_save.F90 @ 5229

Last change on this file since 5229 was 4853, checked in by idelkadi, 8 months ago

Ecrad update in LMDZ / Implementation of Ecrad double call in LMDZ

version 1.6.1 (commit 210d7911380f53a788c3cad73b3cf9b4e022ef87)
interface routines between lmdz and ecrad grouped in a new "lmdz" directory
double call method redesigned so as to go through the Ecrad initialization and configuration part only once for the entire simulation
clean-up in the read.F routine: delete unitules arguments
modification of default gas model in namelist (default: ECCKD)

File size: 63.5 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	! Save final fluxes and save intermediate radiative properties
26	public :: save_fluxes, save_radiative_properties, save_inputs
27	! Save net fluxes IFS style, where upwelling fluxes are actually net down
28	public :: save_net_fluxes
29
30	contains
31
32	!---------------------------------------------------------------------
33	! Save fluxes in "flux" to NetCDF file_name, plus pressure from the
34	! thermodynamics object
35	subroutine save_fluxes(file_name, config, thermodynamics, flux, &
36	& iverbose, is_hdf5_file, experiment_name, &
37	& is_double_precision)
38
39	use yomhook, only : lhook, dr_hook, jphook
40
41	use easy_netcdf
42
43	use radiation_io, only : nulout
44	use radiation_config, only : config_type, IGasModelMonochromatic
45	use radiation_thermodynamics, only : thermodynamics_type
46	use radiation_flux, only : flux_type
47
48	character(len=*), intent(in) :: file_name
49	type(config_type), intent(in) :: config
50	type(thermodynamics_type), intent(in) :: thermodynamics
51	type(flux_type), intent(in) :: flux
52	integer, optional, intent(in) :: iverbose
53	logical, optional, intent(in) :: is_hdf5_file
54	logical, optional, intent(in) :: is_double_precision
55	character(len=*), optional, intent(in) :: experiment_name
56
57	type(netcdf_file) :: out_file
58	integer :: ncol, n_lev_plus1
59	character(5), parameter :: default_lw_units_str = 'W m-2'
60	character(5) :: lw_units_str
61	integer :: i_local_verbose
62
63	real(jphook) :: hook_handle
64
65	if (lhook) call dr_hook('radiation_save:save_fluxes',0,hook_handle)
66
67	if (present(iverbose)) then
68	i_local_verbose = iverbose
69	else
70	i_local_verbose = config%iverbose
71	end if
72
73	! Work out array dimensions
74	if (config%do_sw) then
75	ncol = size(flux%sw_up,1)
76	n_lev_plus1 = size(flux%sw_up,2)
77	elseif (config%do_lw) then
78	ncol = size(flux%lw_up,1)
79	n_lev_plus1 = size(flux%lw_up,2)
80	else
81	if (i_local_verbose >= 1) then
82	write(nulout,'(a,a,a)') 'Warning: neither longwave nor shortwave computed so ', &
83	& trim(file_name),' not written'
84	end if
85	return
86	end if
87
88	if (config%i_gas_model_lw == IGasModelMonochromatic &
89	.and. config%mono_lw_wavelength > 0.0_jprb) then
90	lw_units_str = 'W m-3'
91	else
92	lw_units_str = default_lw_units_str
93	end if
94
95	! Open the file
96	call out_file%create(trim(file_name), iverbose=i_local_verbose, is_hdf5_file=is_hdf5_file)
97
98	! Variables stored internally with column varying fastest, but in
99	! output file column varies most slowly so need to transpose
100	call out_file%transpose_matrices(.true.)
101
102	! Set default precision for file, if specified
103	if (present(is_double_precision)) then
104	call out_file%double_precision(is_double_precision)
105	end if
106
107	! Spectral fluxes in memory are dimensioned (nband,ncol,nlev), but
108	! are reoriented in the output file to be (nband,nlev,ncol), where
109	! the convention here is first dimension varying fastest
110	call out_file%permute_3d_arrays( (/ 1, 3, 2 /) )
111
112	! Define dimensions
113	call out_file%define_dimension("column", ncol)
114	call out_file%define_dimension("half_level", n_lev_plus1)
115
116	if (config%do_save_spectral_flux .or. config%do_toa_spectral_flux) then
117	if (config%do_lw) then
118	call out_file%define_dimension("band_lw", config%n_spec_lw)
119	end if
120	if (config%do_sw) then
121	call out_file%define_dimension("band_sw", config%n_spec_sw)
122	end if
123	else if (config%do_surface_sw_spectral_flux) then
124	if (config%do_sw) then
125	call out_file%define_dimension("band_sw", config%n_bands_sw)
126	end if
127	end if
128
129	if (config%do_lw .and. config%do_canopy_fluxes_lw) then
130	call out_file%define_dimension("canopy_band_lw", &
131	& size(flux%lw_dn_surf_canopy, 1))
132	end if
133	if (config%do_sw .and. config%do_canopy_fluxes_sw) then
134	call out_file%define_dimension("canopy_band_sw", &
135	& size(flux%sw_dn_diffuse_surf_canopy, 1))
136	end if
137
138	! Put global attributes
139	call out_file%put_global_attributes( &
140	& title_str="Radiative flux profiles from the ecRad offline radiation model", &
141	& references_str="Hogan, R. J., and A. Bozzo, 2018: A flexible and efficient radiation " &
142	& //"scheme for the ECMWF model. J. Adv. Modeling Earth Sys., 10, 1990–2008", &
143	& source_str="ecRad offline radiation model")
144
145	! Save "experiment" global attribute if present and not empty
146	if (present(experiment_name)) then
147	if (experiment_name /= " ") then
148	call out_file%put_global_attribute("experiment", experiment_name)
149	end if
150	end if
151
152	! Define variables
153	call out_file%define_variable("pressure_hl", &
154	& dim2_name="column", dim1_name="half_level", &
155	& units_str="Pa", long_name="Pressure", &
156	& standard_name="air_pressure")
157
158	if (config%do_lw) then
159	call out_file%define_variable("flux_up_lw", &
160	& dim2_name="column", dim1_name="half_level", &
161	& units_str=lw_units_str, long_name="Upwelling longwave flux", &
162	& standard_name="upwelling_longwave_flux_in_air")
163	call out_file%define_variable("flux_dn_lw", &
164	& dim2_name="column", dim1_name="half_level", &
165	& units_str=lw_units_str, long_name="Downwelling longwave flux", &
166	& standard_name="downwelling_longwave_flux_in_air")
167	if (config%do_clear) then
168	call out_file%define_variable("flux_up_lw_clear", &
169	& dim2_name="column", dim1_name="half_level", &
170	& units_str=lw_units_str, &
171	& long_name="Upwelling clear-sky longwave flux")
172	call out_file%define_variable("flux_dn_lw_clear", &
173	& dim2_name="column", dim1_name="half_level", &
174	& units_str=lw_units_str, &
175	& long_name="Downwelling clear-sky longwave flux")
176	end if
177
178	if (config%do_lw_derivatives) then
179	call out_file%define_variable("lw_derivative", &
180	& dim2_name="column", dim1_name="half_level", &
181	& units_str="1", &
182	& long_name="Derivative of upwelling LW flux w.r.t. surface value")
183	end if
184
185	if (config%do_save_spectral_flux) then
186	call out_file%define_variable("spectral_flux_up_lw", &
187	& dim3_name="column", dim2_name="half_level", &
188	& dim1_name="band_lw", units_str=lw_units_str, &
189	& long_name="Spectral upwelling longwave flux")
190	call out_file%define_variable("spectral_flux_dn_lw", &
191	& dim3_name="column", dim2_name="half_level", &
192	& dim1_name="band_lw", units_str=lw_units_str, &
193	& long_name="Spectral downwelling longwave flux")
194	if (config%do_clear) then
195	call out_file%define_variable("spectral_flux_up_lw_clear", &
196	& dim3_name="column", dim2_name="half_level", &
197	& dim1_name="band_lw", units_str=lw_units_str, &
198	& long_name="Spectral upwelling clear-sky longwave flux")
199	call out_file%define_variable("spectral_flux_dn_lw_clear", &
200	& dim3_name="column", dim2_name="half_level", &
201	& dim1_name="band_lw", units_str=lw_units_str, &
202	& long_name="Spectral downwelling clear-sky longwave flux")
203	end if
204	end if
205
206	if (config%do_toa_spectral_flux) then
207	call out_file%define_variable("spectral_flux_up_lw_toa", &
208	& dim2_name="column", dim1_name="band_lw", units_str="W m-2", &
209	& long_name="Spectral upwelling longwave flux at top-of-atmosphere")
210	if (config%do_clear) then
211	call out_file%define_variable("spectral_flux_up_lw_toa_clear", &
212	& dim2_name="column", dim1_name="band_lw", units_str="W m-2", &
213	& long_name="Spectral upwelling clear-sky longwave flux at top-of-atmosphere")
214	end if
215	end if
216
217	if (config%do_canopy_fluxes_lw) then
218	call out_file%define_variable("canopy_flux_dn_lw_surf", &
219	& dim2_name="column", dim1_name="canopy_band_lw", units_str=lw_units_str, &
220	& long_name="Surface downwelling longwave flux in canopy bands")
221	end if
222
223	end if
224
225	if (config%do_sw) then
226	call out_file%define_variable("flux_up_sw", &
227	& dim2_name="column", dim1_name="half_level", &
228	& units_str="W m-2", long_name="Upwelling shortwave flux", &
229	& standard_name="upwelling_shortwave_flux_in_air")
230	call out_file%define_variable("flux_dn_sw", &
231	& dim2_name="column", dim1_name="half_level", &
232	& units_str="W m-2", long_name="Downwelling shortwave flux", &
233	& standard_name="downwelling_shortwave_flux_in_air")
234	if (config%do_sw_direct) then
235	call out_file%define_variable("flux_dn_direct_sw", &
236	& dim2_name="column", dim1_name="half_level", &
237	& units_str="W m-2", &
238	& long_name="Downwelling direct shortwave flux")
239	end if
240	if (config%do_clear) then
241	call out_file%define_variable("flux_up_sw_clear", &
242	& dim2_name="column", dim1_name="half_level", &
243	& units_str="W m-2", &
244	& long_name="Upwelling clear-sky shortwave flux")
245	call out_file%define_variable("flux_dn_sw_clear", &
246	& dim2_name="column", dim1_name="half_level", &
247	& units_str="W m-2", &
248	& long_name="Downwelling clear-sky shortwave flux")
249	if (config%do_sw_direct) then
250	call out_file%define_variable("flux_dn_direct_sw_clear", &
251	& dim2_name="column", dim1_name="half_level", &
252	& units_str="W m-2", &
253	& long_name="Downwelling clear-sky direct shortwave flux")
254	end if
255	end if
256
257	if (config%do_save_spectral_flux) then
258	call out_file%define_variable("spectral_flux_up_sw", &
259	& dim3_name="column", dim2_name="half_level", &
260	& dim1_name="band_sw", units_str="W m-2", &
261	& long_name="Spectral upwelling shortwave flux")
262	call out_file%define_variable("spectral_flux_dn_sw", &
263	& dim3_name="column", dim2_name="half_level", &
264	& dim1_name="band_sw", units_str="W m-2", &
265	& long_name="Spectral downwelling shortwave flux")
266	if (config%do_sw_direct) then
267	call out_file%define_variable("spectral_flux_dn_direct_sw", &
268	& dim3_name="column", dim2_name="half_level", &
269	& dim1_name="band_sw", units_str="W m-2", &
270	& long_name="Spectral downwelling direct shortwave flux")
271	end if
272	if (config%do_clear) then
273	call out_file%define_variable("spectral_flux_up_sw_clear", &
274	& dim3_name="column", dim2_name="half_level", &
275	& dim1_name="band_sw", units_str="W m-2", &
276	& long_name="Spectral upwelling clear-sky shortwave flux")
277	call out_file%define_variable("spectral_flux_dn_sw_clear", &
278	& dim3_name="column", dim2_name="half_level", &
279	& dim1_name="band_sw", units_str="W m-2", &
280	& long_name="Spectral downwelling clear-sky shortwave flux")
281	if (config%do_sw_direct) then
282	call out_file%define_variable("spectral_flux_dn_direct_sw_clear", &
283	& dim3_name="column", dim2_name="half_level", &
284	& dim1_name="band_sw", units_str="W m-2", &
285	& long_name="Spectral downwelling clear-sky direct shortwave flux")
286	end if
287	end if
288	else if (config%do_surface_sw_spectral_flux) then
289	call out_file%define_variable("spectral_flux_dn_sw_surf", &
290	& dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
291	& long_name="Spectral downwelling shortwave flux at surface")
292	call out_file%define_variable("spectral_flux_dn_direct_sw_surf", &
293	& dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
294	& long_name="Spectral downwelling direct shortwave flux at surface")
295	if (config%do_clear) then
296	call out_file%define_variable("spectral_flux_dn_sw_surf_clear", &
297	& dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
298	& long_name="Spectral downwelling clear-sky shortwave flux at surface")
299	call out_file%define_variable("spectral_flux_dn_direct_sw_surf_clear", &
300	& dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
301	& long_name="Spectral downwelling clear-sky direct shortwave flux at surface")
302	end if
303	end if
304
305	if (config%do_toa_spectral_flux) then
306	call out_file%define_variable("spectral_flux_dn_sw_toa", &
307	& dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
308	& long_name="Spectral downwelling shortwave flux at top-of-atmosphere")
309	call out_file%define_variable("spectral_flux_up_sw_toa", &
310	& dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
311	& long_name="Spectral upwelling shortwave flux at top-of-atmosphere")
312	if (config%do_clear) then
313	call out_file%define_variable("spectral_flux_up_sw_toa_clear", &
314	& dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
315	& long_name="Spectral upwelling clear-sky shortwave flux at top-of-atmosphere")
316	end if
317	end if
318
319	if (config%do_canopy_fluxes_sw) then
320	call out_file%define_variable("canopy_flux_dn_diffuse_sw_surf", &
321	& dim2_name="column", dim1_name="canopy_band_sw", units_str="W m-2", &
322	& long_name="Surface downwelling diffuse shortwave flux in canopy bands")
323	call out_file%define_variable("canopy_flux_dn_direct_sw_surf", &
324	& dim2_name="column", dim1_name="canopy_band_sw", units_str="W m-2", &
325	& long_name="Surface downwelling direct shortwave flux in canopy bands")
326	end if
327
328	end if
329
330	if (config%do_lw .and. config%do_clouds) then
331	call out_file%define_variable("cloud_cover_lw", &
332	& dim1_name="column", units_str="1", &
333	& long_name="Total cloud cover diagnosed by longwave solver", &
334	& standard_name="cloud_area_fraction")
335	end if
336	if (config%do_sw .and. config%do_clouds) then
337	call out_file%define_variable("cloud_cover_sw", &
338	& dim1_name="column", units_str="1", &
339	& long_name="Total cloud cover diagnosed by shortwave solver", &
340	& standard_name="cloud_area_fraction")
341	end if
342
343	! Write variables
344
345	call out_file%put("pressure_hl", thermodynamics%pressure_hl)
346
347	if (config%do_lw) then
348	call out_file%put("flux_up_lw", flux%lw_up)
349	call out_file%put("flux_dn_lw", flux%lw_dn)
350	if (config%do_clear) then
351	call out_file%put("flux_up_lw_clear", flux%lw_up_clear)
352	call out_file%put("flux_dn_lw_clear", flux%lw_dn_clear)
353	end if
354
355	if (config%do_lw_derivatives) then
356	call out_file%put("lw_derivative", flux%lw_derivatives)
357	end if
358
359	if (config%do_save_spectral_flux) then
360	call out_file%put("spectral_flux_up_lw", flux%lw_up_band)
361	call out_file%put("spectral_flux_dn_lw", flux%lw_dn_band)
362	if (config%do_clear) then
363	call out_file%put("spectral_flux_up_lw_clear", flux%lw_up_clear_band)
364	call out_file%put("spectral_flux_dn_lw_clear", flux%lw_dn_clear_band)
365	end if
366	end if
367
368	if (config%do_toa_spectral_flux) then
369	call out_file%put("spectral_flux_up_lw_toa", flux%lw_up_toa_band, &
370	& do_transp=.false.)
371	if (config%do_clear) then
372	call out_file%put("spectral_flux_up_lw_toa_clear", flux%lw_up_toa_clear_band, &
373	& do_transp=.false.)
374	end if
375	end if
376
377	if (config%do_canopy_fluxes_lw) then
378	call out_file%put("canopy_flux_dn_lw_surf", flux%lw_dn_surf_canopy, &
379	& do_transp = .false.)
380	end if
381
382	end if
383
384	if (config%do_sw) then
385	call out_file%put("flux_up_sw", flux%sw_up)
386	call out_file%put("flux_dn_sw", flux%sw_dn)
387	if (config%do_sw_direct) then
388	call out_file%put("flux_dn_direct_sw", flux%sw_dn_direct)
389	end if
390	if (config%do_clear) then
391	call out_file%put("flux_up_sw_clear", flux%sw_up_clear)
392	call out_file%put("flux_dn_sw_clear", flux%sw_dn_clear)
393	if (config%do_sw_direct) then
394	call out_file%put("flux_dn_direct_sw_clear", flux%sw_dn_direct_clear)
395	end if
396	end if
397
398	if (config%do_save_spectral_flux) then
399	call out_file%put("spectral_flux_up_sw", flux%sw_up_band)
400	call out_file%put("spectral_flux_dn_sw", flux%sw_dn_band)
401	if (config%do_sw_direct) then
402	call out_file%put("spectral_flux_dn_direct_sw", &
403	& flux%sw_dn_direct_band)
404	end if
405	if (config%do_clear) then
406	call out_file%put("spectral_flux_up_sw_clear", flux%sw_up_clear_band)
407	call out_file%put("spectral_flux_dn_sw_clear", flux%sw_dn_clear_band)
408	if (config%do_sw_direct) then
409	call out_file%put("spectral_flux_dn_direct_sw_clear", &
410	& flux%sw_dn_direct_clear_band)
411	end if
412	end if
413	else if (config%do_surface_sw_spectral_flux) then
414	call out_file%put("spectral_flux_dn_sw_surf", flux%sw_dn_surf_band, &
415	& do_transp=.false.)
416	call out_file%put("spectral_flux_dn_direct_sw_surf", flux%sw_dn_direct_surf_band, &
417	& do_transp=.false.)
418	if (config%do_clear) then
419	call out_file%put("spectral_flux_dn_sw_surf_clear", flux%sw_dn_surf_clear_band, &
420	& do_transp=.false.)
421	call out_file%put("spectral_flux_dn_direct_sw_surf_clear", &
422	& flux%sw_dn_direct_surf_clear_band, do_transp=.false.)
423	end if
424	end if
425
426	if (config%do_toa_spectral_flux) then
427	call out_file%put("spectral_flux_dn_sw_toa", flux%sw_dn_toa_band, &
428	& do_transp=.false.)
429	call out_file%put("spectral_flux_up_sw_toa", flux%sw_up_toa_band, &
430	& do_transp=.false.)
431	if (config%do_clear) then
432	call out_file%put("spectral_flux_up_sw_toa_clear", flux%sw_up_toa_clear_band, &
433	& do_transp=.false.)
434	end if
435	end if
436
437	if (config%do_canopy_fluxes_sw) then
438	call out_file%put("canopy_flux_dn_diffuse_sw_surf", flux%sw_dn_diffuse_surf_canopy, &
439	& do_transp = .false.)
440	call out_file%put("canopy_flux_dn_direct_sw_surf", flux%sw_dn_direct_surf_canopy, &
441	& do_transp = .false.)
442	end if
443
444	end if
445
446	if (config%do_lw .and. config%do_clouds) then
447	call out_file%put("cloud_cover_lw", flux%cloud_cover_lw)
448	end if
449	if (config%do_sw .and. config%do_clouds) then
450	call out_file%put("cloud_cover_sw", flux%cloud_cover_sw)
451	end if
452
453	! Close file
454	call out_file%close()
455
456	if (lhook) call dr_hook('radiation_save:save_fluxes',1,hook_handle)
457
458	end subroutine save_fluxes
459
460
461	!---------------------------------------------------------------------
462	! Save IFS-style net fluxes in "flux" to NetCDF file_name, plus
463	! pressure from the thermodynamics object
464	subroutine save_net_fluxes(file_name, config, thermodynamics, flux, &
465	& iverbose, is_hdf5_file, experiment_name, &
466	& is_double_precision)
467
468	use yomhook, only : lhook, dr_hook, jphook
469
470	use easy_netcdf
471
472	use radiation_io, only : nulout
473	use radiation_config, only : config_type, IGasModelMonochromatic
474	use radiation_thermodynamics, only : thermodynamics_type
475	use radiation_flux, only : flux_type
476
477	character(len=*), intent(in) :: file_name
478	type(config_type), intent(in) :: config
479	type(thermodynamics_type), intent(in) :: thermodynamics
480	type(flux_type), intent(in) :: flux
481	integer, optional, intent(in) :: iverbose
482	logical, optional, intent(in) :: is_hdf5_file
483	logical, optional, intent(in) :: is_double_precision
484	character(len=*), optional, intent(in) :: experiment_name
485
486	type(netcdf_file) :: out_file
487	integer :: ncol, n_lev_plus1
488	character(5), parameter :: default_lw_units_str = 'W m-2'
489	character(5) :: lw_units_str
490	integer :: i_local_verbose
491
492	real(jphook) :: hook_handle
493
494	if (lhook) call dr_hook('radiation_save:save_net_fluxes',0,hook_handle)
495
496	if (present(iverbose)) then
497	i_local_verbose = iverbose
498	else
499	i_local_verbose = config%iverbose
500	end if
501
502	! Work out array dimensions
503	if (config%do_sw) then
504	ncol = size(flux%sw_up,1)
505	n_lev_plus1 = size(flux%sw_up,2)
506	elseif (config%do_lw) then
507	ncol = size(flux%lw_up,1)
508	n_lev_plus1 = size(flux%lw_up,2)
509	else
510	if (i_local_verbose >= 1) then
511	write(nulout,'(a,a,a)') 'Warning: neither longwave nor shortwave computed so ', &
512	& file_name,' not written'
513	end if
514	return
515	end if
516
517	if (config%i_gas_model_lw == IGasModelMonochromatic &
518	.and. config%mono_lw_wavelength > 0.0_jprb) then
519	lw_units_str = 'W m-3'
520	else
521	lw_units_str = default_lw_units_str
522	end if
523
524	! Open the file
525	call out_file%create(trim(file_name), iverbose=i_local_verbose, is_hdf5_file=is_hdf5_file)
526
527	! Variables stored internally with column varying fastest, but in
528	! output file column varies most slowly so need to transpose
529	call out_file%transpose_matrices(.true.)
530
531	! Set default precision for file, if specified
532	if (present(is_double_precision)) then
533	call out_file%double_precision(is_double_precision)
534	end if
535
536	! Spectral fluxes in memory are dimensioned (nband,ncol,nlev), but
537	! are reoriented in the output file to be (nband,nlev,ncol), where
538	! the convention here is first dimension varying fastest
539	call out_file%permute_3d_arrays( (/ 1, 3, 2 /) )
540
541	! Define dimensions
542	call out_file%define_dimension("column", ncol)
543	call out_file%define_dimension("half_level", n_lev_plus1)
544
545	if (config%do_lw .and. config%do_canopy_fluxes_lw) then
546	call out_file%define_dimension("canopy_band_lw", &
547	& size(flux%lw_dn_surf_canopy, 1))
548	end if
549	if (config%do_sw .and. config%do_canopy_fluxes_sw) then
550	call out_file%define_dimension("canopy_band_sw", &
551	& size(flux%sw_dn_diffuse_surf_canopy, 1))
552	end if
553
554	! Put global attributes
555	call out_file%put_global_attributes( &
556	& title_str="Radiative flux profiles from the ecRad offline radiation model", &
557	& references_str="Hogan, R. J., and A. Bozzo, 2018: A flexible and efficient radiation " &
558	& //"scheme for the ECMWF model. J. Adv. Modeling Earth Sys., 10, 1990–2008", &
559	& source_str="ecRad offline radiation model")
560
561	! Save "experiment" global attribute if present and not empty
562	if (present(experiment_name)) then
563	if (experiment_name /= " ") then
564	call out_file%put_global_attribute("experiment", experiment_name)
565	end if
566	end if
567
568	! Define variables
569	call out_file%define_variable("pressure_hl", &
570	& dim2_name="column", dim1_name="half_level", &
571	& units_str="Pa", long_name="Pressure", &
572	& standard_name="air_pressure")
573
574	if (config%do_lw) then
575	call out_file%define_variable("flux_net_lw", &
576	& dim2_name="column", dim1_name="half_level", &
577	& units_str=lw_units_str, long_name="Net downward longwave flux", &
578	& standard_name="net_downward_longwave_flux_in_air")
579	call out_file%define_variable("flux_dn_lw_surf", &
580	& dim1_name="column", &
581	& units_str=lw_units_str, long_name="Surface downwelling longwave flux", &
582	& standard_name="surface_downwelling_longwave_flux_in_air")
583	if (config%do_clear) then
584	call out_file%define_variable("flux_net_lw_clear", &
585	& dim2_name="column", dim1_name="half_level", &
586	& units_str=lw_units_str, &
587	& long_name="Net downward clear-sky longwave flux", &
588	& standard_name="net_downward_longwave_flux_in_air_assuming_clear_sky")
589	call out_file%define_variable("flux_dn_lw_clear_surf", &
590	& dim1_name="column", &
591	& units_str=lw_units_str, &
592	& long_name="Surface downwelling clear-sky longwave flux", &
593	& standard_name="surface_downwelling_longwave_flux_in_air_assuming_clear_sky")
594	end if
595
596	if (config%do_lw_derivatives) then
597	call out_file%define_variable("lw_derivative", &
598	& dim2_name="column", dim1_name="half_level", &
599	& units_str="1", &
600	& long_name="Derivative of upwelling LW flux w.r.t. surface value")
601	end if
602
603	if (config%do_canopy_fluxes_lw) then
604	call out_file%define_variable("canopy_flux_dn_lw_surf", &
605	& dim2_name="column", dim1_name="canopy_band_lw", units_str=lw_units_str, &
606	& long_name="Surface downwelling longwave flux in canopy bands")
607	end if
608
609	end if
610
611	if (config%do_sw) then
612	call out_file%define_variable("flux_net_sw", &
613	& dim2_name="column", dim1_name="half_level", &
614	& units_str="W m-2", long_name="Net downward shortwave flux", &
615	& standard_name="net_downward_shortwave_flux_in_air")
616	call out_file%define_variable("flux_dn_sw_surf", &
617	& dim1_name="column", &
618	& units_str="W m-2", long_name="Surface downwelling shortwave flux", &
619	& standard_name="surface_downwelling_shortwave_flux_in_air")
620	call out_file%define_variable("flux_dn_sw_toa", &
621	& dim1_name="column", &
622	& units_str="W m-2", long_name="Top-of-atmosphere downwelling shortwave flux", &
623	& standard_name="toa_incoming_shortwave_flux")
624	if (config%do_sw_direct) then
625	call out_file%define_variable("flux_dn_direct_sw_surf", &
626	& dim1_name="column", &
627	& units_str="W m-2", &
628	& long_name="Surface downwelling direct shortwave flux")
629	end if
630	if (config%do_clear) then
631	call out_file%define_variable("flux_net_sw_clear", &
632	& dim2_name="column", dim1_name="half_level", &
633	& units_str="W m-2", &
634	& long_name="Net downward clear-sky shortwave flux")
635	call out_file%define_variable("flux_dn_sw_clear_surf", &
636	& dim1_name="column", &
637	& units_str="W m-2", &
638	& long_name="Surface downwelling clear-sky shortwave flux")
639	if (config%do_sw_direct) then
640	call out_file%define_variable("flux_dn_direct_sw_clear_surf", &
641	& dim1_name="column", &
642	& units_str="W m-2", &
643	& long_name="Surface downwelling clear-sky direct shortwave flux")
644	end if
645	end if
646
647	if (config%do_canopy_fluxes_sw) then
648	call out_file%define_variable("canopy_flux_dn_diffuse_sw_surf", &
649	& dim2_name="column", dim1_name="canopy_band_sw", units_str="W m-2", &
650	& long_name="Surface downwelling diffuse shortwave flux in canopy bands")
651	call out_file%define_variable("canopy_flux_dn_direct_sw_surf", &
652	& dim2_name="column", dim1_name="canopy_band_sw", units_str="W m-2", &
653	& long_name="Surface downwelling direct shortwave flux in canopy bands")
654	end if
655
656	end if
657
658	! Write variables
659
660	call out_file%put("pressure_hl", thermodynamics%pressure_hl)
661
662	if (config%do_lw) then
663	call out_file%put("flux_net_lw", flux%lw_dn-flux%lw_up)
664	call out_file%put("flux_dn_lw_surf", flux%lw_dn(:,n_lev_plus1))
665	if (config%do_clear) then
666	call out_file%put("flux_net_lw_clear", flux%lw_dn_clear-flux%lw_up_clear)
667	call out_file%put("flux_dn_lw_clear_surf", flux%lw_dn_clear(:,n_lev_plus1))
668	end if
669
670	if (config%do_lw_derivatives) then
671	call out_file%put("lw_derivative", flux%lw_derivatives)
672	end if
673
674	if (config%do_canopy_fluxes_lw) then
675	call out_file%put("canopy_flux_dn_lw_surf", flux%lw_dn_surf_canopy, &
676	& do_transp = .false.)
677	end if
678
679	end if
680
681	if (config%do_sw) then
682	call out_file%put("flux_net_sw", flux%sw_dn-flux%sw_up)
683	call out_file%put("flux_dn_sw_surf", flux%sw_dn(:,n_lev_plus1))
684	call out_file%put("flux_dn_sw_toa", flux%sw_dn(:,1))
685	if (config%do_sw_direct) then
686	call out_file%put("flux_dn_direct_sw_surf", flux%sw_dn_direct(:,n_lev_plus1))
687	end if
688	if (config%do_clear) then
689	call out_file%put("flux_net_sw_clear", flux%sw_dn_clear-flux%sw_up_clear)
690	call out_file%put("flux_dn_sw_clear_surf", flux%sw_dn_clear(:,n_lev_plus1))
691	if (config%do_sw_direct) then
692	call out_file%put("flux_dn_direct_sw_clear_surf", flux%sw_dn_direct_clear(:,n_lev_plus1))
693	end if
694	end if
695
696	if (config%do_canopy_fluxes_sw) then
697	call out_file%put("canopy_flux_dn_diffuse_sw_surf", flux%sw_dn_diffuse_surf_canopy, &
698	& do_transp = .false.)
699	call out_file%put("canopy_flux_dn_direct_sw_surf", flux%sw_dn_direct_surf_canopy, &
700	& do_transp = .false.)
701	end if
702
703	end if
704
705	! Close file
706	call out_file%close()
707
708	if (lhook) call dr_hook('radiation_save:save_net_fluxes',1,hook_handle)
709
710	end subroutine save_net_fluxes
711
712
713	!---------------------------------------------------------------------
714	! Save intermediate radiative properties, specifically the
715	! scattering and absorption properties at each g-point/band
716	subroutine save_radiative_properties(file_name, nlev, &
717	& istartcol, iendcol, &
718	& config, single_level, thermodynamics, cloud, &
719	& planck_hl, lw_emission, lw_albedo, &
720	& sw_albedo_direct, sw_albedo_diffuse, &
721	& incoming_sw, &
722	& od_lw, ssa_lw, g_lw, &
723	& od_sw, ssa_sw, g_sw, &
724	& od_lw_cloud, ssa_lw_cloud, g_lw_cloud, &
725	& od_sw_cloud, ssa_sw_cloud, g_sw_cloud)
726
727	use radiation_config, only : config_type
728	use radiation_single_level, only : single_level_type
729	use radiation_thermodynamics,only : thermodynamics_type
730	use radiation_cloud, only : cloud_type
731	use easy_netcdf
732
733	character(len=*), intent(in) :: file_name
734	type(config_type), intent(in) :: config
735	type(single_level_type), intent(in) :: single_level
736	type(thermodynamics_type),intent(in) :: thermodynamics
737	type(cloud_type), intent(in) :: cloud
738
739	integer, intent(in) :: nlev, istartcol, iendcol
740
741	! Input variables, as defined in radiation_interface.F90
742
743	! Layer optical depth, single scattering albedo and asymmetry factor of
744	! gases and aerosols at each shortwave g-point
745	real(jprb), intent(in), dimension(config%n_g_sw,nlev,istartcol:iendcol) :: od_sw, ssa_sw, g_sw
746
747	! Layer optical depth, single scattering albedo and asymmetry factor of
748	! hydrometeors in each shortwave band
749	real(jprb), intent(in), dimension(config%n_bands_sw,nlev,istartcol:iendcol) :: &
750	& od_sw_cloud, ssa_sw_cloud, g_sw_cloud
751
752	! Direct and diffuse surface albedo, and the incoming shortwave
753	! flux into a plane perpendicular to the incoming radiation at
754	! top-of-atmosphere in each of the shortwave g-points
755	real(jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol) &
756	& :: sw_albedo_direct, sw_albedo_diffuse, incoming_sw
757
758	! Layer optical depth, single scattering albedo and asymmetry factor of
759	! gases and aerosols at each longwave g-point, where the latter
760	! two variables are only defined if aerosol longwave scattering is
761	! enabled (otherwise both are treated as zero).
762	real(jprb), intent(in), dimension(config%n_g_lw,nlev,istartcol:iendcol) :: od_lw
763	real(jprb), intent(in), dimension(config%n_g_lw_if_scattering,nlev,istartcol:iendcol) :: &
764	& ssa_lw, g_lw
765
766	! Layer optical depth, single scattering albedo and asymmetry factor of
767	! hydrometeors in each longwave band, where the latter two
768	! variables are only defined if hydrometeor longwave scattering is
769	! enabled (otherwise both are treated as zero).
770	real(jprb), intent(in), dimension(config%n_bands_lw,nlev,istartcol:iendcol) :: od_lw_cloud
771	real(jprb), intent(in), dimension(config%n_bands_lw_if_scattering,nlev,istartcol:iendcol) :: &
772	& ssa_lw_cloud, g_lw_cloud
773
774	! The Planck function (emitted flux from a black body) at half
775	! levels and at the surface at each longwave g-point
776	real(jprb), intent(in), dimension(config%n_g_lw,nlev+1,istartcol:iendcol) :: planck_hl
777
778	! Emission (Planck*emissivity) and albedo (1-emissivity) at the
779	! surface at each longwave g-point
780	real(jprb), intent(in), dimension(config%n_g_lw, istartcol:iendcol) :: lw_emission, lw_albedo
781
782	! Local variables
783
784	integer :: n_col_local ! Number of columns from istartcol to iendcol
785
786	! Object for output NetCDF file
787	type(netcdf_file) :: out_file
788
789	n_col_local = iendcol + 1 - istartcol
790
791	! Alas the NetCDF library is not thread-safe for writing, so we
792	! must write radiative-property files serially
793
794	!$OMP CRITICAL
795
796	! Open the file
797	call out_file%create(trim(file_name), iverbose=config%iverbose)
798
799	! Configure matrix and 3D-array orientation
800	call out_file%transpose_matrices(.true.)
801
802	! Sometimes the Planck function values are very large or small
803	! even if the fluxes are within a manageable range
804	call out_file%double_precision(.true.)
805
806	! Define dimensions
807	! call out_file%define_dimension("column", n_col_local)
808	call out_file%define_dimension("column", 0) ! "Unlimited" dimension
809	call out_file%define_dimension("level", nlev)
810	call out_file%define_dimension("half_level", nlev+1)
811	if (config%do_clouds) then
812	call out_file%define_dimension("level_interface", nlev-1)
813	end if
814
815	if (config%do_lw) then
816	call out_file%define_dimension("gpoint_lw", config%n_g_lw)
817	if (config%do_clouds) then
818	call out_file%define_dimension("band_lw", config%n_bands_lw)
819	end if
820	end if
821	if (config%do_sw) then
822	call out_file%define_dimension("gpoint_sw", config%n_g_sw)
823	if (config%do_clouds) then
824	call out_file%define_dimension("band_sw", config%n_bands_sw)
825	end if
826	end if
827
828	! Put global attributes
829	call out_file%put_global_attributes( &
830	& title_str="Spectral radiative properties from the ecRad offline radiation model", &
831	& references_str="Hogan, R. J., and A. Bozzo, 2018: A flexible and efficient radiation " &
832	& //"scheme for the ECMWF model. J. Adv. Modeling Earth Sys., 10, 1990–2008", &
833	& source_str="ecRad offline radiation model")
834
835	! Define variables
836	call out_file%define_variable("pressure_hl", &
837	& dim2_name="column", dim1_name="half_level", &
838	& units_str="Pa", long_name="Pressure on half-levels")
839
840	if (allocated(thermodynamics%h2o_sat_liq) .and. config%use_aerosols) then
841	call out_file%define_variable("q_sat_liquid", &
842	& dim2_name="column", dim1_name="level", &
843	& units_str="kg kg-1", long_name="Specific humidity at liquid saturation")
844	end if
845
846	if (config%do_sw) then
847	call out_file%define_variable("cos_solar_zenith_angle", &
848	& dim1_name="column", units_str="1", &
849	& long_name="Cosine of the solar zenith angle")
850	end if
851
852	if (config%do_clouds) then
853	call out_file%define_variable("cloud_fraction", &
854	& dim2_name="column", dim1_name="level", &
855	& units_str="1", long_name="Cloud fraction")
856	call out_file%define_variable("overlap_param", &
857	& dim2_name="column", dim1_name="level_interface", &
858	& units_str="1", long_name="Cloud overlap parameter")
859	end if
860
861	if (config%do_lw) then
862	call out_file%define_variable("planck_hl", &
863	& dim3_name="column", dim2_name="half_level", dim1_name="gpoint_lw", &
864	& units_str="W m-2", long_name="Planck function on half-levels")
865	call out_file%define_variable("lw_emission", &
866	& dim2_name="column", dim1_name="gpoint_lw", &
867	& units_str="W m-2", long_name="Longwave surface emission")
868	call out_file%define_variable("lw_emissivity", &
869	& dim2_name="column", dim1_name="gpoint_lw", &
870	& units_str="1", long_name="Surface longwave emissivity")
871
872	call out_file%define_variable("od_lw", &
873	& dim3_name="column", dim2_name="level", dim1_name="gpoint_lw", &
874	& units_str="1", long_name="Clear-sky longwave optical depth")
875	if (config%do_lw_aerosol_scattering) then
876	call out_file%define_variable("ssa_lw", &
877	& dim3_name="column", dim2_name="level", dim1_name="gpoint_lw", &
878	& units_str="1", long_name="Clear-sky longwave single scattering albedo")
879	call out_file%define_variable("asymmetry_lw", &
880	& dim3_name="column", dim2_name="level", dim1_name="gpoint_lw", &
881	& units_str="1", long_name="Clear-sky longwave asymmetry factor")
882	end if
883
884	if (config%do_clouds) then
885	call out_file%define_variable("od_lw_cloud", &
886	& dim3_name="column", dim2_name="level", dim1_name="band_lw", &
887	& units_str="1", long_name="In-cloud longwave optical depth")
888	if (config%do_lw_cloud_scattering) then
889	call out_file%define_variable("ssa_lw_cloud", &
890	& dim3_name="column", dim2_name="level", dim1_name="band_lw", &
891	& units_str="1", long_name="Cloud longwave single scattering albedo")
892	call out_file%define_variable("asymmetry_lw_cloud", &
893	& dim3_name="column", dim2_name="level", dim1_name="band_lw", &
894	& units_str="1", long_name="Cloud longwave asymmetry factor")
895	end if
896	end if ! do_clouds
897	end if ! do_lw
898
899	if (config%do_sw) then
900	call out_file%define_variable("incoming_sw", &
901	& dim2_name="column", dim1_name="gpoint_sw", &
902	& units_str="W m-2", long_name="Incoming shortwave flux at top-of-atmosphere in direction of sun")
903
904	call out_file%define_variable("sw_albedo", &
905	& dim2_name="column", dim1_name="gpoint_sw", &
906	& units_str="1", long_name="Surface shortwave albedo to diffuse radiation")
907	call out_file%define_variable("sw_albedo_direct", &
908	& dim2_name="column", dim1_name="gpoint_sw", &
909	& units_str="1", long_name="Surface shortwave albedo to direct radiation")
910
911	call out_file%define_variable("od_sw", &
912	& dim3_name="column", dim2_name="level", dim1_name="gpoint_sw", &
913	& units_str="1", long_name="Clear-sky shortwave optical depth")
914	call out_file%define_variable("ssa_sw", &
915	& dim3_name="column", dim2_name="level", dim1_name="gpoint_sw", &
916	& units_str="1", long_name="Clear-sky shortwave single scattering albedo")
917	call out_file%define_variable("asymmetry_sw", &
918	& dim3_name="column", dim2_name="level", dim1_name="gpoint_sw", &
919	& units_str="1", long_name="Clear-sky shortwave asymmetry factor")
920
921	if (config%do_clouds) then
922	call out_file%define_variable("od_sw_cloud", &
923	& dim3_name="column", dim2_name="level", dim1_name="band_sw", &
924	& units_str="1", long_name="In-cloud shortwave optical depth")
925	call out_file%define_variable("ssa_sw_cloud", &
926	& dim3_name="column", dim2_name="level", dim1_name="band_sw", &
927	& units_str="1", long_name="Cloud shortwave single scattering albedo")
928	call out_file%define_variable("asymmetry_sw_cloud", &
929	& dim3_name="column", dim2_name="level", dim1_name="band_sw", &
930	& units_str="1", long_name="Cloud shortwave asymmetry factor")
931	end if
932	end if
933
934	if (config%do_clouds) then
935	if (allocated(cloud%fractional_std)) then
936	call out_file%define_variable("fractional_std", &
937	& dim2_name="column", dim1_name="level", units_str="1", &
938	& long_name="Fractional standard deviation of cloud optical depth")
939	end if
940	if (allocated(cloud%inv_cloud_effective_size)) then
941	call out_file%define_variable("inv_cloud_effective_size", &
942	& dim2_name="column", dim1_name="level", units_str="m-1", &
943	& long_name="Inverse of cloud effective horizontal size")
944	end if
945	if (allocated(cloud%inv_inhom_effective_size)) then
946	call out_file%define_variable("inv_inhom_effective_size", &
947	& dim2_name="column", dim1_name="level", units_str="m-1", &
948	& long_name="Inverse of cloud inhomogeneity effective horizontal size")
949	end if
950	end if
951
952	! Write variables
953	call out_file%put("pressure_hl", thermodynamics%pressure_hl(istartcol:iendcol,:))
954
955	if (allocated(thermodynamics%h2o_sat_liq) .and. config%use_aerosols) then
956	call out_file%put("q_sat_liquid", thermodynamics%h2o_sat_liq(istartcol:iendcol,:))
957	end if
958
959	if (config%do_clouds) then
960	call out_file%put("cloud_fraction", cloud%fraction(istartcol:iendcol,:))
961	call out_file%put("overlap_param", cloud%overlap_param(istartcol:iendcol,:))
962	end if
963
964	if (config%do_sw) then
965	call out_file%put("cos_solar_zenith_angle", single_level%cos_sza(istartcol:iendcol))
966	call out_file%put("sw_albedo", sw_albedo_diffuse, do_transp=.false.)
967	call out_file%put("sw_albedo_direct", sw_albedo_direct, do_transp=.false.)
968	end if
969
970	if (config%do_lw) then
971	call out_file%put("lw_emissivity", 1.0_jprb - lw_albedo, do_transp=.false.)
972	call out_file%put("planck_hl", planck_hl)
973	call out_file%put("lw_emission", lw_emission, do_transp=.false.)
974
975	call out_file%put("od_lw", od_lw)
976	if (config%do_lw_aerosol_scattering) then
977	call out_file%put("ssa_lw", ssa_lw)
978	call out_file%put("asymmetry_lw", g_lw)
979	end if
980
981	if (config%do_clouds) then
982	call out_file%put("od_lw_cloud", od_lw_cloud)
983	if (config%do_lw_cloud_scattering) then
984	call out_file%put("ssa_lw_cloud", ssa_lw_cloud)
985	call out_file%put("asymmetry_lw_cloud", g_lw_cloud)
986	end if
987	end if
988	end if
989
990	if (config%do_sw) then
991	call out_file%put("incoming_sw", incoming_sw, do_transp=.false.)
992
993	call out_file%put("od_sw", od_sw)
994	call out_file%put("ssa_sw", ssa_sw)
995	call out_file%put("asymmetry_sw", g_sw)
996
997	if (config%do_clouds) then
998	call out_file%put("od_sw_cloud", od_sw_cloud)
999	call out_file%put("ssa_sw_cloud", ssa_sw_cloud)
1000	call out_file%put("asymmetry_sw_cloud", g_sw_cloud)
1001	end if
1002	end if
1003
1004	if (config%do_clouds) then
1005	if (allocated(cloud%fractional_std)) then
1006	call out_file%put("fractional_std", cloud%fractional_std(istartcol:iendcol,:))
1007	end if
1008	if (allocated(cloud%inv_cloud_effective_size)) then
1009	call out_file%put("inv_cloud_effective_size", cloud%inv_cloud_effective_size(istartcol:iendcol,:))
1010	end if
1011	if (allocated(cloud%inv_inhom_effective_size)) then
1012	call out_file%put("inv_inhom_effective_size", cloud%inv_inhom_effective_size(istartcol:iendcol,:))
1013	end if
1014	end if
1015
1016	! Close the file
1017	call out_file%close()
1018
1019	!$OMP END CRITICAL
1020
1021	end subroutine save_radiative_properties
1022
1023
1024	!---------------------------------------------------------------------
1025	! Save inputs to the radiation scheme
1026	subroutine save_inputs(file_name, config, single_level, thermodynamics, &
1027	& gas, cloud, aerosol, lat, lon, iverbose)
1028	use yomhook, only : lhook, dr_hook, jphook
1029
1030	use radiation_config, only : config_type
1031	use radiation_single_level, only : single_level_type
1032	use radiation_thermodynamics, only : thermodynamics_type
1033	use radiation_gas
1034	use radiation_cloud, only : cloud_type
1035	use radiation_aerosol, only : aerosol_type
1036	use easy_netcdf
1037
1038	character(len=*), intent(in) :: file_name
1039	type(config_type), intent(in) :: config
1040	type(single_level_type), intent(in) :: single_level
1041	type(thermodynamics_type), intent(in) :: thermodynamics
1042	type(gas_type), intent(inout):: gas
1043	type(cloud_type), intent(in) :: cloud
1044	type(aerosol_type), optional, intent(in) :: aerosol
1045	real(jprb), optional, intent(in) :: lat(:), lon(:)
1046	integer, optional, intent(in) :: iverbose
1047
1048	real(jprb), allocatable :: mixing_ratio(:,:)
1049	real(jprb), allocatable :: aerosol_mmr(:,:,:)
1050	real(jprb), allocatable :: seed(:)
1051	integer :: i_local_verbose
1052	integer :: ncol, nlev
1053	integer :: jgas
1054	character(32) :: var_name, long_name
1055
1056	! Object for output NetCDF file
1057	type(netcdf_file) :: out_file
1058
1059	logical :: do_aerosol
1060
1061	real(jphook) :: hook_handle
1062
1063	if (lhook) call dr_hook('radiation_save:save_inputs',0,hook_handle)
1064
1065	if (present(iverbose)) then
1066	i_local_verbose = iverbose
1067	else
1068	i_local_verbose = config%iverbose
1069	end if
1070
1071	! Work out array dimensions
1072	ncol = size(thermodynamics%pressure_hl,1)
1073	nlev = size(thermodynamics%pressure_hl,2)
1074	nlev = nlev - 1
1075
1076	do_aerosol = config%use_aerosols .and. present(aerosol)
1077
1078	! Open the file
1079	call out_file%create(trim(file_name), iverbose=i_local_verbose)
1080
1081	! Variables stored internally with column varying fastest, but in
1082	! output file column varies most slowly so need to transpose
1083	call out_file%transpose_matrices(.true.)
1084
1085	! In the case of aerosols we convert dimensions (ncol,nlev,ntype)
1086	! in memory to (nlev,ntype,ncol) in file (in both cases the first
1087	! dimension varying fastest).
1088	call out_file%permute_3d_arrays( (/ 2, 3, 1 /) ) ! For aerosols
1089
1090	! Define dimensions
1091	call out_file%define_dimension("column", ncol)
1092	call out_file%define_dimension("level", nlev)
1093	call out_file%define_dimension("half_level", nlev+1)
1094	if (allocated(cloud%overlap_param)) then
1095	call out_file%define_dimension("level_interface", nlev-1)
1096	end if
1097	call out_file%define_dimension("sw_albedo_band", &
1098	& size(single_level%sw_albedo,2))
1099	call out_file%define_dimension("lw_emissivity_band", &
1100	& size(single_level%lw_emissivity,2))
1101
1102	if (do_aerosol) then
1103	call out_file%define_dimension("aerosol_type", size(aerosol%mixing_ratio,3))
1104	end if
1105
1106	! Put global attributes
1107	call out_file%put_global_attributes( &
1108	& title_str="Input profiles to the ecRad offline radiation model", &
1109	& references_str="Hogan, R. J., and A. Bozzo, 2018: A flexible and efficient radiation " &
1110	& //"scheme for the ECMWF model. J. Adv. Modeling Earth Sys., 10, 1990–2008", &
1111	& source_str="ecRad offline radiation model")
1112
1113	! Define single-level variables
1114	call out_file%define_variable("solar_irradiance", &
1115	& units_str="W m-2", long_name="Solar irradiance at Earth's orbit")
1116	if (present(lat)) then
1117	call out_file%define_variable("lat", &
1118	& dim1_name="column", units_str="degrees_north", long_name="Latitude")
1119	end if
1120	if (present(lon)) then
1121	call out_file%define_variable("lon", &
1122	& dim1_name="column", units_str="degrees_east", long_name="Longitude")
1123	end if
1124	call out_file%define_variable("skin_temperature", &
1125	& dim1_name="column", units_str="K", long_name="Skin_temperature")
1126	if (config%do_sw) then
1127	call out_file%define_variable("cos_solar_zenith_angle", &
1128	& dim1_name="column", units_str="1", &
1129	& long_name="Cosine of the solar zenith angle")
1130	end if
1131
1132	if (allocated(single_level%sw_albedo_direct)) then
1133	call out_file%define_variable("sw_albedo", &
1134	& dim2_name="column", dim1_name="sw_albedo_band", &
1135	& units_str="1", long_name="Shortwave surface albedo to diffuse radiation")
1136	call out_file%define_variable("sw_albedo_direct", &
1137	& dim2_name="column", dim1_name="sw_albedo_band", &
1138	& units_str="1", long_name="Shortwave surface albedo to direct radiation")
1139	else
1140	call out_file%define_variable("sw_albedo", &
1141	& dim2_name="column", dim1_name="sw_albedo_band", &
1142	& units_str="1", long_name="Shortwave surface albedo")
1143
1144	end if
1145	call out_file%define_variable("lw_emissivity", &
1146	& dim2_name="column", dim1_name="lw_emissivity_band", &
1147	& units_str="1", long_name="Longwave surface emissivity")
1148
1149	if (allocated(single_level%iseed)) then
1150	call out_file%define_variable("iseed", &
1151	& dim1_name="column", units_str="1", is_double=.true., &
1152	& long_name="Seed for random-number generator")
1153	end if
1154
1155	! Define thermodynamic variables on half levels
1156	call out_file%define_variable("pressure_hl", &
1157	& dim2_name="column", dim1_name="half_level", &
1158	& units_str="Pa", long_name="Pressure")
1159	call out_file%define_variable("temperature_hl", &
1160	& dim2_name="column", dim1_name="half_level", &
1161	& units_str="K", long_name="Temperature")
1162
1163	! Define gas mixing ratios on full levels
1164	call out_file%define_variable("q", &
1165	& dim2_name="column", dim1_name="level", &
1166	& units_str="1", long_name="Specific humidity")
1167	call out_file%define_variable("o3_mmr", &
1168	& dim2_name="column", dim1_name="level", &
1169	& units_str="1", long_name="Ozone mass mixing ratio")
1170	do jgas = 1,NMaxGases
1171	if (gas%is_present(jgas) .and. jgas /= IH2O .and. jgas /= IO3) then
1172	write(var_name,'(a,a)') trim(GasLowerCaseName(jgas)), '_vmr'
1173	write(long_name,'(a,a)') trim(GasName(jgas)), ' volume mixing ratio'
1174	call out_file%define_variable(trim(var_name), &
1175	& dim2_name="column", dim1_name="level", &
1176	& units_str="1", long_name=trim(long_name))
1177	end if
1178	end do
1179
1180	if (config%do_clouds) then
1181	! Define cloud variables on full levels
1182	call out_file%define_variable("cloud_fraction", &
1183	& dim2_name="column", dim1_name="level", &
1184	& units_str="1", long_name="Cloud fraction")
1185	call out_file%define_variable("q_liquid", &
1186	& dim2_name="column", dim1_name="level", &
1187	& units_str="1", long_name="Gridbox-mean liquid water mixing ratio")
1188	call out_file%define_variable("q_ice", &
1189	& dim2_name="column", dim1_name="level", &
1190	& units_str="1", long_name="Gridbox-mean ice water mixing ratio")
1191	call out_file%define_variable("re_liquid", &
1192	& dim2_name="column", dim1_name="level", &
1193	& units_str="m", long_name="Ice effective radius")
1194	if (associated(cloud%re_ice)) then
1195	call out_file%define_variable("re_ice", &
1196	& dim2_name="column", dim1_name="level", &
1197	& units_str="m", long_name="Ice effective radius")
1198	end if
1199	if (allocated(cloud%overlap_param)) then
1200	call out_file%define_variable("overlap_param", &
1201	& dim2_name="column", dim1_name="level_interface", &
1202	& units_str="1", long_name="Cloud overlap parameter")
1203	end if
1204	if (allocated(cloud%fractional_std)) then
1205	call out_file%define_variable("fractional_std", &
1206	& dim2_name="column", dim1_name="level", units_str="1", &
1207	& long_name="Fractional standard deviation of cloud optical depth")
1208	end if
1209	if (allocated(cloud%inv_cloud_effective_size)) then
1210	call out_file%define_variable("inv_cloud_effective_size", &
1211	& dim2_name="column", dim1_name="level", units_str="m-1", &
1212	& long_name="Inverse of cloud effective horizontal size")
1213	end if
1214	if (allocated(cloud%inv_inhom_effective_size)) then
1215	call out_file%define_variable("inv_inhom_effective_size", &
1216	& dim2_name="column", dim1_name="level", units_str="m-1", &
1217	& long_name="Inverse of cloud inhomogeneity effective horizontal size")
1218	end if
1219	end if ! do_clouds
1220
1221	! Define aerosol mass mixing ratio
1222	if (do_aerosol) then
1223	call out_file%define_variable("aerosol_mmr", &
1224	& dim3_name="column", dim2_name="aerosol_type", &
1225	& dim1_name="level", units_str="kg kg-1", &
1226	& long_name="Aerosol mass mixing ratio")
1227	end if
1228
1229	! Write variables
1230	call out_file%put("solar_irradiance", single_level%solar_irradiance)
1231	if (present(lat)) then
1232	call out_file%put("lat", lat)
1233	end if
1234	if (present(lon)) then
1235	call out_file%put("lon", lon)
1236	end if
1237	call out_file%put("skin_temperature", single_level%skin_temperature)
1238	if (config%do_sw) then
1239	call out_file%put("cos_solar_zenith_angle", single_level%cos_sza)
1240	end if
1241	call out_file%put("sw_albedo", single_level%sw_albedo)
1242	if (allocated(single_level%sw_albedo_direct)) then
1243	call out_file%put("sw_albedo_direct", single_level%sw_albedo_direct)
1244	end if
1245	call out_file%put("lw_emissivity", single_level%lw_emissivity)
1246	if (config%do_clouds .and. allocated(single_level%iseed)) then
1247	allocate(seed(ncol))
1248	seed = single_level%iseed
1249	call out_file%put("iseed", seed)
1250	deallocate(seed)
1251	end if
1252
1253	call out_file%put("pressure_hl", thermodynamics%pressure_hl)
1254	call out_file%put("temperature_hl", thermodynamics%temperature_hl)
1255
1256	allocate(mixing_ratio(ncol,nlev))
1257	call gas%get(IH2O, IMassMixingRatio, mixing_ratio)
1258	call out_file%put("q", mixing_ratio)
1259	call gas%get(IO3, IMassMixingRatio, mixing_ratio)
1260	call out_file%put("o3_mmr", mixing_ratio)
1261	do jgas = 1,NMaxGases
1262	if (gas%is_present(jgas) .and. jgas /= IH2O .and. jgas /= IO3) then
1263	write(var_name,'(a,a)') trim(GasLowerCaseName(jgas)), '_vmr'
1264	call gas%get(jgas, IVolumeMixingRatio, mixing_ratio)
1265	call out_file%put(trim(var_name), mixing_ratio)
1266	end if
1267	end do
1268	deallocate(mixing_ratio)
1269
1270	if (config%do_clouds) then
1271	call out_file%put("cloud_fraction", cloud%fraction)
1272	call out_file%put("q_liquid", cloud%q_liq)
1273	call out_file%put("q_ice", cloud%q_ice)
1274	call out_file%put("re_liquid", cloud%re_liq)
1275	if (associated(cloud%re_ice)) then
1276	call out_file%put("re_ice", cloud%re_ice)
1277	end if
1278	if (allocated(cloud%overlap_param)) then
1279	call out_file%put("overlap_param", cloud%overlap_param)
1280	end if
1281	if (allocated(cloud%fractional_std)) then
1282	call out_file%put("fractional_std", cloud%fractional_std)
1283	end if
1284	if (allocated(cloud%inv_cloud_effective_size)) then
1285	call out_file%put("inv_cloud_effective_size", cloud%inv_cloud_effective_size)
1286	end if
1287	if (allocated(cloud%inv_inhom_effective_size)) then
1288	call out_file%put("inv_inhom_effective_size", cloud%inv_inhom_effective_size)
1289	end if
1290	end if
1291
1292	if (do_aerosol) then
1293	allocate(aerosol_mmr(ncol, nlev, size(aerosol%mixing_ratio,3)))
1294	aerosol_mmr = 0.0_jprb
1295	aerosol_mmr(:,aerosol%istartlev:aerosol%iendlev,:) = aerosol%mixing_ratio
1296	call out_file%put("aerosol_mmr", aerosol_mmr)
1297	deallocate(aerosol_mmr)
1298	end if
1299
1300	! Close the file
1301	call out_file%close()
1302
1303	if (lhook) call dr_hook('radiation_save:save_inputs',1,hook_handle)
1304
1305	end subroutine save_inputs
1306
1307
1308	!---------------------------------------------------------------------
1309	! Save shortwave diagnostics computed from "flux" to NetCDF
1310	! file_name. The "mapping" matrix maps from fluxes in bands or
1311	! g-points to user specified spectral intervals, and should have
1312	! been produced by config%get_sw_mapping. See the example in
1313	! ecrad_driver.F90.
1314	subroutine save_sw_diagnostics(file_name, config, wavelength_bound, mapping, flux, &
1315	& iverbose, is_hdf5_file, experiment_name, &
1316	& is_double_precision)
1317
1318	use yomhook, only : lhook, dr_hook, jphook
1319
1320	use easy_netcdf
1321
1322	use radiation_io, only : nulout
1323	use radiation_flux, only : flux_type
1324	use radiation_config, only : config_type
1325	use radiation_matrix, only : sparse_x_dense
1326
1327	character(len=*), intent(in) :: file_name
1328	type(config_type), intent(in) :: config
1329	real(jprb), intent(in) :: wavelength_bound(:) ! m
1330	real(jprb), intent(in) :: mapping(:,:)
1331	type(flux_type), intent(in) :: flux
1332	integer, optional, intent(in) :: iverbose
1333	logical, optional, intent(in) :: is_hdf5_file
1334	logical, optional, intent(in) :: is_double_precision
1335	character(len=*), optional, intent(in) :: experiment_name
1336
1337	integer :: nwav ! Number of wavelength intervals
1338	real(jprb), allocatable :: flux_out(:,:)
1339	type(netcdf_file) :: out_file
1340	integer :: ncol, n_lev_plus1
1341	integer :: i_local_verbose
1342
1343	real(jphook) :: hook_handle
1344
1345	if (lhook) call dr_hook('radiation_save:save_sw_diagnostics',0,hook_handle)
1346
1347	if (present(iverbose)) then
1348	i_local_verbose = iverbose
1349	else
1350	i_local_verbose = config%iverbose
1351	end if
1352
1353	! Open the file
1354	call out_file%create(trim(file_name), iverbose=i_local_verbose, is_hdf5_file=is_hdf5_file)
1355
1356	! Set default precision for file, if specified
1357	if (present(is_double_precision)) then
1358	call out_file%double_precision(is_double_precision)
1359	end if
1360
1361	! Define dimensions
1362	ncol = size(flux%sw_up,1)
1363	call out_file%define_dimension("column", ncol)
1364	nwav = size(wavelength_bound)-1
1365	call out_file%define_dimension("wavelength", nwav)
1366
1367	! Put global attributes
1368	call out_file%put_global_attributes( &
1369	& title_str="Shortwave spectral diagnostics from the ecRad offline radiation model", &
1370	& references_str="Hogan, R. J., and A. Bozzo, 2018: A flexible and efficient radiation " &
1371	& //"scheme for the ECMWF model. J. Adv. Modeling Earth Sys., 10, 1990–2008", &
1372	& source_str="ecRad offline radiation model")
1373
1374	! Save "experiment" global attribute if present and not empty
1375	if (present(experiment_name)) then
1376	if (experiment_name /= " ") then
1377	call out_file%put_global_attribute("experiment", experiment_name)
1378	end if
1379	end if
1380
1381	! Define variables
1382	call out_file%define_variable("wavelength1", dim1_name="wavelength", &
1383	units_str="m", long_name="Wavelength lower bound")
1384	call out_file%define_variable("wavelength2", dim1_name="wavelength", &
1385	units_str="m", long_name="Wavelength upper bound")
1386
1387	! These are always available if
1388	! config%do_surface_sw_spectral_flux=true, which is required for
1389	! this routine
1390	call out_file%define_variable("flux_dn_sw_surf", &
1391	& dim2_name="column", dim1_name="wavelength", &
1392	& units_str="W m-2", long_name="Surface downwelling shortwave flux")
1393	call out_file%define_variable("flux_dn_direct_sw_surf", &
1394	& dim2_name="column", dim1_name="wavelength", &
1395	& units_str="W m-2", long_name="Surface downwelling direct shortwave flux")
1396	if (config%do_clear) then
1397	call out_file%define_variable("flux_dn_sw_surf_clear", &
1398	& dim2_name="column", dim1_name="wavelength", &
1399	& units_str="W m-2", long_name="Surface downwelling clear-sky shortwave flux")
1400	call out_file%define_variable("flux_dn_direct_sw_surf_clear", &
1401	& dim2_name="column", dim1_name="wavelength", &
1402	& units_str="W m-2", long_name="Surface downwelling clear-sky direct shortwave flux")
1403	end if
1404
1405	! The following are only availble if
1406	! config%do_save_spectral_flux=true
1407	if (allocated(flux%sw_up_band)) then
1408	call out_file%define_variable("flux_up_sw_toa", &
1409	& dim2_name="column", dim1_name="wavelength", &
1410	& units_str="W m-2", long_name="Top-of-atmosphere upwelling shortwave flux")
1411	call out_file%define_variable("flux_dn_sw_toa", &
1412	& dim2_name="column", dim1_name="wavelength", &
1413	& units_str="W m-2", long_name="Top-of-atmosphere downwelling shortwave flux")
1414	call out_file%define_variable("flux_up_sw_surf", &
1415	& dim2_name="column", dim1_name="wavelength", &
1416	& units_str="W m-2", long_name="Surface upwelling shortwave flux")
1417	if (allocated(flux%sw_up_clear_band)) then
1418	call out_file%define_variable("flux_up_sw_toa_clear", &
1419	& dim2_name="column", dim1_name="wavelength", &
1420	& units_str="W m-2", long_name="Top-of-atmosphere upwelling clear-sky shortwave flux")
1421	call out_file%define_variable("flux_up_sw_surf_clear", &
1422	& dim2_name="column", dim1_name="wavelength", &
1423	& units_str="W m-2", long_name="Surface upwelling clear-sky shortwave flux")
1424	end if
1425	end if
1426
1427	! Write variables
1428	call out_file%put("wavelength1", wavelength_bound(1:nwav))
1429	call out_file%put("wavelength2", wavelength_bound(2:nwav+1))
1430
1431	n_lev_plus1 = size(flux%sw_up,2)
1432
1433	! The mapping matrix is usually sparse, in which case we can check
1434	! its elements before multiplying a column of bandwise fluxes by
1435	! it.
1436	#define USE_SPARSE_MATMUL 1
1437	#ifdef USE_SPARSE_MATMUL
1438	#define MY_MATMUL sparse_x_dense
1439	#else
1440	#define MY_MATMUL matmul
1441	#endif
1442
1443	flux_out = MY_MATMUL(mapping, flux%sw_dn_surf_band)
1444	call out_file%put("flux_dn_sw_surf", flux_out)
1445	flux_out = MY_MATMUL(mapping, flux%sw_dn_direct_surf_band)
1446	call out_file%put("flux_dn_direct_sw_surf", flux_out)
1447	if (config%do_clear) then
1448	flux_out = MY_MATMUL(mapping, flux%sw_dn_surf_clear_band)
1449	call out_file%put("flux_dn_sw_surf_clear", flux_out)
1450	flux_out = MY_MATMUL(mapping, flux%sw_dn_direct_surf_clear_band)
1451	call out_file%put("flux_dn_direct_sw_surf_clear", flux_out)
1452	end if
1453
1454	if (allocated(flux%sw_up_band)) then
1455	flux_out = MY_MATMUL(mapping, flux%sw_up_band(:,:,n_lev_plus1))
1456	call out_file%put("flux_up_sw_surf", flux_out)
1457	flux_out = MY_MATMUL(mapping, flux%sw_up_band(:,:,1))
1458	call out_file%put("flux_up_sw_toa", flux_out)
1459	flux_out = MY_MATMUL(mapping, flux%sw_dn_band(:,:,1))
1460	call out_file%put("flux_dn_sw_toa", flux_out)
1461	if (allocated(flux%sw_up_clear_band)) then
1462	flux_out = MY_MATMUL(mapping, flux%sw_up_clear_band(:,:,1))
1463	call out_file%put("flux_up_sw_toa_clear", flux_out)
1464	flux_out = MY_MATMUL(mapping, flux%sw_up_clear_band(:,:,n_lev_plus1))
1465	call out_file%put("flux_up_sw_surf_clear", flux_out)
1466	end if
1467	end if
1468
1469	call out_file%close()
1470
1471	if (lhook) call dr_hook('radiation_save:save_sw_diagnostics',1,hook_handle)
1472
1473	end subroutine save_sw_diagnostics
1474
1475	end module radiation_save

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

Download in other formats:

Original Format