Context Navigation

radiation_save.F90 @ 5229

Last change on this file since 5229 was 4489, checked in by lguez, 20 months ago
Merge LMDZ_ECRad branch back into trunk!
File size: 43.3 KB

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

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

Download in other formats:

Original Format