source: LMDZ6/trunk/libf/phylmd/ecrad/radsurf_intermediate.F90 @ 4586

Last change on this file since 4586 was 3908, checked in by idelkadi, 3 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
16module 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
108contains
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             &  = StefanBoltzmann*surface%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_wall*this%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_wall*this%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_wall*this%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_scaling*sum(this%lw_wall_abs_frac(:,iregion,jcol)*lw_abs_g)
1292
1293                surface_flux%lw_dn_facet(jcol,iwall) &
1294                     &  = wall_scaling*sum(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
1316end module radsurf_intermediate
Note: See TracBrowser for help on using the repository browser.