Context Navigation

radiation_flux.F90 @ 5418

Last change on this file since 5418 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: 36.1 KB

Rev	Line
[4773]	1	! radiation_flux.F90 - Derived type to store the output fluxes
	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-09-08 R. Hogan Store g-point fluxes
	17	! 2017-10-23 R. Hogan Renamed single-character variables
	18	! 2019-01-08 R. Hogan Added "indexed_sum_profile"
	19	! 2019-01-14 R. Hogan out_of_physical_bounds calls routine in radiation_config
	20	! 2021-01-20 R. Hogan Added heating_rate_out_of_physical_bounds function
	21	! 2022-12-07 R. Hogan Added top-of-atmosphere spectral output
	22
[4853]	23	#include "ecrad_config.h"
	24
[4773]	25	module radiation_flux
	26
	27	use parkind1, only : jprb
	28
	29	implicit none
	30	public
	31
	32	!---------------------------------------------------------------------
	33	! This derived type contains the output from the radiation
	34	! calculation. Currently this is solely flux profiles, but in
	35	! future surface fluxes in each band may be stored in order that the
	36	! calling program can compute surface-radiation such as
	37	! photosynthetically active radiation and UV index.
	38	type flux_type
	39	! All the following are broad-band fluxes in W m-2 with
	40	! dimensions (ncol,nlev+1). Note that only those fluxes that are
	41	! requested will be used, so clear-sky and direct-beam arrays may
	42	! not be allocated
	43	real(jprb), allocatable, dimension(:,:) :: &
	44	& lw_up, lw_dn, & ! Upwelling and downwelling longwave
	45	& sw_up, sw_dn, & ! Upwelling and downwelling shortwave
	46	& sw_dn_direct, & ! Direct-beam shortwave into a horizontal plane
	47	& lw_up_clear, lw_dn_clear, & ! Clear-sky quantities...
	48	& sw_up_clear, sw_dn_clear, &
	49	& sw_dn_direct_clear
	50	! As above but fluxes in each spectral band in W m-2 with
	51	! dimensions (nband,ncol,nlev+1). These are only allocated if
	52	! config%do_save_spectral_flux==.true.
	53	real(jprb), allocatable, dimension(:,:,:) :: &
	54	& lw_up_band, lw_dn_band, & ! Upwelling and downwelling longwave
	55	& sw_up_band, sw_dn_band, & ! Upwelling and downwelling shortwave
	56	& sw_dn_direct_band, & ! Direct-beam shortwave
	57	& lw_up_clear_band, lw_dn_clear_band, & ! Clear-sky quantities...
	58	& sw_up_clear_band, sw_dn_clear_band, &
	59	& sw_dn_direct_clear_band
	60	! Surface downwelling quantities at each g point, dimensioned
	61	! (ng,ncol), that are always saved by the solver, except for the
	62	! clear-sky ones that are only produced if
	63	! config%do_clear==.true.
	64	real(jprb), allocatable, dimension(:,:) :: &
	65	& lw_dn_surf_g, lw_dn_surf_clear_g, &
	66	& sw_dn_diffuse_surf_g, sw_dn_direct_surf_g, &
	67	& sw_dn_diffuse_surf_clear_g, sw_dn_direct_surf_clear_g
	68	! Top-of-atmosphere quantities at each g point, dimensioned
	69	! (ng,ncol), that are always saved by the solver, except for the
	70	! clear-sky ones that are only produced if
	71	! config%do_clear==.true.
	72	real(jprb), allocatable, dimension(:,:) :: &
	73	& lw_up_toa_g, lw_up_toa_clear_g, &
	74	& sw_dn_toa_g, sw_up_toa_g, sw_up_toa_clear_g
	75	! Shortwave downwelling spectral fluxes in W m-2 at the surface,
	76	! from which quantities such as photosynthetically active and UV
	77	! radiation can be computed. Only allocated if
	78	! config%do_surface_sw_spectral_flux==.true. Note that the
	79	! clear-sky quantities are only computed if
	80	! config%do_clear==.true., but direct fluxes are computed whether
	81	! or not do_direct==.true.. The dimensions are (nband,ncol).
	82	real(jprb), allocatable, dimension(:,:) :: &
	83	& sw_dn_surf_band, sw_dn_direct_surf_band, &
	84	& sw_dn_surf_clear_band, sw_dn_direct_surf_clear_band
	85	! Top-of-atmosphere spectral fluxes in W m-2. Only allocated if
	86	! config%do_toa_spectral_flux=.true.. Note that the clear-sky
	87	! quantities are only computed if config%do_clear==.true.. The
	88	! dimensions are (nband,ncol).
	89	real(jprb), allocatable, dimension(:,:) :: &
	90	& lw_up_toa_band, lw_up_toa_clear_band, &
	91	& sw_dn_toa_band, sw_up_toa_band, sw_up_toa_clear_band
	92	! Surface downwelling fluxes in W m-2 at the spectral resolution
	93	! needed by any subsequent canopy radiative transfer. If
	94	! config%use_canopy_full_spectrum_[sw\|lw] then these will be at
	95	! g-point resolution; otherwise they will be at
	96	! config%n_albedo_bands and config%n_emiss_bands resolution.
	97	real(jprb), allocatable, dimension(:,:) :: &
	98	& lw_dn_surf_canopy, &
	99	& sw_dn_diffuse_surf_canopy, sw_dn_direct_surf_canopy
	100
	101	! Diagnosed cloud cover from the short- and long-wave solvers
	102	real(jprb), allocatable, dimension(:) :: &
	103	& cloud_cover_lw, cloud_cover_sw
	104	! Longwave derivatives needed by Hogan and Bozzo (2015) method
	105	! for approximate longwave updates in between the full radiation
	106	! calls: rate of change of upwelling broad-band flux with respect
	107	! to surface value, dimensioned (ncol,nlev+1)
	108	real(jprb), allocatable, dimension(:,:) :: &
	109	& lw_derivatives
	110
	111	contains
	112	procedure :: allocate => allocate_flux_type
	113	procedure :: deallocate => deallocate_flux_type
	114	procedure :: calc_surface_spectral
	115	procedure :: calc_toa_spectral
	116	procedure :: out_of_physical_bounds
	117	procedure :: heating_rate_out_of_physical_bounds
	118	end type flux_type
	119
	120	! Added for DWD (2020)
[4853]	121	#ifdef DWD_VECTOR_OPTIMIZATIONS
[4773]	122	logical, parameter :: use_indexed_sum_vec = .true.
	123	#else
	124	logical, parameter :: use_indexed_sum_vec = .false.
	125	#endif
	126
	127	contains
	128
	129	!---------------------------------------------------------------------
	130	! Allocate arrays for flux profiles, using config to define which
	131	! fluxes are needed. The arrays are dimensioned for columns between
	132	! istartcol, iendcol and levels from 1 to nlev+1
	133	subroutine allocate_flux_type(this, config, istartcol, iendcol, nlev)
	134
	135	use yomhook, only : lhook, dr_hook, jphook
	136	use radiation_io, only : nulerr, radiation_abort
	137	use radiation_config, only : config_type
	138
	139	integer, intent(in) :: istartcol, iendcol, nlev
	140	class(flux_type), intent(inout) :: this
	141	type(config_type), intent(in) :: config
	142
	143	real(jphook) :: hook_handle
	144
	145	if (lhook) call dr_hook('radiation_flux:allocate',0,hook_handle)
	146
	147	! Allocate longwave arrays
	148	if (config%do_lw) then
	149	allocate(this%lw_up(istartcol:iendcol,nlev+1))
	150	allocate(this%lw_dn(istartcol:iendcol,nlev+1))
	151	if (config%do_clear) then
	152	allocate(this%lw_up_clear(istartcol:iendcol,nlev+1))
	153	allocate(this%lw_dn_clear(istartcol:iendcol,nlev+1))
	154	end if
	155
	156	if (config%do_save_spectral_flux) then
	157	if (config%n_spec_lw == 0) then
	158	write(nulerr,'(a)') '*** Error: number of LW spectral points to save not yet defined ' &
	159	& // 'so cannot allocate spectral flux arrays'
	160	call radiation_abort()
	161	end if
	162
	163	allocate(this%lw_up_band(config%n_spec_lw,istartcol:iendcol,nlev+1))
	164	allocate(this%lw_dn_band(config%n_spec_lw,istartcol:iendcol,nlev+1))
	165	if (config%do_clear) then
	166	allocate(this%lw_up_clear_band(config%n_spec_lw, &
	167	& istartcol:iendcol,nlev+1))
	168	allocate(this%lw_dn_clear_band(config%n_spec_lw, &
	169	& istartcol:iendcol,nlev+1))
	170	end if
	171	end if
	172
	173	if (config%do_lw_derivatives) then
	174	allocate(this%lw_derivatives(istartcol:iendcol,nlev+1))
	175	end if
	176
	177	if (config%do_toa_spectral_flux) then
	178	if (config%n_bands_lw == 0) then
	179	write(nulerr,'(a)') '*** Error: number of LW bands not yet defined ' &
	180	& // 'so cannot allocate TOA spectral flux arrays'
	181	call radiation_abort()
	182	end if
	183	allocate(this%lw_up_toa_band(config%n_bands_lw, istartcol:iendcol))
	184	if (config%do_clear) then
	185	allocate(this%lw_up_toa_clear_band(config%n_bands_lw, istartcol:iendcol))
	186	end if
	187	end if
	188
	189	! Allocate g-point downwelling fluxes at surface, and TOA fluxes
	190	allocate(this%lw_dn_surf_g(config%n_g_lw,istartcol:iendcol))
	191	allocate(this%lw_up_toa_g (config%n_g_lw,istartcol:iendcol))
	192	if (config%do_clear) then
	193	allocate(this%lw_dn_surf_clear_g(config%n_g_lw,istartcol:iendcol))
	194	allocate(this%lw_up_toa_clear_g (config%n_g_lw,istartcol:iendcol))
	195	end if
	196
	197	if (config%do_canopy_fluxes_lw) then
	198	! Downward fluxes at top of canopy at the spectral resolution
	199	! used in the canopy radiative transfer scheme
	200	allocate(this%lw_dn_surf_canopy(config%n_canopy_bands_lw,istartcol:iendcol))
	201	end if
	202	end if
	203
	204	! Allocate shortwave arrays
	205	if (config%do_sw) then
	206	allocate(this%sw_up(istartcol:iendcol,nlev+1))
	207	allocate(this%sw_dn(istartcol:iendcol,nlev+1))
	208	if (config%do_sw_direct) then
	209	allocate(this%sw_dn_direct(istartcol:iendcol,nlev+1))
	210	end if
	211	if (config%do_clear) then
	212	allocate(this%sw_up_clear(istartcol:iendcol,nlev+1))
	213	allocate(this%sw_dn_clear(istartcol:iendcol,nlev+1))
	214	if (config%do_sw_direct) then
	215	allocate(this%sw_dn_direct_clear(istartcol:iendcol,nlev+1))
	216	end if
	217	end if
	218
	219	if (config%do_save_spectral_flux) then
	220	if (config%n_spec_sw == 0) then
	221	write(nulerr,'(a)') '*** Error: number of SW spectral points to save not yet defined ' &
	222	& // 'so cannot allocate spectral flux arrays'
	223	call radiation_abort()
	224	end if
	225
	226	allocate(this%sw_up_band(config%n_spec_sw,istartcol:iendcol,nlev+1))
	227	allocate(this%sw_dn_band(config%n_spec_sw,istartcol:iendcol,nlev+1))
	228
	229	if (config%do_sw_direct) then
	230	allocate(this%sw_dn_direct_band(config%n_spec_sw, &
	231	& istartcol:iendcol,nlev+1))
	232	end if
	233	if (config%do_clear) then
	234	allocate(this%sw_up_clear_band(config%n_spec_sw, &
	235	& istartcol:iendcol,nlev+1))
	236	allocate(this%sw_dn_clear_band(config%n_spec_sw, &
	237	& istartcol:iendcol,nlev+1))
	238	if (config%do_sw_direct) then
	239	allocate(this%sw_dn_direct_clear_band(config%n_spec_sw, &
	240	& istartcol:iendcol, nlev+1))
	241	end if
	242	end if
	243	end if
	244
	245	if (config%do_surface_sw_spectral_flux) then
	246	if (config%n_bands_sw == 0) then
	247	write(nulerr,'(a)') '*** Error: number of SW bands not yet defined ' &
	248	& // 'so cannot allocate TOA spectral flux arrays'
	249	call radiation_abort()
	250	end if
	251	allocate(this%sw_dn_surf_band(config%n_bands_sw,istartcol:iendcol))
	252	allocate(this%sw_dn_direct_surf_band(config%n_bands_sw,istartcol:iendcol))
	253	if (config%do_clear) then
	254	allocate(this%sw_dn_surf_clear_band(config%n_bands_sw, &
	255	& istartcol:iendcol))
	256	allocate(this%sw_dn_direct_surf_clear_band(config%n_bands_sw, &
	257	& istartcol:iendcol))
	258	end if
	259	end if
	260
	261	if (config%do_toa_spectral_flux) then
	262	if (config%n_bands_sw == 0) then
	263	write(nulerr,'(a)') '*** Error: number of SW bands not yet defined ' &
	264	& // 'so cannot allocate surface spectral flux arrays'
	265	call radiation_abort()
	266	end if
	267	allocate(this%sw_dn_toa_band(config%n_bands_sw, istartcol:iendcol))
	268	allocate(this%sw_up_toa_band(config%n_bands_sw, istartcol:iendcol))
	269	if (config%do_clear) then
	270	allocate(this%sw_up_toa_clear_band(config%n_bands_sw, istartcol:iendcol))
	271	end if
	272	end if
	273
	274	! Allocate g-point downwelling fluxes at surface, and TOA fluxes
	275	allocate(this%sw_dn_diffuse_surf_g(config%n_g_sw,istartcol:iendcol))
	276	allocate(this%sw_dn_direct_surf_g (config%n_g_sw,istartcol:iendcol))
	277	allocate(this%sw_dn_toa_g (config%n_g_sw,istartcol:iendcol))
	278	allocate(this%sw_up_toa_g (config%n_g_sw,istartcol:iendcol))
	279	if (config%do_clear) then
	280	allocate(this%sw_dn_diffuse_surf_clear_g(config%n_g_sw,istartcol:iendcol))
	281	allocate(this%sw_dn_direct_surf_clear_g (config%n_g_sw,istartcol:iendcol))
	282	allocate(this%sw_up_toa_clear_g (config%n_g_sw,istartcol:iendcol))
	283	end if
	284
	285	if (config%do_canopy_fluxes_sw) then
	286	! Downward fluxes at top of canopy at the spectral resolution
	287	! used in the canopy radiative transfer scheme
	288	allocate(this%sw_dn_diffuse_surf_canopy(config%n_canopy_bands_sw,istartcol:iendcol))
	289	allocate(this%sw_dn_direct_surf_canopy (config%n_canopy_bands_sw,istartcol:iendcol))
	290	end if
	291	end if
	292
	293	! Allocate cloud cover arrays
	294	allocate(this%cloud_cover_lw(istartcol:iendcol))
	295	allocate(this%cloud_cover_sw(istartcol:iendcol))
	296
	297	! Some solvers may not write to cloud cover, so we initialize to
	298	! an unphysical value
	299	this%cloud_cover_lw = -1.0_jprb
	300	this%cloud_cover_sw = -1.0_jprb
	301
	302	if (lhook) call dr_hook('radiation_flux:allocate',1,hook_handle)
	303
	304	end subroutine allocate_flux_type
	305
	306
	307	!---------------------------------------------------------------------
	308	! Deallocate flux arrays
	309	subroutine deallocate_flux_type(this)
	310
	311	use yomhook, only : lhook, dr_hook, jphook
	312
	313	class(flux_type), intent(inout) :: this
	314	real(jphook) :: hook_handle
	315
	316	if (lhook) call dr_hook('radiation_flux:deallocate',0,hook_handle)
	317
	318	if (allocated(this%lw_up)) then
	319	deallocate(this%lw_up)
	320	if (allocated(this%lw_dn)) deallocate(this%lw_dn)
	321	if (allocated(this%lw_up_clear)) deallocate(this%lw_up_clear)
	322	if (allocated(this%lw_dn_clear)) deallocate(this%lw_dn_clear)
	323	end if
	324
	325	if (allocated(this%sw_up)) then
	326	deallocate(this%sw_up)
	327	if (allocated(this%sw_dn)) deallocate(this%sw_dn)
	328	if (allocated(this%sw_up_clear)) deallocate(this%sw_up_clear)
	329	if (allocated(this%sw_dn_clear)) deallocate(this%sw_dn_clear)
	330	if (allocated(this%sw_dn_direct)) deallocate(this%sw_dn_direct)
	331	if (allocated(this%sw_dn_direct_clear)) &
	332	& deallocate(this%sw_dn_direct_clear)
	333	end if
	334
	335	if (allocated(this%lw_up_band)) then
	336	deallocate(this%lw_up_band)
	337	if (allocated(this%lw_dn_band)) deallocate(this%lw_dn_band)
	338	if (allocated(this%lw_up_clear_band)) deallocate(this%lw_up_clear_band)
	339	if (allocated(this%lw_dn_clear_band)) deallocate(this%lw_dn_clear_band)
	340	end if
	341
	342	if (allocated(this%sw_up_band)) then
	343	deallocate(this%sw_up_band)
	344	if (allocated(this%sw_dn_band)) deallocate(this%sw_dn_band)
	345	if (allocated(this%sw_up_clear_band)) deallocate(this%sw_up_clear_band)
	346	if (allocated(this%sw_dn_clear_band)) deallocate(this%sw_dn_clear_band)
	347	if (allocated(this%sw_dn_direct_band)) deallocate(this%sw_dn_direct_band)
	348	if (allocated(this%sw_dn_direct_clear_band)) &
	349	& deallocate(this%sw_dn_direct_clear_band)
	350	end if
	351
	352	if (allocated(this%sw_dn_surf_band)) then
	353	deallocate(this%sw_dn_surf_band)
	354	deallocate(this%sw_dn_direct_surf_band)
	355	end if
	356	if (allocated(this%sw_dn_surf_clear_band)) then
	357	deallocate(this%sw_dn_surf_clear_band)
	358	deallocate(this%sw_dn_direct_surf_clear_band)
	359	end if
	360
	361	if (allocated(this%lw_dn_surf_canopy)) deallocate(this%lw_dn_surf_canopy)
	362	if (allocated(this%sw_dn_diffuse_surf_canopy)) deallocate(this%sw_dn_diffuse_surf_canopy)
	363	if (allocated(this%sw_dn_direct_surf_canopy)) deallocate(this%sw_dn_direct_surf_canopy)
	364
	365	if (allocated(this%cloud_cover_sw)) then
	366	deallocate(this%cloud_cover_sw)
	367	end if
	368	if (allocated(this%cloud_cover_lw)) then
	369	deallocate(this%cloud_cover_lw)
	370	end if
	371
	372	if (allocated(this%lw_derivatives)) then
	373	deallocate(this%lw_derivatives)
	374	end if
	375
	376	if (allocated(this%lw_dn_surf_g)) deallocate(this%lw_dn_surf_g)
	377	if (allocated(this%lw_dn_surf_clear_g)) deallocate(this%lw_dn_surf_clear_g)
	378	if (allocated(this%sw_dn_diffuse_surf_g)) deallocate(this%sw_dn_diffuse_surf_g)
	379	if (allocated(this%sw_dn_direct_surf_g)) deallocate(this%sw_dn_direct_surf_g)
	380	if (allocated(this%sw_dn_diffuse_surf_clear_g)) deallocate(this%sw_dn_diffuse_surf_clear_g)
	381	if (allocated(this%sw_dn_direct_surf_clear_g)) deallocate(this%sw_dn_direct_surf_clear_g)
	382
	383	if (allocated(this%lw_up_toa_g)) deallocate(this%lw_up_toa_g)
	384	if (allocated(this%sw_up_toa_g)) deallocate(this%sw_up_toa_g)
	385	if (allocated(this%sw_dn_toa_g)) deallocate(this%sw_dn_toa_g)
	386	if (allocated(this%lw_up_toa_clear_g)) deallocate(this%lw_up_toa_clear_g)
	387	if (allocated(this%sw_up_toa_clear_g)) deallocate(this%sw_up_toa_clear_g)
	388
	389	if (lhook) call dr_hook('radiation_flux:deallocate',1,hook_handle)
	390
	391	end subroutine deallocate_flux_type
	392
	393
	394	!---------------------------------------------------------------------
	395	! Calculate surface downwelling fluxes in each band using the
	396	! downwelling surface fluxes at each g point
	397	subroutine calc_surface_spectral(this, config, istartcol, iendcol)
	398
	399	use yomhook, only : lhook, dr_hook, jphook
	400	use radiation_config, only : config_type
	401
	402	class(flux_type), intent(inout) :: this
	403	type(config_type), intent(in) :: config
	404	integer, intent(in) :: istartcol, iendcol
	405
	406	integer :: jcol, jband, jalbedoband, nalbedoband
	407
	408	! Longwave surface downwelling in each band needed to compute
	409	! canopy fluxes
	410	real(jprb) :: lw_dn_surf_band(config%n_bands_lw,istartcol:iendcol)
	411
	412	real(jphook) :: hook_handle
	413
	414	if (lhook) call dr_hook('radiation_flux:calc_surface_spectral',0,hook_handle)
	415
	416	if (config%do_sw .and. config%do_surface_sw_spectral_flux) then
	417
	418	if (use_indexed_sum_vec) then
	419	call indexed_sum_vec(this%sw_dn_direct_surf_g, &
	420	& config%i_band_from_reordered_g_sw, &
	421	& this%sw_dn_direct_surf_band, istartcol, iendcol)
	422	call indexed_sum_vec(this%sw_dn_diffuse_surf_g, &
	423	& config%i_band_from_reordered_g_sw, &
	424	& this%sw_dn_surf_band, istartcol, iendcol)
	425	do jcol = istartcol,iendcol
	426	this%sw_dn_surf_band(:,jcol) &
	427	& = this%sw_dn_surf_band(:,jcol) &
	428	& + this%sw_dn_direct_surf_band(:,jcol)
	429	end do
	430	else
	431	do jcol = istartcol,iendcol
	432	call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), &
	433	& config%i_band_from_reordered_g_sw, &
	434	& this%sw_dn_direct_surf_band(:,jcol))
	435	call indexed_sum(this%sw_dn_diffuse_surf_g(:,jcol), &
	436	& config%i_band_from_reordered_g_sw, &
	437	& this%sw_dn_surf_band(:,jcol))
	438	this%sw_dn_surf_band(:,jcol) &
	439	& = this%sw_dn_surf_band(:,jcol) &
	440	& + this%sw_dn_direct_surf_band(:,jcol)
	441	end do
	442	end if
	443
	444	if (config%do_clear) then
	445	if (use_indexed_sum_vec) then
	446	call indexed_sum_vec(this%sw_dn_direct_surf_clear_g, &
	447	& config%i_band_from_reordered_g_sw, &
	448	& this%sw_dn_direct_surf_clear_band, istartcol, iendcol)
	449	call indexed_sum_vec(this%sw_dn_diffuse_surf_clear_g, &
	450	& config%i_band_from_reordered_g_sw, &
	451	& this%sw_dn_surf_clear_band, istartcol, iendcol)
	452	do jcol = istartcol,iendcol
	453	this%sw_dn_surf_clear_band(:,jcol) &
	454	& = this%sw_dn_surf_clear_band(:,jcol) &
	455	& + this%sw_dn_direct_surf_clear_band(:,jcol)
	456	end do
	457	else
	458	do jcol = istartcol,iendcol
	459	call indexed_sum(this%sw_dn_direct_surf_clear_g(:,jcol), &
	460	& config%i_band_from_reordered_g_sw, &
	461	& this%sw_dn_direct_surf_clear_band(:,jcol))
	462	call indexed_sum(this%sw_dn_diffuse_surf_clear_g(:,jcol), &
	463	& config%i_band_from_reordered_g_sw, &
	464	& this%sw_dn_surf_clear_band(:,jcol))
	465	this%sw_dn_surf_clear_band(:,jcol) &
	466	& = this%sw_dn_surf_clear_band(:,jcol) &
	467	& + this%sw_dn_direct_surf_clear_band(:,jcol)
	468	end do
	469	end if
	470	end if
	471
	472	end if ! do_surface_sw_spectral_flux
	473
	474	! Fluxes in bands required for canopy radiative transfer
	475	if (config%do_sw .and. config%do_canopy_fluxes_sw) then
	476	if (config%use_canopy_full_spectrum_sw) then
	477	this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) = this%sw_dn_diffuse_surf_g(:,istartcol:iendcol)
	478	this%sw_dn_direct_surf_canopy (:,istartcol:iendcol) = this%sw_dn_direct_surf_g (:,istartcol:iendcol)
	479	else if (config%do_nearest_spectral_sw_albedo) then
	480	if (use_indexed_sum_vec) then
	481	call indexed_sum_vec(this%sw_dn_direct_surf_g, &
	482	& config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
	483	& this%sw_dn_direct_surf_canopy, istartcol, iendcol)
	484	call indexed_sum_vec(this%sw_dn_diffuse_surf_g, &
	485	& config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
	486	& this%sw_dn_diffuse_surf_canopy, istartcol, iendcol)
	487	else
	488	do jcol = istartcol,iendcol
	489	call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), &
	490	& config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
	491	& this%sw_dn_direct_surf_canopy(:,jcol))
	492	call indexed_sum(this%sw_dn_diffuse_surf_g(:,jcol), &
	493	& config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
	494	& this%sw_dn_diffuse_surf_canopy(:,jcol))
	495	end do
	496	end if
	497	else
	498	! More accurate calculations using weights, but requires
	499	! this%sw_dn_[direct_]surf_band to be defined, i.e.
	500	! config%do_surface_sw_spectral_flux == .true.
	501	nalbedoband = size(config%sw_albedo_weights,1)
	502	this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) = 0.0_jprb
	503	this%sw_dn_direct_surf_canopy (:,istartcol:iendcol) = 0.0_jprb
	504	do jband = 1,config%n_bands_sw
	505	do jalbedoband = 1,nalbedoband
	506	if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then
	507	! Initially, "diffuse" is actually "total"
	508	this%sw_dn_diffuse_surf_canopy(jalbedoband,istartcol:iendcol) &
	509	& = this%sw_dn_diffuse_surf_canopy(jalbedoband,istartcol:iendcol) &
	510	& + config%sw_albedo_weights(jalbedoband,jband) &
	511	& * this%sw_dn_surf_band(jband,istartcol:iendcol)
	512	this%sw_dn_direct_surf_canopy(jalbedoband,istartcol:iendcol) &
	513	& = this%sw_dn_direct_surf_canopy(jalbedoband,istartcol:iendcol) &
	514	& + config%sw_albedo_weights(jalbedoband,jband) &
	515	& * this%sw_dn_direct_surf_band(jband,istartcol:iendcol)
	516	end if
	517	end do
	518	end do
	519	! Subtract the direct from total to get diffuse
	520	this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) &
	521	& = this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) &
	522	& - this%sw_dn_direct_surf_canopy(:,istartcol:iendcol)
	523	end if
	524
	525	end if ! do_canopy_fluxes_sw
	526
	527	if (config%do_lw .and. config%do_canopy_fluxes_lw) then
	528	if (config%use_canopy_full_spectrum_lw) then
	529	this%lw_dn_surf_canopy(:,istartcol:iendcol) = this%lw_dn_surf_g(:,istartcol:iendcol)
	530	else if (config%do_nearest_spectral_lw_emiss) then
	531	if (use_indexed_sum_vec) then
	532	call indexed_sum_vec(this%lw_dn_surf_g, &
	533	& config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), &
	534	& this%lw_dn_surf_canopy, istartcol, iendcol)
	535	else
	536	do jcol = istartcol,iendcol
	537	call indexed_sum(this%lw_dn_surf_g(:,jcol), &
	538	& config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), &
	539	& this%lw_dn_surf_canopy(:,jcol))
	540	end do
	541	end if
	542	else
	543	! Compute fluxes in each longwave emissivity interval using
	544	! weights; first sum over g points to get the values in bands
	545	if (use_indexed_sum_vec) then
	546	call indexed_sum_vec(this%lw_dn_surf_g, &
	547	& config%i_band_from_reordered_g_lw, &
	548	& lw_dn_surf_band, istartcol, iendcol)
	549	else
	550	do jcol = istartcol,iendcol
	551	call indexed_sum(this%lw_dn_surf_g(:,jcol), &
	552	& config%i_band_from_reordered_g_lw, &
	553	& lw_dn_surf_band(:,jcol))
	554	end do
	555	end if
	556	nalbedoband = size(config%lw_emiss_weights,1)
	557	this%lw_dn_surf_canopy(:,istartcol:iendcol) = 0.0_jprb
	558	do jband = 1,config%n_bands_lw
	559	do jalbedoband = 1,nalbedoband
	560	if (config%lw_emiss_weights(jalbedoband,jband) /= 0.0_jprb) then
	561	this%lw_dn_surf_canopy(jalbedoband,istartcol:iendcol) &
	562	& = this%lw_dn_surf_canopy(jalbedoband,istartcol:iendcol) &
	563	& + config%lw_emiss_weights(jalbedoband,jband) &
	564	& * lw_dn_surf_band(jband,istartcol:iendcol)
	565	end if
	566	end do
	567	end do
	568	end if
	569	end if
	570
	571	if (lhook) call dr_hook('radiation_flux:calc_surface_spectral',1,hook_handle)
	572
	573	end subroutine calc_surface_spectral
	574
	575
	576	!---------------------------------------------------------------------
	577	! Calculate top-of-atmosphere fluxes in each band using the fluxes
	578	! at each g point
	579	subroutine calc_toa_spectral(this, config, istartcol, iendcol)
	580
	581	use yomhook, only : lhook, dr_hook, jphook
	582	use radiation_config, only : config_type
	583
	584	class(flux_type), intent(inout) :: this
	585	type(config_type), intent(in) :: config
	586	integer, intent(in) :: istartcol, iendcol
	587
	588	integer :: jcol, jband
	589
	590	real(jphook) :: hook_handle
	591
	592	if (lhook) call dr_hook('radiation_flux:calc_toa_spectral',0,hook_handle)
	593
	594	if (config%do_sw .and. config%do_toa_spectral_flux) then
	595
	596	if (use_indexed_sum_vec) then
	597	call indexed_sum_vec(this%sw_dn_toa_g, &
	598	& config%i_band_from_reordered_g_sw, &
	599	& this%sw_dn_toa_band, istartcol, iendcol)
	600	call indexed_sum_vec(this%sw_up_toa_g, &
	601	& config%i_band_from_reordered_g_sw, &
	602	& this%sw_up_toa_band, istartcol, iendcol)
	603	else
	604	do jcol = istartcol,iendcol
	605	call indexed_sum(this%sw_dn_toa_g(:,jcol), &
	606	& config%i_band_from_reordered_g_sw, &
	607	& this%sw_dn_toa_band(:,jcol))
	608	call indexed_sum(this%sw_up_toa_g(:,jcol), &
	609	& config%i_band_from_reordered_g_sw, &
	610	& this%sw_up_toa_band(:,jcol))
	611	end do
	612	end if
	613
	614	if (config%do_clear) then
	615	if (use_indexed_sum_vec) then
	616	call indexed_sum_vec(this%sw_up_toa_clear_g, &
	617	& config%i_band_from_reordered_g_sw, &
	618	& this%sw_up_toa_clear_band, istartcol, iendcol)
	619	else
	620	do jcol = istartcol,iendcol
	621	call indexed_sum(this%sw_up_toa_clear_g(:,jcol), &
	622	& config%i_band_from_reordered_g_sw, &
	623	& this%sw_up_toa_clear_band(:,jcol))
	624	end do
	625	end if
	626	end if
	627	end if
	628
	629	if (config%do_lw .and. config%do_toa_spectral_flux) then
	630
	631	if (use_indexed_sum_vec) then
	632	call indexed_sum_vec(this%lw_up_toa_g, &
	633	& config%i_band_from_reordered_g_lw, &
	634	& this%lw_up_toa_band, istartcol, iendcol)
	635	else
	636	do jcol = istartcol,iendcol
	637	call indexed_sum(this%lw_up_toa_g(:,jcol), &
	638	& config%i_band_from_reordered_g_lw, &
	639	& this%lw_up_toa_band(:,jcol))
	640	end do
	641	end if
	642
	643	if (config%do_clear) then
	644	if (use_indexed_sum_vec) then
	645	call indexed_sum_vec(this%lw_up_toa_clear_g, &
	646	& config%i_band_from_reordered_g_lw, &
	647	& this%lw_up_toa_clear_band, istartcol, iendcol)
	648	else
	649	do jcol = istartcol,iendcol
	650	call indexed_sum(this%lw_up_toa_clear_g(:,jcol), &
	651	& config%i_band_from_reordered_g_lw, &
	652	& this%lw_up_toa_clear_band(:,jcol))
	653	end do
	654	end if
	655	end if
	656	end if
	657
	658	if (lhook) call dr_hook('radiation_flux:calc_toa_spectral',1,hook_handle)
	659
	660	end subroutine calc_toa_spectral
	661
	662
	663	!---------------------------------------------------------------------
	664	! Return .true. if the most important flux variables are out of a
	665	! physically sensible range, optionally only considering columns
	666	! between istartcol and iendcol
	667	function out_of_physical_bounds(this, istartcol, iendcol) result(is_bad)
	668
	669	use yomhook, only : lhook, dr_hook, jphook
	670	use radiation_check, only : out_of_bounds_2d
	671
	672	class(flux_type), intent(inout) :: this
	673	integer, optional,intent(in) :: istartcol, iendcol
	674	logical :: is_bad
	675
	676	real(jphook) :: hook_handle
	677
	678	if (lhook) call dr_hook('radiation_flux:out_of_physical_bounds',0,hook_handle)
	679
	680	is_bad = out_of_bounds_2d(this%lw_up, 'lw_up', 10.0_jprb, 900.0_jprb, .false., i1=istartcol, i2=iendcol) &
	681	& .or. out_of_bounds_2d(this%lw_dn, 'lw_dn', 0.0_jprb, 800.0_jprb, .false., i1=istartcol, i2=iendcol) &
	682	& .or. out_of_bounds_2d(this%sw_up, 'sw_up', 0.0_jprb, 1500.0_jprb, .false., i1=istartcol, i2=iendcol) &
	683	& .or. out_of_bounds_2d(this%sw_dn, 'sw_dn', 0.0_jprb, 1500.0_jprb, .false., i1=istartcol, i2=iendcol) &
	684	& .or. out_of_bounds_2d(this%sw_dn_direct, 'sw_dn_direct', 0.0_jprb, 1500.0_jprb, .false., i1=istartcol, i2=iendcol) &
	685	& .or. out_of_bounds_2d(this%lw_derivatives, 'lw_derivatives', 0.0_jprb, 1.0_jprb, .false., i1=istartcol, i2=iendcol) &
	686	& .or. out_of_bounds_2d(this%sw_dn_surf_band, 'sw_dn_surf_band', 0.0_jprb, 1500.0_jprb, &
	687	& .false., j1=istartcol, j2=iendcol) &
	688	& .or. out_of_bounds_2d(this%sw_dn_surf_clear_band, 'sw_dn_surf_clear_band', 0.0_jprb, 1500.0_jprb, &
	689	& .false., j1=istartcol, j2=iendcol)
	690
	691	if (lhook) call dr_hook('radiation_flux:out_of_physical_bounds',1,hook_handle)
	692
	693	end function out_of_physical_bounds
	694
	695	!---------------------------------------------------------------------
	696	! Return .true. if the heating rates are out of a physically
	697	! sensible range, optionally only considering columns between
	698	! istartcol and iendcol. This function allocates and deallocates
	699	! memory due to the requirements for inputs of out_of_bounds_2d.
	700	function heating_rate_out_of_physical_bounds(this, nlev, istartcol, iendcol, pressure_hl) result(is_bad)
	701
	702	use radiation_check, only : out_of_bounds_2d
	703	use radiation_constants, only : AccelDueToGravity
	704
	705	! "Cp" (J kg-1 K-1)
	706	real(jprb), parameter :: SpecificHeatDryAir = 1004.0
	707
	708	class(flux_type), intent(inout) :: this
	709	integer, intent(in) :: istartcol, iendcol, nlev
	710	logical :: is_bad
	711
	712	real(jprb), intent(in) :: pressure_hl(:,:)
	713
	714	real(jprb), allocatable :: hr_K_day(:,:)
	715
	716	real(jprb) :: scaling(istartcol:iendcol,nlev)
	717
	718	allocate(hr_K_day(istartcol:iendcol,nlev))
	719
	720	scaling = -(24.0_jprb * 3600.0_jprb * AccelDueToGravity / SpecificHeatDryAir) &
	721	& / (pressure_hl(istartcol:iendcol,2:nlev+1) - pressure_hl(istartcol:iendcol,1:nlev))
	722	! Shortwave
	723	hr_K_day = scaling * (this%sw_dn(istartcol:iendcol,2:nlev+1) - this%sw_up(istartcol:iendcol,2:nlev+1) &
	724	& -this%sw_dn(istartcol:iendcol,1:nlev) + this%sw_up(istartcol:iendcol,1:nlev))
	725	is_bad = out_of_bounds_2d(hr_K_day, 'sw_heating_rate_K_day', 0.0_jprb, 200.0_jprb, &
	726	& .false., i1=istartcol, i2=iendcol)
	727
	728	! Longwave
	729	hr_K_day = scaling * (this%lw_dn(istartcol:iendcol,2:nlev+1) - this%lw_up(istartcol:iendcol,2:nlev+1) &
	730	& -this%lw_dn(istartcol:iendcol,1:nlev) + this%lw_up(istartcol:iendcol,1:nlev))
	731	is_bad = is_bad .or. out_of_bounds_2d(hr_K_day, 'lw_heating_rate_K_day', -250.0_jprb, 150.0_jprb, &
	732	& .false., i1=istartcol, i2=iendcol)
	733
	734	deallocate(hr_K_day)
	735
	736	end function heating_rate_out_of_physical_bounds
	737
	738
	739	!---------------------------------------------------------------------
	740	! Sum elements of "source" into "dest" according to index "ind".
	741	! "source" and "ind" should have the same size and bounds, and no
	742	! element of "ind" should refer outside the bounds of "dest". This
	743	! version increments existing contents of "dest".
	744	pure subroutine add_indexed_sum(source, ind, dest)
	745
	746	real(jprb), intent(in) :: source(:)
	747	integer, intent(in) :: ind(:)
	748	real(jprb), intent(inout) :: dest(:)
	749
	750	integer :: ig, jg, istart, iend
	751
	752	istart = lbound(source,1)
	753	iend = ubound(source,1)
	754
	755	do jg = istart, iend
	756	ig = ind(jg)
	757	dest(ig) = dest(ig) + source(jg)
	758	end do
	759
	760	end subroutine add_indexed_sum
	761
	762
	763	!---------------------------------------------------------------------
	764	! As "add_indexed_sum" but this version overwrites existing contents
	765	! of "dest"
	766	pure subroutine indexed_sum(source, ind, dest)
	767
	768	real(jprb), intent(in) :: source(:)
	769	integer, intent(in) :: ind(:)
	770	real(jprb), intent(out) :: dest(:)
	771
	772	integer :: ig, jg, istart, iend
	773
	774	dest = 0.0
	775
	776	istart = lbound(source,1)
	777	iend = ubound(source,1)
	778
	779	do jg = istart, iend
	780	ig = ind(jg)
	781	dest(ig) = dest(ig) + source(jg)
	782	end do
	783
	784	end subroutine indexed_sum
	785
	786	!---------------------------------------------------------------------
	787	! Vectorized version of "add_indexed_sum"
	788	subroutine indexed_sum_vec(source, ind, dest, ist, iend)
	789
	790	real(jprb), intent(in) :: source(:,:)
	791	integer, intent(in) :: ind(:)
	792	real(jprb), intent(out) :: dest(:,:)
	793	integer, intent(in) :: ist, iend
	794
	795	integer :: ig, jg, jc
	796
	797	dest = 0.0
	798
	799	do jg = lbound(source,1), ubound(source,1)
	800	ig = ind(jg)
	801	do jc = ist, iend
	802	dest(ig,jc) = dest(ig,jc) + source(jg,jc)
	803	end do
	804	end do
	805
	806	end subroutine indexed_sum_vec
	807
	808	!---------------------------------------------------------------------
	809	! As "add_indexed_sum" but a whole vertical profiles
	810	pure subroutine add_indexed_sum_profile(source, ind, dest)
	811
	812	real(jprb), intent(in) :: source(:,:)
	813	integer, intent(in) :: ind(:)
	814	real(jprb), intent(out) :: dest(:,:)
	815
	816	integer :: ig, jg, istart, iend, jlev, nlev
	817
	818	istart = lbound(source,1)
	819	iend = ubound(source,1)
	820	nlev = size(source,2)
	821
	822	do jlev = 1,nlev
	823	do jg = istart, iend
	824	ig = ind(jg)
	825	dest(ig,jlev) = dest(ig,jlev) + source(jg,jlev)
	826	end do
	827	end do
	828
	829	end subroutine add_indexed_sum_profile
	830
	831
	832	!---------------------------------------------------------------------
	833	! As "indexed_sum" but a whole vertical profiles
	834	pure subroutine indexed_sum_profile(source, ind, dest)
	835
	836	real(jprb), intent(in) :: source(:,:)
	837	integer, intent(in) :: ind(:)
	838	real(jprb), intent(out) :: dest(:,:)
	839
	840	integer :: ig, jg, istart, iend, jlev, nlev
	841
	842	dest = 0.0
	843
	844	istart = lbound(source,1)
	845	iend = ubound(source,1)
	846	nlev = size(source,2)
	847
	848	do jlev = 1,nlev
	849	do jg = istart, iend
	850	ig = ind(jg)
	851	dest(ig,jlev) = dest(ig,jlev) + source(jg,jlev)
	852	end do
	853	end do
	854
	855	end subroutine indexed_sum_profile
	856
	857	end module radiation_flux

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

Download in other formats:

Original Format