Context Navigation

radiation_save.F90 @ 5433

Last change on this file since 5433 was 4853, checked in by idelkadi, 9 months ago

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

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

File size: 63.5 KB

Rev	Line
[4773]	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
[4853]	88	if (config%i_gas_model_lw == IGasModelMonochromatic &
[4773]	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
[4853]	517	if (config%i_gas_model_lw == IGasModelMonochromatic &
[4773]	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