Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

radiation_save.F90 @ 5407

Last change on this file since 5407 was 4728, checked in by idelkadi, 14 months ago

Update of ecrad in the LMDZ_ECRad branch of LMDZ:

version 1.6.1 of ecrad
files are no longer grouped in the same ecrad directory.
the structure of ecrad offline is preserved to facilitate updating in LMDZ
cfg.bld modified to take into account the new added subdirectories.
the interface routines and those added in ecrad are moved to the phylmd directory

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 == 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 == 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