Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

radsurf_intermediate.F90 @ 4013

Last change on this file since 4013 was 3908, checked in by idelkadi, 4 years ago

Online implementation of the radiative transfer code ECRAD in the LMDZ model.

Inclusion of the ecrad directory containing the sources of the ECRAD code
- interface routine : radiation_scheme.F90
Adaptation of compilation scripts :
- compilation under CPP key CPP_ECRAD
- compilation with option "-rad ecard" or "-ecard true"
- The "-rad old/rtm/ecran" build option will need to replace the "-rrtm true" and "-ecrad true" options in the future.
Runing LMDZ simulations with ecrad, you need :
- logical key iflag_rrtm = 2 in physiq.def
- namelist_ecrad (DefLists?)
- the directory "data" containing the configuration files is temporarily placed in ../libfphylmd/ecrad/
Compilation and execution are tested in the 1D case. The repository under svn would allow to continue the implementation work: tests, verification of the results, ...

File size: 65.1 KB

Line
1	! radsurf_intermediate.f90 - Derived type for intermediate radiative properties of surface
2	!
3	! (C) Copyright 2017- 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
16	module radsurf_intermediate
17
18	use parkind1, only : jpim, jprb
19
20	implicit none
21
22	public
23
24	!---------------------------------------------------------------------
25	! Derived type storing a description of radiative properties at the
26	! level of individual facets
27	type surface_intermediate_type
28
29	! Surface "facet" properties, dimensioned (ng, nfacet, istartcol:iendcol)
30
31	! Longwave blackbody emission (W m-2) and emissivity from facets
32	real(kind=jprb), allocatable, dimension(:,:,:) :: planck_facet
33	real(kind=jprb), allocatable, dimension(:,:,:) :: lw_emissivity
34	! Shortwave direct and diffuse albedo
35	real(kind=jprb), allocatable, dimension(:,:,:) :: sw_albedo_direct
36	real(kind=jprb), allocatable, dimension(:,:,:) :: sw_albedo_diffuse
37
38	! Longwave blackbody emission (W m-2) from regions (volumes
39	! e.g. vegetation and urban canopies)
40	real(kind=jprb), allocatable, dimension(:,:,:) :: planck_region
41
42	! Volumetric "region" properties of the canopy, e.g. vegetation or
43	! the space between buildings, dimensioned
44	! (ng,nregion,istartcol:iendcol)
45
46	! Shortwave reflectance and transmittance to diffuse radiation
47	real(kind=jprb), allocatable, dimension(:,:,:) :: sw_ref_dif
48	real(kind=jprb), allocatable, dimension(:,:,:) :: sw_tra_dif
49	! Shortwave reflectance and transmittance to direct radiation,
50	! where the transmittance has separate direct-to-diffuse
51	! transmittance (including scattering) and direct-to-direct
52	! transmittance (without scattering)
53	real(kind=jprb), allocatable, dimension(:,:,:) :: sw_ref_dir
54	real(kind=jprb), allocatable, dimension(:,:,:) :: sw_tra_dir_dif
55	real(kind=jprb), allocatable, dimension(:,:,:) :: sw_tra_dir_dir
56	! Fraction of direct radiation at canyon top that is absorbed by
57	! the wall and the atmosphere
58	real(kind=jprb), allocatable, dimension(:,:,:) :: sw_wall_abs_dir
59	real(kind=jprb), allocatable, dimension(:,:,:) :: sw_air_abs_dir
60	! Ratio of diffuse absorption in a street canyon by the wall
61	! (rather than the air or vegetation in the canyon)
62	real(kind=jprb), allocatable, dimension(:,:,:) :: sw_wall_abs_frac_dif
63	real(kind=jprb), allocatable, dimension(:,:,:) :: lw_wall_abs_frac
64	! Shortwave direct and diffuse albedo at the top of a region
65	real(kind=jprb), allocatable, dimension(:,:,:) :: sw_albedo_direct_reg
66	real(kind=jprb), allocatable, dimension(:,:,:) :: sw_albedo_diffuse_reg
67	! Longwave reflectance, transmittance and source (note that source
68	! is the same up and down because the canopy temperature is
69	! assumed constant with height)
70	real(kind=jprb), allocatable, dimension(:,:,:) :: lw_reflectance
71	real(kind=jprb), allocatable, dimension(:,:,:) :: lw_transmittance
72	real(kind=jprb), allocatable, dimension(:,:,:) :: lw_source
73	! Longwave emission and emissivity passed to the atmosphere scheme
74	! from the top of a region
75	real(kind=jprb), allocatable, dimension(:,:,:) :: lw_emissivity_region
76	real(kind=jprb), allocatable, dimension(:,:,:) :: lw_emission_region
77	! Total emission from wall and the air in an urban canopy, needed
78	! at the final partitioning stage to determine net fluxes into
79	! wall and air
80	real(kind=jprb), allocatable, dimension(:,:,:) :: lw_total_wall_emission
81	real(kind=jprb), allocatable, dimension(:,:,:) :: lw_total_canopy_emission
82
83	! Number of bands in which calculations are to be performed (can
84	! match either the spectral resolution of the atmosphere or that
85	! of the surface data)
86	integer :: nswbands, nlwbands
87
88	! Do we represent gas radiative transfer in street/vegetation
89	! canopies?
90	!logical :: do_canopy_gases_sw = .false.
91	!logical :: do_canopy_gases_lw = .false.
92
93	! Do we use the full spectrum?
94	!logical :: use_full_spectrum_sw = .false.
95	!logical :: use_full_spectrum_lw = .true.
96
97	contains
98
99	procedure :: allocate => allocate_surface_intermediate
100	procedure :: deallocate => deallocate_surface_intermediate
101	procedure :: calc_boundary_conditions_sw
102	procedure :: calc_boundary_conditions_lw
103	procedure :: calc_boundary_conditions
104	procedure :: partition_fluxes
105
106	end type surface_intermediate_type
107
108	contains
109
110	!---------------------------------------------------------------------
111	subroutine allocate_surface_intermediate(this, istartcol, iendcol, &
112	& config, surface)
113
114	use radiation_config, only : config_type
115	use radsurf_properties, only : surface_type
116
117	class(surface_intermediate_type), intent(inout) :: this
118	integer(kind=jpim), intent(in) :: istartcol, iendcol
119	type(config_type), intent(in) :: config
120	type(surface_type), intent(in) :: surface
121
122	call this%deallocate
123
124	!this%use_full_spectrum_sw = config%use_canopy_full_spectrum_sw
125	!this%use_full_spectrum_lw = config%use_canopy_full_spectrum_lw
126	! Assume that canopy gases will only be used if we do surface
127	! calculations at the full spectral resolution
128	!this%do_canopy_gases_sw = config%do_canopy_gases_sw
129	!this%do_canopy_gases_lw = config%do_canopy_gases_lw
130
131	if (config%use_canopy_full_spectrum_sw) then
132	! Calculations will be performed at the same spectral resolution
133	! as in the atmosphere
134	this%nswbands = config%n_g_sw
135	else
136	! Calculations will be performed at the same spectral resolution
137	! as the input surface-property data
138	this%nswbands = surface%nalbedobands
139	end if
140
141	if (config%use_canopy_full_spectrum_lw) then
142	! Calculations will be performed at the same spectral resolution
143	! as in the atmosphere
144	this%nlwbands = config%n_g_lw
145	else
146	! Calculations will be performed at the same spectral resolution
147	! as the input surface-property data
148	this%nlwbands = surface%nemissbands
149	end if
150
151	if (config%do_lw) then
152	allocate(this%planck_facet (this%nlwbands,surface%nfacet,istartcol:iendcol))
153	allocate(this%lw_emissivity (this%nlwbands,surface%nfacet,istartcol:iendcol))
154	end if
155	if (config%do_sw) then
156	allocate(this%sw_albedo_direct (this%nswbands,surface%nfacet,istartcol:iendcol))
157	allocate(this%sw_albedo_diffuse(this%nswbands,surface%nfacet,istartcol:iendcol))
158	end if
159
160	if (surface%nregion > 0) then
161	if (config%do_lw) then
162	allocate(this%planck_region (this%nlwbands,surface%nregion,istartcol:iendcol))
163	allocate(this%lw_reflectance (this%nlwbands,surface%nregion,istartcol:iendcol))
164	allocate(this%lw_transmittance(this%nlwbands,surface%nregion,istartcol:iendcol))
165	allocate(this%lw_source (this%nlwbands,surface%nregion,istartcol:iendcol))
166	allocate(this%lw_emission_region(this%nlwbands,surface%nregion,istartcol:iendcol))
167	allocate(this%lw_emissivity_region(this%nlwbands,surface%nregion,istartcol:iendcol))
168	allocate(this%lw_wall_abs_frac(this%nlwbands,surface%nregion,istartcol:iendcol))
169	allocate(this%lw_total_wall_emission(this%nlwbands,surface%nregion,istartcol:iendcol))
170	allocate(this%lw_total_canopy_emission(this%nlwbands,surface%nregion,istartcol:iendcol))
171	end if
172	if (config%do_sw) then
173	allocate(this%sw_ref_dif (this%nswbands,surface%nregion,istartcol:iendcol))
174	allocate(this%sw_tra_dif (this%nswbands,surface%nregion,istartcol:iendcol))
175	allocate(this%sw_ref_dir (this%nswbands,surface%nregion,istartcol:iendcol))
176	allocate(this%sw_tra_dir_dif(this%nswbands,surface%nregion,istartcol:iendcol))
177	allocate(this%sw_tra_dir_dir(this%nswbands,surface%nregion,istartcol:iendcol))
178	allocate(this%sw_wall_abs_dir(this%nswbands,surface%nregion,istartcol:iendcol))
179	allocate(this%sw_air_abs_dir(this%nswbands,surface%nregion,istartcol:iendcol))
180	allocate(this%sw_wall_abs_frac_dif(this%nswbands,surface%nregion,istartcol:iendcol))
181	allocate(this%sw_albedo_diffuse_reg(this%nswbands,surface%nregion,istartcol:iendcol))
182	allocate(this%sw_albedo_direct_reg(this%nswbands,surface%nregion,istartcol:iendcol))
183	end if
184	end if
185
186	end subroutine allocate_surface_intermediate
187
188	!---------------------------------------------------------------------
189	subroutine deallocate_surface_intermediate(this)
190
191	class(surface_intermediate_type), intent(inout) :: this
192
193	if (allocated(this%planck_facet)) then
194	deallocate(this%planck_facet)
195	end if
196	if (allocated(this%planck_region)) then
197	deallocate(this%planck_region)
198	end if
199	if (allocated(this%lw_emissivity)) then
200	deallocate(this%lw_emissivity)
201	end if
202	if (allocated(this%sw_albedo_direct)) then
203	deallocate(this%sw_albedo_direct)
204	end if
205	if (allocated(this%sw_albedo_diffuse)) then
206	deallocate(this%sw_albedo_diffuse)
207	end if
208
209	if (allocated(this%lw_reflectance)) then
210	deallocate(this%lw_reflectance)
211	end if
212	if (allocated(this%lw_transmittance)) then
213	deallocate(this%lw_transmittance)
214	end if
215	if (allocated(this%lw_source)) then
216	deallocate(this%lw_source)
217	end if
218	if (allocated(this%lw_emissivity_region)) then
219	deallocate(this%lw_emissivity_region)
220	end if
221	if (allocated(this%lw_emission_region)) then
222	deallocate(this%lw_emission_region)
223	end if
224	if (allocated(this%lw_wall_abs_frac)) then
225	deallocate(this%lw_wall_abs_frac)
226	end if
227	if (allocated(this%lw_total_wall_emission)) then
228	deallocate(this%lw_total_wall_emission)
229	end if
230	if (allocated(this%lw_total_canopy_emission)) then
231	deallocate(this%lw_total_canopy_emission)
232	end if
233
234	if (allocated(this%sw_ref_dif)) then
235	deallocate(this%sw_ref_dif)
236	end if
237	if (allocated(this%sw_tra_dif)) then
238	deallocate(this%sw_tra_dif)
239	end if
240	if (allocated(this%sw_ref_dir)) then
241	deallocate(this%sw_ref_dir)
242	end if
243	if (allocated(this%sw_tra_dir_dif)) then
244	deallocate(this%sw_tra_dir_dif)
245	end if
246	if (allocated(this%sw_tra_dir_dir)) then
247	deallocate(this%sw_tra_dir_dir)
248	end if
249	if (allocated(this%sw_wall_abs_frac_dif)) then
250	deallocate(this%sw_wall_abs_frac_dif)
251	end if
252	if (allocated(this%sw_wall_abs_dir)) then
253	deallocate(this%sw_wall_abs_dir)
254	end if
255	if (allocated(this%sw_air_abs_dir)) then
256	deallocate(this%sw_air_abs_dir)
257	end if
258	if (allocated(this%sw_albedo_direct_reg)) then
259	deallocate(this%sw_albedo_direct_reg)
260	end if
261	if (allocated(this%sw_albedo_diffuse_reg)) then
262	deallocate(this%sw_albedo_diffuse_reg)
263	end if
264
265	end subroutine deallocate_surface_intermediate
266
267
268
269	!---------------------------------------------------------------------
270	! Use the detailed physical properties of the surface in "surface",
271	! and optionally the atmospheric optical properties of the lowest
272	! atmospheric level, to work out the shortwave direct/diffuse albedo
273	! that is presented to the rest of the radiation scheme, and stored
274	! in "single_level". Also, store the necessary information so that
275	! after the atmospheric radiation scheme has been run, the net
276	! fluxes at each.
277	subroutine calc_boundary_conditions_sw(this, istartcol, iendcol, & ! in
278	& config, surface, & ! in
279	& single_level, & ! out
280	& ext_sw_air, ssa_sw_air, g_sw_air) ! in, optional
281
282	use yomhook, only : lhook, dr_hook
283
284	use radiation_io, only : nulerr, radiation_abort
285	use radiation_config, only : config_type
286	use radiation_single_level, only : single_level_type
287	use radsurf_properties, only : surface_type, ITileFlat,ITileVegetation,ITileUrban3D
288	use radiation_two_stream, only : calc_two_stream_gammas_sw, &
289	& calc_reflectance_transmittance_sw, &
290	& calc_reflectance_transmittance_z_sw
291
292	use radiation_constants, only: Pi, GasConstantDryAir, AccelDueToGravity
293
294	class(surface_intermediate_type), intent(inout) :: this
295	integer(kind=jpim), intent(in) :: istartcol, iendcol
296	type(config_type), intent(in) :: config
297	type(surface_type), intent(in) :: surface
298	type(single_level_type), intent(inout) :: single_level
299
300	! Input properties of the air in the lowest model level:
301	! extinction coefficient (m-1), single-scattering albedo and
302	! asymmetry factor
303	real(kind=jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol), optional &
304	& :: ext_sw_air, ssa_sw_air, g_sw_air
305
306	! Shortwave region properties
307	real(kind=jprb), dimension(this%nswbands) :: od_sw_region, ssa_sw_region, g_sw_region
308	! Extinction coefficient better for urban areas
309	real(kind=jprb), dimension(this%nswbands) :: ext_sw_region
310	! Optical depth of air
311	real(kind=jprb), dimension(this%nswbands) :: od_sw_air
312
313	real(kind=jprb) :: tile_fraction, cos_sza
314
315	! Exchange coefficients (m-1) for direct and diffuse radiation
316	! into building walls
317	real(kind=jprb) :: fdiff, fdir
318
319	! One minus the building fraction
320	real(kind=jprb) :: canyon_fraction
321
322	! Tangent of solar zenith angle
323	real(kind=jprb) :: tan_sza
324
325	! Two-stream coefficients
326	real(jprb), dimension(this%nswbands) :: gamma1_sw, gamma2_sw
327	real(jprb), dimension(this%nswbands) :: gamma0_sw, gamma3_sw, gamma4_sw
328
329	real(jprb), dimension(this%nswbands) :: inv_denominator_sw
330
331	integer(kind=jpim) :: jcol,jtile,iregion
332
333	! Mapping from albedo bands to reordered shortwave g points
334	integer(kind=jpim) :: i_albedo_from_g(config%n_g_sw)
335
336	! Indices to different facets
337	integer(kind=jpim) :: iground, iwall, iroof
338
339	real(jprb) :: hook_handle
340
341	if (lhook) call dr_hook('radsurf_intermediate:calc_boundary_conditions_sw',0,hook_handle)
342
343	if (present(ext_sw_air)) then
344	if (.not. present(ssa_sw_air) .or. .not. present(g_sw_air)) then
345	write(nulerr,'(a)') '*** Error: if ext_sw_air is provided then ssa_sw_air and g_sw_air must also be provided'
346	call radiation_abort()
347	end if
348	end if
349
350	if (size(single_level%sw_albedo,2) /= this%nswbands) then
351	write(nulerr,'(a,i0,a,i0)') '*** Error: single-level albedo has ', &
352	& size(single_level%sw_albedo,2), ' bands, needs ', this%nswbands
353	call radiation_abort()
354	end if
355
356	if (config%use_canopy_full_spectrum_sw) then
357	! Put shortwave albedo on g-point grid and permute
358	i_albedo_from_g = config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw)
359	do jcol = istartcol,iendcol
360	this%sw_albedo_diffuse(:,:,jcol) = surface%sw_albedo(jcol,i_albedo_from_g,:)
361	if (allocated(surface%sw_albedo_direct)) then
362	this%sw_albedo_direct(:,:,jcol) = surface%sw_albedo_direct(jcol,i_albedo_from_g,:)
363	else
364	this%sw_albedo_direct = this%sw_albedo_diffuse
365	end if
366	end do
367	else
368	do jcol = istartcol,iendcol
369	! No change to spectral resolution: simply permute
370	this%sw_albedo_diffuse(:,:,jcol) = surface%sw_albedo(jcol,:,:)
371	if (allocated(surface%sw_albedo_direct)) then
372	this%sw_albedo_direct(:,:,jcol) = surface%sw_albedo_direct(jcol,:,:)
373	else
374	this%sw_albedo_direct = this%sw_albedo_diffuse
375	end if
376	end do
377	end if
378
379	if (surface%is_simple) then
380	! We have a "traditional" representation: one flat tile
381	single_level%sw_albedo_direct(istartcol:iendcol,:) &
382	& = transpose(this%sw_albedo_direct(:,1,istartcol:iendcol))
383	single_level%sw_albedo (istartcol:iendcol,:) &
384	& = transpose(this%sw_albedo_diffuse(:,1,istartcol:iendcol))
385	else
386	! More complex description of surface
387
388	! Firstly initialize outputs to zero
389	single_level%sw_albedo_direct(istartcol:iendcol,:) = 0.0_jprb
390	single_level%sw_albedo (istartcol:iendcol,:) = 0.0_jprb
391
392	! Loop over column and tile
393	do jcol = istartcol,iendcol
394	cos_sza = single_level%cos_sza(jcol)
395
396	do jtile = 1,surface%ntile
397	tile_fraction = surface%tile_fraction(jcol,jtile)
398
399	if (tile_fraction > 0.0_jprb) then
400
401	select case (surface%i_representation(jtile))
402	case (ITileFlat)
403	! SIMPLE FLAT TILE
404
405	! Add the contribution from this simple flat tile to the
406	! accumulated values for the column
407	iground = surface%i_ground_facet(jcol,jtile)
408	single_level%sw_albedo_direct(jcol,:) = single_level%sw_albedo_direct(jcol,:) &
409	& + tile_fraction * this%sw_albedo_direct(:,iground,jcol)
410	single_level%sw_albedo(jcol,:) = single_level%sw_albedo(jcol,:) &
411	& + tile_fraction * this%sw_albedo_diffuse(:,iground,jcol)
412
413	case (ITileVegetation)
414	! HORIZONTALLY HOMOGENEOUS VEGETATION CANOPY WITH
415	! SELLERS-LIKE FORMULATION
416
417	iground = surface%i_ground_facet(jcol,jtile)
418	iregion = surface%i_region_1(jcol,jtile)
419
420	! Shortwave calculation
421	if (present(ext_sw_air)) then
422	od_sw_air = surface%canopy_depth(jcol,jtile)*ext_sw_air(:,jcol)
423	od_sw_region = od_sw_air + surface%vegetation_optical_depth(jcol,jtile)
424	ssa_sw_region= (ssa_sw_air(:,jcol)*od_sw_air &
425	& + surface%vegetation_optical_depth(jcol,jtile) &
426	& * surface%vegetation_sw_albedo(jcol,:,jtile)) &
427	& / od_sw_region
428	! Assume asymmetry factor of vegetation is zero
429	g_sw_region = g_sw_air(:,jcol)ssa_sw_air(:,jcol)od_sw_air &
430	& / (ssa_sw_region*od_sw_region)
431	else
432	od_sw_region = surface%vegetation_optical_depth(jcol,jtile)
433	ssa_sw_region= surface%vegetation_sw_albedo(jcol,:,jtile)
434	g_sw_region = 0.0_jprb
435	end if
436	call calc_two_stream_gammas_sw(this%nswbands, &
437	& cos_sza, ssa_sw_region, g_sw_region, &
438	& gamma1_sw, gamma2_sw, gamma3_sw)
439	call calc_reflectance_transmittance_sw(this%nswbands, &
440	& cos_sza, od_sw_region, ssa_sw_region, &
441	& gamma1_sw, gamma2_sw, gamma3_sw, &
442	& this%sw_ref_dif(:,iregion,jcol), this%sw_tra_dif(:,iregion,jcol), &
443	& this%sw_ref_dir(:,iregion,jcol), this%sw_tra_dir_dif(:,iregion,jcol), &
444	& this%sw_tra_dir_dir(:,iregion,jcol))
445
446	! Shortwave adding method for a single layer
447	inv_denominator_sw(:) = 1.0_jprb &
448	& / (1.0_jprb - this%sw_albedo_diffuse(:,iground,jcol) &
449	& *this%sw_ref_dif(:,iregion,jcol))
450	this%sw_albedo_diffuse_reg(:,iregion,jcol) = this%sw_ref_dif(:,iregion,jcol) &
451	& + this%sw_tra_dif(:,iregion,jcol)**2 &
452	& * this%sw_albedo_diffuse(:,iground,jcol) * inv_denominator_sw(:)
453	this%sw_albedo_direct_reg(:,iregion,jcol) = this%sw_ref_dir(:,iregion,jcol) &
454	& + (this%sw_tra_dir_dir(:,iregion,jcol)*this%sw_albedo_direct (:,iground,jcol) &
455	& +this%sw_tra_dir_dif(:,iregion,jcol)*this%sw_albedo_diffuse(:,iground,jcol)) &
456	& * this%sw_tra_dif(:,iregion,jcol) * inv_denominator_sw(:)
457	single_level%sw_albedo(jcol,:) = single_level%sw_albedo(jcol,:) &
458	& + tile_fraction * this%sw_albedo_diffuse_reg(:,iregion,jcol)
459	single_level%sw_albedo_direct(jcol,:) = single_level%sw_albedo_direct(jcol,:) &
460	& + tile_fraction * this%sw_albedo_direct_reg(:,iregion,jcol)
461
462	case (ITileUrban3D)
463	! URBAN CANOPY WITH NO VEGETATION, TREATED USING
464	! SPARTACUS METHODOLOGY
465
466	iground = surface%i_ground_facet(jcol,jtile)
467	iwall = surface%i_wall_facet(jcol,jtile)
468	iroof = surface%i_roof_facet(jcol,jtile)
469	iregion = surface%i_region_1(jcol,jtile)
470
471	canyon_fraction = 1.0_jprb - surface%building_fraction(jcol,jtile)
472
473	fdiff = 0.5_jprb * surface%building_normalized_perimeter(jcol,jtile) / canyon_fraction
474
475	tan_sza = sqrt(1.0_jprb / (cos_sza*cos_sza) - 1.0_jprb)
476	fdir = surface%building_normalized_perimeter(jcol,jtile) * tan_sza / (Pi * canyon_fraction)
477
478	if (present(ext_sw_air)) then
479	ext_sw_region = ext_sw_air(:,jcol)
480	ssa_sw_region = ssa_sw_air(:,jcol)
481	g_sw_region = g_sw_air(:,jcol)
482	else
483	ext_sw_region = 0.0_jprb
484	ssa_sw_region = 0.0_jprb
485	g_sw_region = 0.0_jprb
486	end if
487
488	! Get gammas for atmosphere only
489	call calc_two_stream_gammas_sw(this%nswbands, &
490	& cos_sza, ssa_sw_region, g_sw_region, &
491	& gamma1_sw, gamma2_sw, gamma3_sw)
492
493	! At this point gamma1_sw-gamma2_sw is the rate of
494	! absorption per unit optical depth of the air in the
495	! canyon
496	this%sw_wall_abs_frac_dif(:,iregion,jcol) &
497	& = fdiff * (1.0_jprb - this%sw_albedo_diffuse(:,iwall,jcol))
498	this%sw_wall_abs_frac_dif(:,iregion,jcol) = this%sw_wall_abs_frac_dif(:,iregion,jcol) &
499	& / max(1.0e-8_jprb,ext_sw_region*(gamma1_sw-gamma2_sw) &
500	& + this%sw_wall_abs_frac_dif(:,iregion,jcol))
501
502	gamma4_sw = 1.0_jprb - gamma3_sw
503	gamma0_sw = ext_sw_region / cos_sza + fdir
504	gamma1_sw = ext_sw_region * gamma1_sw &
505	& + fdiff * (1.0_jprb - 0.5_jprb*this%sw_albedo_diffuse(:,iwall,jcol))
506	gamma2_sw = ext_sw_region * gamma2_sw &
507	& + fdiff * 0.5_jprb*this%sw_albedo_diffuse(:,iwall,jcol)
508	gamma3_sw = ext_sw_region * ssa_sw_region * gamma3_sw &
509	& + 0.5_jprb * fdir * this%sw_albedo_direct(:,iwall,jcol)
510	gamma4_sw = ext_sw_region * ssa_sw_region * gamma4_sw &
511	& + 0.5_jprb * fdir * this%sw_albedo_direct(:,iwall,jcol)
512
513	call calc_reflectance_transmittance_z_sw(this%nswbands, &
514	& cos_sza, surface%canopy_depth(jcol,jtile), gamma0_sw, &
515	& gamma1_sw, gamma2_sw, gamma3_sw, gamma3_sw, &
516	& this%sw_ref_dif(:,iregion,jcol), this%sw_tra_dif(:,iregion,jcol), &
517	& this%sw_ref_dir(:,iregion,jcol), this%sw_tra_dir_dif(:,iregion,jcol), &
518	& this%sw_tra_dir_dir(:,iregion,jcol))
519
520	! Compute fraction of direct at canyon top that is absorbed by wall and air
521	this%sw_wall_abs_dir(:,iregion,jcol) &
522	& = (1.0_jprb - this%sw_tra_dir_dir(:,iregion,jcol)) &
523	& * fdir * (1.0_jprb - this%sw_albedo_direct(:,iwall,jcol)) * cos_sza &
524	& / max(1.0e-8_jprb, fdir*cos_sza + ext_sw_region)
525	this%sw_air_abs_dir(:,iregion,jcol)&
526	& = (1.0_jprb - this%sw_tra_dir_dir(:,iregion,jcol)) &
527	& * ext_sw_region * (1.0_jprb - ssa_sw_region) &
528	& / max(1.0e-8_jprb, fdir*cos_sza + ext_sw_region)
529
530	! Add roof component
531	single_level%sw_albedo(jcol,:) = single_level%sw_albedo(jcol,:) &
532	& + tile_fraction * surface%building_fraction(jcol,jtile) &
533	& * this%sw_albedo_diffuse(:,iroof,jcol)
534	single_level%sw_albedo_direct(jcol,:) = single_level%sw_albedo_direct(jcol,:) &
535	& + tile_fraction * surface%building_fraction(jcol,jtile) &
536	& * this%sw_albedo_direct(:,iroof,jcol)
537
538	! Add canyon component using shortwave adding method for a single layer
539	inv_denominator_sw(:) = 1.0_jprb &
540	& / (1.0_jprb - this%sw_albedo_diffuse(:,iground,jcol) &
541	& *this%sw_ref_dif(:,iregion,jcol))
542	this%sw_albedo_diffuse_reg(:,iregion,jcol) = this%sw_ref_dif(:,iregion,jcol) &
543	& + this%sw_tra_dif(:,iregion,jcol)**2 &
544	& * this%sw_albedo_diffuse(:,iground,jcol) * inv_denominator_sw(:)
545	this%sw_albedo_direct_reg(:,iregion,jcol) = this%sw_ref_dir(:,iregion,jcol) &
546	& + (this%sw_tra_dir_dir(:,iregion,jcol)*this%sw_albedo_direct (:,iground,jcol) &
547	& +this%sw_tra_dir_dif(:,iregion,jcol)*this%sw_albedo_diffuse(:,iground,jcol)) &
548	& * this%sw_tra_dif(:,iregion,jcol) * inv_denominator_sw(:)
549	single_level%sw_albedo(jcol,:) = single_level%sw_albedo(jcol,:) &
550	& + tile_fraction * canyon_fraction * this%sw_albedo_diffuse_reg(:,iregion,jcol)
551	single_level%sw_albedo_direct(jcol,:) = single_level%sw_albedo_direct(jcol,:) &
552	& + tile_fraction * canyon_fraction * this%sw_albedo_direct_reg(:,iregion,jcol)
553
554	end select
555	end if
556	end do
557	end do
558	end if
559
560	if (lhook) call dr_hook('radsurf_intermediate:calc_boundary_conditions_sw',1,hook_handle)
561
562	end subroutine calc_boundary_conditions_sw
563
564
565	!---------------------------------------------------------------------
566	! As calc_boundary_conditions_sw, but for the longwave.
567	subroutine calc_boundary_conditions_lw(this, istartcol, iendcol, & ! in
568	& config, surface, & ! in
569	& single_level, & ! out
570	& ext_lw_air, ssa_lw_air, g_lw_air) ! in, optional
571
572	use yomhook, only : lhook, dr_hook
573
574	use radiation_io, only : nulerr, radiation_abort
575	use radiation_config, only : config_type
576	use radiation_single_level, only : single_level_type
577	use radsurf_properties, only : surface_type, ITileFlat,ITileVegetation,ITileUrban3D
578	use radiation_two_stream, only : calc_two_stream_gammas_lw, &
579	& calc_reflectance_transmittance_isothermal_lw, &
580	& LwDiffusivity
581	use radiation_ifs_rrtm, only : planck_function
582	use radiation_constants, only: Pi, GasConstantDryAir, AccelDueToGravity, StefanBoltzmann
583
584	class(surface_intermediate_type), intent(inout) :: this
585	integer(kind=jpim), intent(in) :: istartcol, iendcol
586	type(config_type), intent(in) :: config
587	type(surface_type), intent(in) :: surface
588	type(single_level_type), intent(inout) :: single_level
589
590	! Input properties of the air in the lowest model level:
591	! extinction coefficient (m-1), single-scattering albedo and
592	! asymmetry factor
593	real(kind=jprb), intent(in), dimension(config%n_g_lw,istartcol:iendcol), optional &
594	& :: ext_lw_air, ssa_lw_air, g_lw_air
595
596	! Longwave region properties
597	real(kind=jprb), dimension(this%nlwbands) :: od_lw_region, ssa_lw_region, g_lw_region
598
599	! Optical depth of air
600	real(kind=jprb), dimension(this%nlwbands) :: od_lw_air
601
602	! Effective Planck function of urban canopy as a weighted average
603	! of wall and air emission
604	real(kind=jprb), dimension(this%nlwbands) :: planck_canopy
605
606	! Vegetation emissivity using local spectral representation
607	real(kind=jprb), dimension(this%nlwbands) :: vegetation_lw_emissivity
608
609	real(kind=jprb) :: tile_fraction
610
611	! Exchange coefficient (m-1) for diffuse radiation into building
612	! walls, multiplied by canyon depth (m) to get a dimensionless
613	! analogue for optical depth
614	real(kind=jprb) :: od_lw_wall
615
616	! One minus the building fraction
617	real(kind=jprb) :: canyon_fraction
618
619	! Two-stream coefficients
620	real(jprb), dimension(this%nlwbands) :: gamma1_lw, gamma2_lw
621
622	real(jprb), dimension(this%nlwbands) :: inv_denominator_lw
623
624	integer(kind=jpim) :: jcol,jtile,jfacet,iregion
625
626	! Indices to different facets
627	integer(kind=jpim) :: iground, iwall, iroof
628
629	real(jprb) :: hook_handle
630
631	if (lhook) call dr_hook('radsurf_intermediate:calc_boundary_conditions_lw',0,hook_handle)
632
633	if (present(ssa_lw_air)) then
634	if (.not. present(g_lw_air)) then
635	write(nulerr,'(a)') '*** Error: if ssa_lw_air is provided then g_lw_air must also be provided'
636	call radiation_abort()
637	end if
638	end if
639
640	if (size(single_level%lw_emissivity,2) /= this%nlwbands) then
641	write(nulerr,'(a,i0,a,i0)') '*** Error: single-level emissivity has ', &
642	& size(single_level%lw_emissivity,2), ' bands, needs ', this%nlwbands
643	call radiation_abort()
644	end if
645
646	! FIX: This section ought to check what tiles are present in each column
647	if (config%use_canopy_full_spectrum_lw) then
648	! Put longwave emissivity on g-point grid and permute
649	do jcol = istartcol,iendcol
650	this%lw_emissivity(:,:,jcol) = surface%lw_emissivity(jcol,config%i_emiss_from_band_lw( &
651	& config%i_band_from_reordered_g_lw),:)
652	end do
653	! Compute Planck function at each facet
654	do jfacet = 1,surface%nfacet
655	do jcol = istartcol,iendcol
656	! FIX this function assumes contiguous data for planck_facet
657	! so copies into a temporary
658	call planck_function(config, &
659	& surface%skin_temperature(jcol,jfacet), &
660	& this%planck_facet(:,jfacet,jcol))
661	end do
662	end do
663	do jtile = 1,surface%ntile
664	do jcol = istartcol,iendcol
665	iregion = surface%i_region_1(jcol,jtile)
666	if (iregion > 0) then
667	call planck_function(config, &
668	& surface%canopy_temperature(jcol,jtile), &
669	& this%planck_region(:,iregion,jcol:jcol))
670	end if
671	end do
672	end do
673	else
674	if (surface%nemissbands /= 1) then
675	write(nulerr,'(a)') '*** Error: insufficient information to compute Planck function in emissivity bands'
676	call radiation_abort()
677	end if
678	do jcol = istartcol,iendcol
679	! No change to spectral resolution: simply permute
680	this%lw_emissivity(:,:,jcol) = surface%lw_emissivity(jcol,:,:)
681	end do
682	! Broadband calculation: Stefan-Boltzmann law
683	do jfacet = 1,surface%nfacet
684	this%planck_facet(1,jfacet,istartcol:iendcol) &
685	& = StefanBoltzmannsurface%skin_temperature(istartcol:iendcol,jfacet)*4
686	end do
687	do jtile = 1,surface%ntile
688	do jcol = istartcol,iendcol
689	iregion = surface%i_region_1(jcol,jtile)
690	if (iregion > 0) then
691	this%planck_region(1,iregion,jcol) &
692	& = StefanBoltzmann*surface%canopy_temperature(jcol,jtile)
693	end if
694	end do
695	end do
696	end if
697
698	if (surface%is_simple) then
699	! We have a "traditional" representation: one flat tile
700	single_level%lw_emissivity(istartcol:iendcol,:) &
701	& = transpose(this%lw_emissivity(:,1,istartcol:iendcol))
702	single_level%lw_emission(istartcol:iendcol,:) &
703	& = transpose(this%planck_facet(:,1,istartcol:iendcol) &
704	& * this%lw_emissivity(:,1,istartcol:iendcol))
705
706	else
707	! More complex description of surface
708
709	! Firstly initialize outputs to zero
710	single_level%lw_emissivity(istartcol:iendcol,:) = 0.0_jprb
711	single_level%lw_emission (istartcol:iendcol,:) = 0.0_jprb
712
713	! Loop over column and tile
714	do jcol = istartcol,iendcol
715
716	do jtile = 1,surface%ntile
717	tile_fraction = surface%tile_fraction(jcol,jtile)
718
719	if (tile_fraction > 0.0_jprb) then
720
721	select case (surface%i_representation(jtile))
722	case (ITileFlat)
723	! SIMPLE FLAT TILE
724
725	! Add the contribution from this simple flat tile to the
726	! accumulated values for the column
727	iground = surface%i_ground_facet(jcol,jtile)
728	single_level%lw_emissivity(jcol,:) = single_level%lw_emissivity(jcol,:) &
729	& + tile_fraction*this%lw_emissivity(:,iground,jcol)
730	single_level%lw_emission(jcol,:) = single_level%lw_emission(jcol,:) + tile_fraction &
731	& * (this%planck_facet(:,iground,jcol)*this%lw_emissivity(:,iground,jcol))
732
733	case (ITileVegetation)
734	! HORIZONTALLY HOMOGENEOUS VEGETATION CANOPY
735
736	iground = surface%i_ground_facet(jcol,jtile)
737	iregion = surface%i_region_1(jcol,jtile)
738
739	! Convert vegetation emissivity to local spectral representation
740	if (config%use_canopy_full_spectrum_lw) then
741	vegetation_lw_emissivity = surface%vegetation_lw_emissivity(jcol, &
742	& config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), &
743	& jtile)
744	else
745	vegetation_lw_emissivity = surface%vegetation_lw_emissivity(jcol,:,jtile)
746	end if
747
748	if (present(ext_lw_air)) then
749	od_lw_air = surface%canopy_depth(jcol,jtile)*ext_lw_air(:,jcol)
750	od_lw_region = od_lw_air + surface%vegetation_optical_depth(jcol,jtile)
751	if (present(ssa_lw_air)) then
752	ssa_lw_region= (ssa_lw_air(:,jcol)*od_lw_air &
753	& + surface%vegetation_optical_depth(jcol,jtile) &
754	& * (1.0_jprb-vegetation_lw_emissivity)) &
755	& / od_lw_region
756	! Assume asymmetry factor of vegetation is zero
757	g_lw_region = g_lw_air(:,jcol)ssa_lw_air(:,jcol)od_lw_air &
758	& / (ssa_lw_region * od_lw_region)
759	else
760	! No longwave scattering properties for air
761	ssa_lw_region= surface%vegetation_optical_depth(jcol,jtile) &
762	& * (1.0_jprb-vegetation_lw_emissivity) &
763	& / od_lw_region
764	! Assume asymmetry factor of vegetation is zero
765	g_lw_region = 0.0_jprb
766
767	end if
768	else
769	! No gases
770	od_lw_region = surface%vegetation_optical_depth(jcol,jtile)
771	ssa_lw_region= 1.0_jprb-vegetation_lw_emissivity
772	g_lw_region = 0.0_jprb
773	end if
774
775	call calc_two_stream_gammas_lw(this%nlwbands, &
776	& ssa_lw_region, g_lw_region, &
777	& gamma1_lw, gamma2_lw)
778	call calc_reflectance_transmittance_isothermal_lw(this%nlwbands, &
779	& od_lw_region, gamma1_lw, gamma2_lw, this%planck_region(:,iregion,jcol), &
780	& this%lw_reflectance(:,iregion,jcol), this%lw_transmittance(:,iregion,jcol), &
781	& this%lw_source(:,iregion,jcol))
782
783	! Longwave adding method for a single layer
784	inv_denominator_lw(:) = 1.0_jprb &
785	& / (1.0_jprb - (1.0_jprb - this%lw_emissivity(:,iground,jcol)) &
786	& *this%lw_reflectance(:,iregion,jcol))
787	single_level%lw_emissivity(jcol,:) = single_level%lw_emissivity(jcol,:) &
788	& + tile_fraction * (1.0_jprb - (this%lw_reflectance(:,iregion,jcol) &
789	& + this%lw_transmittance(:,iregion,jcol)**2 &
790	& * (1.0_jprb - this%lw_emissivity(:,iground,jcol)) * inv_denominator_lw(:)))
791	single_level%lw_emission(jcol,:) = single_level%lw_emission(jcol,:) + tile_fraction &
792	& * (this%lw_source(:,iregion,jcol) * (1.0_jprb + inv_denominator_lw(:) &
793	& * (1.0_jprb - this%lw_emissivity(:,iground,jcol)) &
794	& * this%lw_transmittance(:,iregion,jcol)) &
795	& +this%planck_facet(:,iground,jcol)*this%lw_emissivity(:,iground,jcol) &
796	& this%lw_transmittance(:,iregion,jcol)inv_denominator_lw(:))
797
798	case (ITileUrban3D)
799	! URBAN CANOPY WITH NO VEGETATION, TREATED USING
800	! SPARTACUS METHODOLOGY
801
802	iground = surface%i_ground_facet(jcol,jtile)
803	iwall = surface%i_wall_facet(jcol,jtile)
804	iroof = surface%i_roof_facet(jcol,jtile)
805	iregion = surface%i_region_1(jcol,jtile)
806
807	canyon_fraction = 1.0_jprb - surface%building_fraction(jcol,jtile)
808
809	! fdiff woould be 0.5 * building_normalized_perimeter /
810	! canyon_fraction, but to get equivalent zenith optical
811	! depth we multiply by the building height and divide by
812	! the longwave diffusivity angle
813	! od_lw_wall = 0.5_jprb * surface%building_normalized_perimeter(jcol,jtile) &
814	! & * surface%canopy_depth(jcol,jtile) / (LwDiffusivity * canyon_fraction)
815
816	! Or first compute H/W
817	od_lw_wall = 0.5_jprb * surface%building_normalized_perimeter(jcol,jtile) &
818	& * surface%canopy_depth(jcol,jtile) / canyon_fraction
819	! And then use the Harman et al. (2004) formula for
820	! street-to-sky transmittance T=sqrt[(H/W)^2+1]-H/W, to
821	! compute equivalent optical depth knowing that the
822	! two-stream scheme treatment will be T=exp(-D*od)
823	od_lw_wall = -log(sqrt(od_lw_wall*od_lw_wall + 1) - od_lw_wall) / LwDiffusivity
824
825	if (present(ext_lw_air)) then
826	od_lw_air = ext_lw_air(:,jcol)surface%canopy_depth(jcol,jtile) !10.0_jprb
827	od_lw_region = od_lw_air + od_lw_wall
828	if (present(ssa_lw_air)) then
829	ssa_lw_region = (od_lw_air * ssa_lw_air(:,jcol) &
830	& + od_lw_wall * (1.0_jprb - this%lw_emissivity(:,iwall,jcol))) &
831	& / max(od_lw_region,1.0e-6_jprb)
832	! We assume that any scattering off the walls is
833	! equally likely to be up or down, so the effective
834	! asymmetry factor is zero
835	g_lw_region = (od_lw_air * ssa_lw_air(:,jcol)*g_lw_air(:,jcol)) &
836	& / max(od_lw_region*ssa_lw_region,1.0e-6_jprb)
837	! Effective Planck function of the canopy is the
838	! weighted average of wall and air emission
839	this%lw_total_wall_emission(:,iregion,jcol) = LwDiffusivity &
840	& * od_lw_wallthis%lw_emissivity(:,iwall,jcol)this%planck_facet(:,iwall,jcol)
841	this%lw_total_canopy_emission(:,iregion,jcol) = LwDiffusivity &
842	& * od_lw_air(1.0_jprb-ssa_lw_air(:,jcol))this%planck_region(:,iregion,jcol)
843	planck_canopy = (this%lw_total_wall_emission(:,iregion,jcol) &
844	& +this%lw_total_canopy_emission(:,iregion,jcol)) &
845	& / max(od_lw_region(1.0_jprb-ssa_lw_region)LwDiffusivity,1.0e-6_jprb)
846	else
847	ssa_lw_region = od_lw_wall * (1.0_jprb - this%lw_emissivity(:,iwall,jcol)) &
848	& / max(od_lw_region,1.0e-6_jprb)
849	! We assume that any scattering off the walls is
850	! equally likely to be up or down, so the effective
851	! asymmetry factor is zero
852	g_lw_region = 0.0_jprb
853	! Effective Planck function of the canopy is the
854	! weighted average of wall and air emission
855	this%lw_total_wall_emission(:,iregion,jcol) = LwDiffusivity &
856	& * od_lw_wallthis%lw_emissivity(:,iwall,jcol)this%planck_facet(:,iwall,jcol)
857	this%lw_total_canopy_emission(:,iregion,jcol) = LwDiffusivity &
858	& * od_lw_air*this%planck_region(:,iregion,jcol)
859	planck_canopy = (this%lw_total_wall_emission(:,iregion,jcol) &
860	& +this%lw_total_canopy_emission(:,iregion,jcol)) &
861	& / max(od_lw_region(1.0_jprb-ssa_lw_region)LwDiffusivity,1.0e-6_jprb)
862	end if
863
864	! Compute fraction of canyon absorption by the wall
865	this%lw_wall_abs_frac(:,iregion,jcol) = od_lw_wall * this%lw_emissivity(:,iwall,jcol) &
866	& / max(od_lw_region*(1.0_jprb-ssa_lw_region),1.0e-6_jprb)
867
868
869	else
870	od_lw_region = od_lw_wall
871	ssa_lw_region = 1.0_jprb - this%lw_emissivity(:,iwall,jcol);
872	g_lw_region = 0.0_jprb
873
874	! All absorption and emission is by the wall
875	this%lw_wall_abs_frac(:,iregion,jcol) = 1.0_jprb
876	this%lw_total_wall_emission(:,iregion,jcol) = LwDiffusivity &
877	& * od_lw_wallthis%lw_emissivity(:,iwall,jcol)this%planck_facet(:,iwall,jcol)
878	this%lw_total_canopy_emission(:,iregion,jcol) = 0.0_jprb
879	planck_canopy = this%planck_facet(:,iwall,jcol)
880	end if
881
882	call calc_two_stream_gammas_lw(this%nlwbands, &
883	& ssa_lw_region, g_lw_region, &
884	& gamma1_lw, gamma2_lw)
885	call calc_reflectance_transmittance_isothermal_lw(this%nlwbands, &
886	& od_lw_region, gamma1_lw, gamma2_lw, planck_canopy, &
887	& this%lw_reflectance(:,iregion,jcol), this%lw_transmittance(:,iregion,jcol), &
888	& this%lw_source(:,iregion,jcol))
889
890	! Add roof component
891	single_level%lw_emissivity(jcol,:) = single_level%lw_emissivity(jcol,:) &
892	& + tile_fraction * surface%building_fraction(jcol,jtile) &
893	& * this%lw_emissivity(:,iroof,jcol)
894	single_level%lw_emission(jcol,:) = single_level%lw_emission(jcol,:) &
895	& + tile_fraction * surface%building_fraction(jcol,jtile) &
896	& * this%lw_emissivity(:,iroof,jcol)*this%planck_facet(:,iroof,jcol)
897
898	! Add canyon component: longwave adding method for a single layer
899	inv_denominator_lw(:) = 1.0_jprb &
900	& / (1.0_jprb - (1.0_jprb - this%lw_emissivity(:,iground,jcol)) &
901	& *this%lw_reflectance(:,iregion,jcol))
902	this%lw_emissivity_region(:,iregion,jcol) &
903	& = (1.0_jprb - (this%lw_reflectance(:,iregion,jcol) &
904	& + this%lw_transmittance(:,iregion,jcol)**2 &
905	& * (1.0_jprb - this%lw_emissivity(:,iground,jcol)) * inv_denominator_lw(:)))
906	single_level%lw_emissivity(jcol,:) = single_level%lw_emissivity(jcol,:) &
907	& + tile_fraction * canyon_fraction * this%lw_emissivity_region(:,iregion,jcol)
908	this%lw_emission_region(:,iregion,jcol) &
909	& = this%lw_source(:,iregion,jcol) * (1.0_jprb + inv_denominator_lw(:) &
910	& * (1.0_jprb - this%lw_emissivity(:,iground,jcol)) &
911	& * this%lw_transmittance(:,iregion,jcol)) &
912	& +this%planck_facet(:,iground,jcol)*this%lw_emissivity(:,iground,jcol) &
913	& this%lw_transmittance(:,iregion,jcol)inv_denominator_lw(:)
914	single_level%lw_emission(jcol,:) = single_level%lw_emission(jcol,:) &
915	& + tile_fraction * canyon_fraction * this%lw_emission_region(:,iregion,jcol)
916
917	end select
918	end if
919	end do
920	end do
921	end if
922
923	if (lhook) call dr_hook('radsurf_intermediate:calc_boundary_conditions_lw',1,hook_handle)
924
925	end subroutine calc_boundary_conditions_lw
926
927
928	!---------------------------------------------------------------------
929	! Compute both shortwave and longwave boundary conditions neglecting
930	! gases within vegetation and urban canopies
931	subroutine calc_boundary_conditions_vacuum(this, istartcol, iendcol, & ! in
932	& config, surface, & ! in
933	& single_level) ! out
934
935	use radiation_config, only : config_type
936	use radiation_single_level, only : single_level_type
937	use radsurf_properties, only : surface_type
938
939	class(surface_intermediate_type), intent(inout) :: this
940	integer(kind=jpim), intent(in) :: istartcol, iendcol
941	type(config_type), intent(in) :: config
942	type(surface_type), intent(in) :: surface
943	type(single_level_type), intent(inout) :: single_level
944
945	call this%calc_boundary_conditions_sw(istartcol, iendcol, &
946	& config, surface, single_level)
947	call this%calc_boundary_conditions_lw(istartcol, iendcol, &
948	& config, surface, single_level)
949
950	end subroutine calc_boundary_conditions_vacuum
951
952
953	!---------------------------------------------------------------------
954	! Compute both shortwave and longwave boundary conditions
955	subroutine calc_boundary_conditions(this, istartcol, iendcol, & ! in
956	& config, surface, & ! in
957	& thermodynamics, gas, & ! in
958	& single_level) ! out
959
960	use radiation_config, only : config_type
961	use radiation_thermodynamics,only: thermodynamics_type
962	use radiation_gas, only : gas_type
963	use radiation_single_level, only : single_level_type
964	use radiation_ifs_rrtm, only : gas_optics
965	use radsurf_properties, only : surface_type
966	use radiation_constants, only : GasConstantDryAir, &
967	& AccelDueToGravity
968
969	class(surface_intermediate_type), intent(inout) :: this
970	integer(kind=jpim), intent(in) :: istartcol, iendcol
971	type(config_type), intent(in) :: config
972	type(surface_type), intent(in) :: surface
973	type(gas_type), intent(in) :: gas
974	type(thermodynamics_type), intent(in) :: thermodynamics
975	type(single_level_type), intent(inout) :: single_level
976
977	! Ratio of gas constant for dry air to acceleration due to gravity
978	real(jprb), parameter :: R_over_g = GasConstantDryAir / AccelDueToGravity
979
980	! Number of leves to request gas optical depths for; we only need
981	! one but gas optics scheme assumes more
982	integer, parameter :: NGasLevels = 1
983
984	! Canopy longwave optical depth
985	real(jprb), dimension(config%n_g_lw,NGasLevels,istartcol:iendcol) :: od_lw
986
987	! Canopy longwave extinction coefficient
988	real(jprb), dimension(config%n_g_lw,istartcol:iendcol) :: ext_lw
989
990	! Canopy shortwave optical depth, single scattering albedo and
991	! asymmetry factor of gases and aerosols at each shortwave g-point
992	real(jprb), dimension(config%n_g_sw,NGasLevels,istartcol:iendcol) :: od_sw, ssa_sw, g_sw
993
994	! Thickness of lowest model level
995	real(jprb) :: layer_depth
996
997	! Index to final level of full model grid
998	integer :: iendlev
999
1000	! Column loop counter, number of columns
1001	integer :: jcol, ncol
1002
1003	ncol = ubound(thermodynamics%pressure_hl,1)
1004	iendlev = ubound(thermodynamics%pressure_hl,2)-1
1005
1006	call gas_optics(ncol, NGasLevels, istartcol, iendcol, &
1007	& config, single_level, thermodynamics, gas, &
1008	& od_lw, od_sw, ssa_sw)
1009
1010	! Scale optical depths to extinction
1011	do jcol = istartcol,iendcol
1012	layer_depth = R_over_g &
1013	& * (thermodynamics%pressure_hl(jcol,iendlev+1) &
1014	& - thermodynamics%pressure_hl(jcol,iendlev)) &
1015	& * (thermodynamics%temperature_hl(jcol,iendlev) &
1016	& + thermodynamics%temperature_hl(jcol,iendlev+1)) &
1017	& / (thermodynamics%pressure_hl(jcol,iendlev) &
1018	& + thermodynamics%pressure_hl(jcol,iendlev+1))
1019	!if (config%do_sw) then
1020	! ext_sw(jcol) = ext_sw(jcol) / layer_depth
1021	!end if
1022	if (config%do_lw) then
1023	ext_lw(:,jcol) = od_lw(:,NGasLevels,jcol) / layer_depth
1024	end if
1025	end do
1026
1027	call this%calc_boundary_conditions_sw(istartcol, iendcol, &
1028	& config, surface, single_level)
1029	if (config%do_canopy_gases_lw) then
1030	call this%calc_boundary_conditions_lw(istartcol, iendcol, &
1031	& config, surface, single_level, &
1032	& ext_lw_air=ext_lw)
1033	else
1034	call this%calc_boundary_conditions_lw(istartcol, iendcol, &
1035	& config, surface, single_level)
1036	end if
1037
1038	end subroutine calc_boundary_conditions
1039
1040
1041	!---------------------------------------------------------------------
1042	subroutine partition_fluxes(this, istartcol, iendcol, config, surface, flux, &
1043	& surface_flux)
1044
1045	use yomhook, only : lhook, dr_hook
1046
1047	use radiation_flux, only : flux_type
1048	use radsurf_flux, only : surface_flux_type
1049	use radiation_config, only : config_type
1050	use radsurf_properties,only : surface_type, ITileFlat,ITileVegetation,ITileUrban3D
1051
1052	class(surface_intermediate_type), intent(in) :: this
1053	integer(kind=jpim), intent(in) :: istartcol, iendcol
1054	type(config_type), intent(in) :: config
1055	type(surface_type), intent(in) :: surface
1056	type(flux_type), intent(in) :: flux
1057	type(surface_flux_type), intent(inout) :: surface_flux
1058
1059	real(kind=jprb), dimension(config%n_g_lw) :: lw_dn_g, lw_up_g
1060	real(kind=jprb), dimension(config%n_canopy_bands_sw) &
1061	& :: sw_dn_diffuse_g, sw_dn_direct_g, sw_up_g, sw_abs_g
1062	real(kind=jprb), dimension(config%n_canopy_bands_lw) :: lw_abs_g
1063
1064	! Ratio of street planar area to wall frontal area
1065	real(kind=jprb) :: wall_scaling
1066
1067	real(kind=jprb) :: tile_fraction
1068	integer(kind=jpim) :: iground, iroof, iwall, iregion
1069	integer(kind=jpim) :: isurf ! index to lowest flux level (=nlev+1)
1070	integer(kind=jpim) :: jcol, jtile
1071
1072	real(kind=jprb) :: hook_handle
1073
1074	if (lhook) call dr_hook('radsurf_intermediate:partition_fluxes',0,hook_handle)
1075
1076	if (.not. surface%is_simple) then
1077	isurf = size(flux%lw_dn,2)
1078
1079	if (config%do_lw) then
1080	surface_flux%lw_dn_facet(istartcol:iendcol,:) = 0.0_jprb
1081	surface_flux%lw_up_facet(istartcol:iendcol,:) = 0.0_jprb
1082	surface_flux%lw_abs_canopy(istartcol:iendcol,:)=0.0_jprb
1083	end if
1084	if (config%do_sw) then
1085	surface_flux%sw_dn_facet(istartcol:iendcol,:) = 0.0_jprb
1086	surface_flux%sw_dn_direct_facet(istartcol:iendcol,:) = 0.0_jprb
1087	surface_flux%sw_up_facet(istartcol:iendcol,:) = 0.0_jprb
1088	surface_flux%sw_abs_canopy(istartcol:iendcol,:)=0.0_jprb
1089	end if
1090
1091	! Loop over column and tile
1092	do jcol = istartcol,iendcol
1093	do jtile = 1,surface%ntile
1094	tile_fraction = surface%tile_fraction(jcol,jtile)
1095
1096	if (tile_fraction > 0.0_jprb) then
1097
1098	select case (surface%i_representation(jtile))
1099	case (ITileFlat)
1100	iground = surface%i_ground_facet(jcol,jtile)
1101	if (config%do_lw) then
1102	surface_flux%lw_dn_facet(jcol,iground) = flux%lw_dn(jcol,isurf)
1103	surface_flux%lw_up_facet(jcol,iground) &
1104	& = sum(this%lw_emissivity(:,iground,jcol)*this%planck_facet(:,iground,jcol) &
1105	& + (1.0_jprb - this%lw_emissivity(:,iground,jcol)) &
1106	& * flux%lw_dn_surf_canopy(:,jcol))
1107	end if
1108	if (config%do_sw) then
1109	surface_flux%sw_dn_facet(jcol,iground) = flux%sw_dn(jcol,isurf)
1110	surface_flux%sw_dn_direct_facet(jcol,iground) = flux%sw_dn_direct(jcol,isurf)
1111	surface_flux%sw_up_facet(jcol,iground) &
1112	& = sum(this%sw_albedo_diffuse(:,iground,jcol) &
1113	& * flux%sw_dn_diffuse_surf_canopy(:,jcol) &
1114	& + this%sw_albedo_direct (:,iground,jcol) &
1115	& * flux%sw_dn_direct_surf_canopy (:,jcol))
1116	end if
1117
1118	case (ITileVegetation)
1119	iground = surface%i_ground_facet(jcol,jtile)
1120	iregion = surface%i_region_1(jcol,jtile)
1121	! Adding method in longwave and shortwave
1122	if (config%do_lw) then
1123	! Surface downwelling fluxes at each g point
1124	lw_dn_g = (this%lw_transmittance(:,iregion,jcol)*flux%lw_dn_surf_canopy(:,jcol) &
1125	& +this%lw_reflectance (:,iregion,jcol) &
1126	& this%lw_emissivity(:,iground,jcol)this%planck_facet(:,iground,jcol) &
1127	& +this%lw_source(:,iregion,jcol)) &
1128	& / (1.0_jprb - (1.0_jprb - this%lw_emissivity(:,iground,jcol)) &
1129	& *this%lw_reflectance(:,iregion,jcol))
1130	surface_flux%lw_dn_facet(jcol,iground) = sum(lw_dn_g)
1131	surface_flux%lw_up_facet(jcol,iground) &
1132	& = sum((1.0_jprb-this%lw_emissivity(:,iground,jcol))*lw_dn_g &
1133	& +this%lw_emissivity(:,iground,jcol)*this%planck_facet(:,iground,jcol))
1134	surface_flux%lw_abs_canopy(jcol,jtile) = flux%lw_dn(jcol,isurf) &
1135	& - flux%lw_up(jcol,isurf) - surface_flux%lw_dn_facet(jcol,iground) &
1136	& + surface_flux%lw_up_facet(jcol,iground)
1137
1138	end if
1139	if (config%do_sw) then
1140	sw_dn_direct_g = this%sw_tra_dir_dir(:,iregion,jcol) &
1141	& * flux%sw_dn_direct_surf_canopy(:,jcol)
1142	! Note that the following is initially just the
1143	! upwelling due to scattering of the direct beam
1144	sw_up_g = sw_dn_direct_g * this%sw_albedo_direct(:,iground,jcol)
1145	sw_dn_diffuse_g &
1146	& = (this%sw_tra_dif(:,iregion,jcol)*flux%sw_dn_diffuse_surf_canopy(:,jcol) &
1147	& +this%sw_ref_dif(:,iregion,jcol)*sw_up_g &
1148	& +this%sw_tra_dir_dif(:,iregion,jcol)*flux%sw_dn_direct_surf_canopy(:,jcol)) &
1149	& / (1.0_jprb - this%sw_albedo_diffuse(:,iground,jcol) &
1150	& *this%sw_ref_dif(:,iregion,jcol))
1151	sw_up_g = sw_up_g + sw_dn_diffuse_g * this%sw_albedo_diffuse(:,iground,jcol)
1152	surface_flux%sw_dn_direct_facet(jcol,iground) = sum(sw_dn_direct_g)
1153	surface_flux%sw_dn_facet(jcol,iground) = surface_flux%sw_dn_direct_facet(jcol,iground) &
1154	& + sum(sw_dn_diffuse_g)
1155	surface_flux%sw_up_facet(jcol,iground) = sum(sw_up_g)
1156	surface_flux%sw_abs_canopy(jcol,jtile) &
1157	& = flux%sw_dn(jcol,isurf) - flux%sw_up(jcol,isurf) &
1158	& - surface_flux%sw_dn_facet(jcol,iground) &
1159	& + surface_flux%sw_up_facet(jcol,iground)
1160	end if
1161	case (ITileUrban3D)
1162	iground = surface%i_ground_facet(jcol,jtile)
1163	iroof = surface%i_roof_facet(jcol,jtile)
1164	iwall = surface%i_wall_facet(jcol,jtile)
1165	iregion = surface%i_region_1(jcol,jtile)
1166
1167	! We want wall fluxes per unit area of wall, not per
1168	! unit area of street
1169	wall_scaling = (1.0_jprb - surface%building_fraction(jcol,jtile)) &
1170	& / max(1.0e-4_jprb, surface%building_normalized_perimeter(jcol,jtile) &
1171	& * surface%canopy_depth(jcol,jtile))
1172
1173	if (config%do_sw) then
1174	! Roof fluxes
1175	surface_flux%sw_dn_facet(jcol,iroof) = flux%sw_dn(jcol,isurf)
1176	surface_flux%sw_dn_direct_facet(jcol,iroof) = flux%sw_dn_direct(jcol,isurf)
1177	surface_flux%sw_up_facet(jcol,iroof) &
1178	& = sum(this%sw_albedo_diffuse(:,iroof,jcol) &
1179	& * flux%sw_dn_diffuse_surf_canopy(:,jcol) &
1180	& + this%sw_albedo_direct (:,iroof,jcol) &
1181	& * flux%sw_dn_direct_surf_canopy (:,jcol))
1182
1183	! Ground fluxes
1184	sw_dn_direct_g = this%sw_tra_dir_dir(:,iregion,jcol) &
1185	& * flux%sw_dn_direct_surf_canopy(:,jcol)
1186	! Note that the following is initially just the
1187	! upwelling due to scattering of the direct beam
1188	sw_up_g = sw_dn_direct_g * this%sw_albedo_direct(:,iground,jcol)
1189	sw_dn_diffuse_g &
1190	& = (this%sw_tra_dif(:,iregion,jcol)*flux%sw_dn_diffuse_surf_canopy(:,jcol) &
1191	& +this%sw_ref_dif(:,iregion,jcol)*sw_up_g &
1192	& +this%sw_tra_dir_dif(:,iregion,jcol)*flux%sw_dn_direct_surf_canopy(:,jcol)) &
1193	& / (1.0_jprb - this%sw_albedo_diffuse(:,iground,jcol) &
1194	& *this%sw_ref_dif(:,iregion,jcol))
1195	sw_up_g = sw_up_g + sw_dn_diffuse_g * this%sw_albedo_diffuse(:,iground,jcol)
1196	surface_flux%sw_dn_direct_facet(jcol,iground) = sum(sw_dn_direct_g)
1197	surface_flux%sw_dn_facet(jcol,iground) = surface_flux%sw_dn_direct_facet(jcol,iground) &
1198	& + sum(sw_dn_diffuse_g)
1199	surface_flux%sw_up_facet(jcol,iground) = sum(sw_up_g)
1200
1201	! Wall fluxes
1202	sw_abs_g = flux%sw_dn_direct_surf_canopy(:,jcol)*this%sw_wall_abs_dir(:,iregion,jcol)
1203	surface_flux%sw_dn_direct_facet(jcol,iwall) &
1204	& = wall_scaling * sum(sw_abs_g / (1.0_jprb - this%sw_albedo_direct(:,iwall,jcol)))
1205	! Initially just the direct reflection
1206	surface_flux%sw_up_facet(jcol,iwall) = wall_scaling &
1207	& * sum(sw_abs_g * this%sw_albedo_direct(:,iwall,jcol) &
1208	& / (1.0_jprb - this%sw_albedo_direct(:,iwall,jcol)))
1209	! Initially just the direct absorption
1210	surface_flux%sw_abs_canopy(jcol,jtile) &
1211	& = sum(flux%sw_dn_direct_surf_canopy(:,jcol)*this%sw_air_abs_dir(:,iregion,jcol))
1212
1213	! Diffuse absorption within canopy
1214	sw_abs_g = flux%sw_dn_direct_surf_canopy(:,jcol) &
1215	& * (1.0-this%sw_albedo_direct_reg (:,iregion,jcol)) &
1216	& + flux%sw_dn_diffuse_surf_canopy(:,jcol) &
1217	& * (1.0-this%sw_albedo_diffuse_reg(:,iregion,jcol)) &
1218	& - sw_dn_direct_g - sw_dn_diffuse_g + sw_up_g - sw_abs_g
1219	! Add diffuse absorption
1220	surface_flux%sw_abs_canopy(jcol,jtile) = surface_flux%sw_abs_canopy(jcol,jtile) &
1221	& + sum(sw_abs_g * (1.0_jprb - this%sw_wall_abs_frac_dif(:,iregion,jcol)))
1222	! Add diffuse reflection
1223	surface_flux%sw_up_facet(jcol,iwall) = surface_flux%sw_up_facet(jcol,iwall) &
1224	& + wall_scaling * sum(sw_abs_g * this%sw_wall_abs_frac_dif(:,iregion,jcol) &
1225	& * this%sw_albedo_diffuse(:,iwall,jcol) &
1226	& / (1.0_jprb - this%sw_albedo_diffuse(:,iwall,jcol)))
1227	! Add diffuse into wall
1228	surface_flux%sw_dn_facet(jcol,iwall) = surface_flux%sw_dn_direct_facet(jcol,iwall) &
1229	& + wall_scaling * sum(sw_abs_g * this%sw_wall_abs_frac_dif(:,iregion,jcol) &
1230	& / (1.0_jprb - this%sw_albedo_diffuse(:,iwall,jcol)))
1231	!!$
1232	!!$ sw_direct_abs_g = flux%sw_dn_direct_surf_g (:,jcol) - sw_dn_direct_g
1233	!!$ sw_diffuse_abs_g = flux%sw_dn_diffuse_surf_g(:,jcol)*(1.0_jprb - this%sw_albedo_diffuse_reg(:,iregion,jcol)) &
1234	!!$ & + flux%sw_dn_direct_surf_g (:,jcol)*(1.0_jprb - this%sw_albedo_direct_reg (:,iregion,jcol)) &
1235	!!$ & - sw_dn_direct_g - sw_dn_diffuse_g + sw_up_g
1236	!!$
1237	!!$
1238	!!$ surface_flux%sw_abs_canopy(jcol,jtile) = sum((1.0_jprb - this%sw_wall_abs_frac(:,iregion,jcol))*sw_abs_g)
1239	!!$ surface_flux%sw_dn_facet(jcol,iwall) = sum(this%sw_wall_abs_frac(:,iregion,jcol)*sw_abs_g &
1240	!!$ & / (1.0_jprb - this%sw_albedo_diffuse(:,iwall,jcol)))
1241	!!$ surface_flux%sw_up_facet(jcol,iwall) = sum(this%sw_albedo_diffuse(:,iwall,jcol) &
1242	!!$ & this%sw_wall_abs_frac(:,iregion,jcol)sw_abs_g &
1243	!!$ & / (1.0_jprb - this%sw_albedo_diffuse(:,iwall,jcol)))
1244	end if
1245
1246	if (config%do_lw) then
1247	surface_flux%lw_dn_facet(jcol,iroof) = flux%lw_dn(jcol,isurf)
1248	surface_flux%lw_up_facet(jcol,iroof) &
1249	& = sum(this%lw_emissivity(:,iroof,jcol)*this%planck_facet(:,iroof,jcol) &
1250	& + (1.0_jprb - this%lw_emissivity(:,iroof,jcol))*flux%lw_dn_surf_canopy(:,jcol))
1251	lw_dn_g = (this%lw_transmittance(:,iregion,jcol)*flux%lw_dn_surf_canopy(:,jcol) &
1252	& +this%lw_reflectance (:,iregion,jcol) &
1253	& this%lw_emissivity(:,iground,jcol)this%planck_facet(:,iground,jcol) &
1254	& +this%lw_source(:,iregion,jcol)) &
1255	& / (1.0_jprb - (1.0_jprb - this%lw_emissivity(:,iground,jcol)) &
1256	& *this%lw_reflectance(:,iregion,jcol))
1257	lw_up_g = (1.0_jprb-this%lw_emissivity(:,iground,jcol))*lw_dn_g &
1258	& +this%lw_emissivity(:,iground,jcol)*this%planck_facet(:,iground,jcol)
1259	surface_flux%lw_dn_facet(jcol,iground) = sum(lw_dn_g)
1260	surface_flux%lw_up_facet(jcol,iground) = sum(lw_up_g)
1261	! lw_abs_g = flux%lw_dn(jcol,isurf)-flux%lw_up(jcol,isurf) &
1262	! & - surface_flux%lw_dn_facet(jcol,iground)+surface_flux%lw_up_facet(jcol,iground)
1263	! lw_abs_g = flux%lw_dn_surf_canopy(:,jcol) &
1264	! & * (1.0_jprb-this%lw_emissivity_region(:,iregion,jcol)) &
1265	! & - this%lw_emission_region(:,iregion,jcol) &
1266	! & - lw_dn_g*(1.0_jprb-this%lw_emissivity(:,iground,jcol)) &
1267	! & + this%lw_emissivity(:,iground,jcol)*this%planck_facet(:,iground,jcol)
1268	! lw_abs_g = flux%lw_dn_surf_canopy(:,jcol) &
1269	! & * this%lw_emissivity_region(:,iregion,jcol) &
1270	! & - this%lw_emission_region(:,iregion,jcol) &
1271	! & - lw_dn_g*this%lw_emissivity(:,iground,jcol) &
1272	! & + this%lw_emissivity(:,iground,jcol)*this%planck_facet(:,iground,jcol)
1273	lw_abs_g = (flux%lw_dn_surf_canopy(:,jcol) + lw_up_g) &
1274	& * (1.0_jprb - this%lw_reflectance(:,iregion,jcol) &
1275	& - this%lw_transmittance(:,iregion,jcol)) &
1276	& + this%lw_total_wall_emission(:,iregion,jcol) &
1277	& + this%lw_total_canopy_emission(:,iregion,jcol) &
1278	& - 2.0_jprb * this%lw_source(:,iregion,jcol)
1279
1280	! surface_flux%lw_up_facet(jcol,iwall) = sum(this%planck_facet(:,iwall,jcol) &
1281	! & + wall_scaling*(1.0_jprb-this%lw_emissivity(:,iwall,jcol)) &
1282	! & this%lw_wall_abs_frac(:,iregion,jcol)lw_abs_g &
1283	! & / this%lw_emissivity(:,iwall,jcol))
1284	! surface_flux%lw_dn_facet(jcol,iwall) = sum((this%planck_facet(:,iwall,jcol) &
1285	! & + this%lw_total_wall_emission(:,iregion,jcol)) / this%lw_emissivity(:,iwall,jcol))
1286
1287	! surface_flux%lw_dn_facet(jcol,iwall) = wall_scaling*sum((this%lw_total_wall_emission(:,iregion,jcol) &
1288	! & + this%lw_wall_abs_frac(:,iregion,jcol)*lw_abs_g) &
1289	! & / this%lw_emissivity(:,iwall,jcol))
1290	! surface_flux%lw_up_facet(jcol,iwall) = surface_flux%lw_dn_facet(jcol,iwall) &
1291	! & - wall_scalingsum(this%lw_wall_abs_frac(:,iregion,jcol)lw_abs_g)
1292
1293	surface_flux%lw_dn_facet(jcol,iwall) &
1294	& = wall_scalingsum(this%lw_wall_abs_frac(:,iregion,jcol)lw_abs_g &
1295	& / this%lw_emissivity(:,iwall,jcol))
1296	surface_flux%lw_up_facet(jcol,iwall) = surface_flux%lw_dn_facet(jcol,iwall) &
1297	& + wall_scaling * sum(this%lw_total_wall_emission(:,iregion,jcol) &
1298	& - this%lw_wall_abs_frac(:,iregion,jcol)*lw_abs_g)
1299
1300	surface_flux%lw_abs_canopy(jcol,jtile) &
1301	& = sum(lw_abs_g*(1.0_jprb-this%lw_wall_abs_frac(:,iregion,jcol)) &
1302	& -this%lw_total_canopy_emission(:,iregion,jcol))
1303	end if
1304	end select
1305
1306	end if
1307
1308	end do
1309	end do
1310	end if
1311
1312	if (lhook) call dr_hook('radsurf_intermediate:partition_fluxes',1,hook_handle)
1313
1314	end subroutine partition_fluxes
1315
1316	end module radsurf_intermediate

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

Download in other formats:

Original Format