source: LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_overlap.F90 @ 5408

Last change on this file since 5408 was 4773, checked in by idelkadi, 12 months ago
  • Update of Ecrad in LMDZ The same organization of the Ecrad offline version is retained in order to facilitate the updating of Ecrad in LMDZ and the comparison between online and offline results. version 1.6.1 of Ecrad (https://github.com/lguez/ecrad.git)
  • Implementation of the double call of Ecrad in LMDZ


File size: 17.0 KB
Line 
1! radiation_overlap.F90 - Module to compute cloud overlap quantities
2!
3! (C) Copyright 2014- ECMWF.
4!
5! This software is licensed under the terms of the Apache Licence Version 2.0
6! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
7!
8! In applying this licence, ECMWF does not waive the privileges and immunities
9! granted to it by virtue of its status as an intergovernmental organisation
10! nor does it submit to any jurisdiction.
11!
12! Author:  Robin Hogan
13! Email:   r.j.hogan@ecmwf.int
14!
15! Modifications
16!   2017-10-23  R. Hogan  Renamed single-character variables
17!   2018-10-05  R. Hogan  Generalized alpha overlap for non-equal regions
18!   2018-10-08  R. Hogan  Removed calc_region_fractions
19
20module radiation_overlap
21
22  implicit none
23
24  public :: calc_overlap_matrices
25
26contains
27
28
29  ! This function now superceded by calc_region_properties in module
30  ! radiation_regions
31  ! !---------------------------------------------------------------------
32  ! ! Return an array of length nreg containing the fraction of the
33  ! ! gridbox occupied by each region for the specified cloud fraction.
34  ! pure function calc_region_fractions(nreg, cloud_fraction)
35
36  !   use parkind1, only : jprb
37
38  !   integer,    intent(in)      :: nreg
39  !   real(jprb), intent(in)      :: cloud_fraction
40
41  !   real(jprb), dimension(nreg) :: calc_region_fractions
42  !   integer :: jreg
43
44  !   if (nreg == 1) then
45  !     ! Only one region: must occupy all of gridbox
46  !     calc_region_fractions(1) = 1.0_jprb
47  !   else
48  !     ! Two or more regions: the first is the cloud-free region
49  !     calc_region_fractions(1) = 1.0_jprb - cloud_fraction
50
51  !     do jreg = 2,nreg
52  !       ! The cloudy regions are assumed to each have the same
53  !       ! fraction - see Shonk and Hogan (2008) for justification
54  !       calc_region_fractions(jreg) = cloud_fraction / (nreg - 1.0_jprb)
55  !     end do
56  !   end if
57
58  ! end function calc_region_fractions
59
60  !---------------------------------------------------------------------
61  ! Calculate a matrix expressing the overlap of regions in adjacent
62  ! layers, using the method of Shonk et al. (2010) in terms of their
63  ! "beta" overlap parameter
64  pure function calc_beta_overlap_matrix(nreg, op, frac_upper, frac_lower, &
65       &  frac_threshold) result(overlap_matrix)
66
67    use parkind1, only : jprb
68   
69    integer, intent(in) :: nreg ! Number of regions
70
71    ! Overlap parameter for each region, and fraction of the gridbox
72    ! occupied by each region in the upper and lower layers
73    real(jprb), intent(in), dimension(nreg) :: op, frac_upper, frac_lower
74
75    ! Cloud-fraction threshold below which cloud is deemed not to be
76    ! present
77    real(jprb), intent(in) :: frac_threshold
78
79    ! Output overlap matrix
80    real(jprb) :: overlap_matrix(nreg,nreg)
81
82    ! Denominator and its reciprocal in computing the random part of
83    ! the overlap matrix
84    real(jprb) :: denominator, factor
85
86    ! Beta overlap parameter multiplied by the minimum region fraction
87    ! of the upper and lower layers
88    real(jprb) :: op_x_frac_min(nreg)
89
90    integer :: jupper, jlower, jreg
91
92    ! In computing the random part of the overlap matrix we need
93    ! to divide all elements by "denominator", or for efficiency
94    ! multiply by "factor"
95    denominator = 1.0_jprb
96    do jreg = 1,nreg
97      op_x_frac_min(jreg) = op(jreg) &
98           &  * min(frac_upper(jreg), frac_lower(jreg))
99      denominator = denominator - op_x_frac_min(jreg)
100    end do
101    ! In principle the denominator can be zero
102    if (denominator >= frac_threshold) then
103      factor = 1.0_jprb / denominator
104      ! Create the random part of the overlap matrix
105      do jupper = 1,nreg
106        do jlower = 1,nreg
107          overlap_matrix(jupper,jlower) = factor &
108               &  * (frac_lower(jlower)-op_x_frac_min(jlower)) &
109               &  * (frac_upper(jupper)-op_x_frac_min(jupper))
110        end do
111      end do
112    else
113      overlap_matrix = 0.0_jprb
114    end if
115   
116    ! Add on the maximum part of the overlap matrix
117    do jreg = 1,nreg
118      overlap_matrix(jreg,jreg) = overlap_matrix(jreg,jreg) &
119           &  + op_x_frac_min(jreg)
120    end do
121
122  end function calc_beta_overlap_matrix
123
124
125  !---------------------------------------------------------------------
126  ! Calculate a matrix expressing the overlap of regions in adjacent
127  ! layers, using the Hogan and Illingworth (2000) "alpha" overlap
128  ! parameter, but allowing for the two cloudy regions in the
129  ! Tripleclouds assumption to have different areas
130  pure function calc_alpha_overlap_matrix(nreg, op, op_inhom, &
131       &  frac_upper, frac_lower) result(overlap_matrix)
132
133    use parkind1, only : jprb
134   
135    integer, intent(in) :: nreg ! Number of regions
136
137    ! Overlap parameter for cloud boundaries and for internal
138    ! inhomogeneities
139    real(jprb), intent(in) :: op, op_inhom
140
141    ! Fraction of the gridbox occupied by each region in the upper and
142    ! lower layers
143    real(jprb), intent(in), dimension(nreg) :: frac_upper, frac_lower
144
145    ! Output overlap matrix
146    real(jprb) :: overlap_matrix(nreg,nreg)
147
148    ! Combined cloud cover of pair of layers
149    real(jprb) :: pair_cloud_cover
150
151    ! Cloud fraction of upper and lower layers
152    real(jprb) :: cf_upper, cf_lower
153
154    ! One divided by cloud fraction
155    real(jprb) :: one_over_cf
156
157    ! Fraction of domain with cloud in both layers
158    real(jprb) :: frac_both
159
160    cf_upper = sum(frac_upper(2:nreg))
161    cf_lower = sum(frac_lower(2:nreg))
162
163    pair_cloud_cover = op*max(cf_upper,cf_lower) &
164           &  + (1.0_jprb - op) &
165           &  * (cf_upper+cf_lower-cf_upper*cf_lower)
166
167    ! Clear in both layers
168    overlap_matrix(1,1) = 1.0_jprb - pair_cloud_cover
169    if (nreg == 2) then
170      ! Clear in upper layer, cloudy in lower layer
171      overlap_matrix(1,2) = pair_cloud_cover - cf_upper
172      ! Clear in lower layer, cloudy in upper layer
173      overlap_matrix(2,1) = pair_cloud_cover - cf_lower
174      ! Cloudy in both layers
175      overlap_matrix(2,2) = cf_upper + cf_lower - pair_cloud_cover
176    else
177       ! Clear in upper layer, cloudy in lower layer
178      one_over_cf = 1.0_jprb / max(cf_lower, 1.0e-6_jprb)
179      overlap_matrix(1,2) = (pair_cloud_cover - cf_upper) &
180           &              * frac_lower(2) * one_over_cf
181      overlap_matrix(1,3) = (pair_cloud_cover - cf_upper) &
182           &              * frac_lower(3) * one_over_cf
183      ! Clear in lower layer, cloudy in upper layer
184      one_over_cf = 1.0_jprb / max(cf_upper, 1.0e-6_jprb)
185      overlap_matrix(2,1) = (pair_cloud_cover - cf_lower) &
186           &              * frac_upper(2) * one_over_cf
187      overlap_matrix(3,1) = (pair_cloud_cover - cf_lower) &
188           &              * frac_upper(3) * one_over_cf
189      ! Cloudy in both layers: frac_both is the fraction of the
190      ! gridbox with cloud in both layers
191      frac_both = cf_upper + cf_lower - pair_cloud_cover
192      ! Treat low and high optical-depth regions within frac_both as
193      ! one treats clear and cloudy skies in the whole domain;
194      ! redefine the following variables treating the high
195      ! optical-depth region as the cloud
196      cf_upper = frac_upper(3) / max(cf_upper, 1.0e-6_jprb)
197      cf_lower = frac_lower(3) / max(cf_lower, 1.0e-6_jprb)
198      pair_cloud_cover = op_inhom*max(cf_upper,cf_lower) &
199           &  + (1.0_jprb - op_inhom) &
200           &  * (cf_upper+cf_lower-cf_upper*cf_lower)
201      ! Assign overlaps for this 2x2 section of the 3x3 matrix as for
202      ! the 2-region case above, but multiplied by frac_both
203      overlap_matrix(2,2) = frac_both * (1.0_jprb - pair_cloud_cover)
204      overlap_matrix(2,3) = frac_both * (pair_cloud_cover - cf_upper)
205      overlap_matrix(3,2) = frac_both * (pair_cloud_cover - cf_lower)
206      overlap_matrix(3,3) = frac_both * (cf_upper+cf_lower-pair_cloud_cover)
207    end if
208
209  end function calc_alpha_overlap_matrix
210
211
212  !---------------------------------------------------------------------
213  ! Calculate a matrix expressing the overlap of regions in adjacent
214  ! layers, using the Hogan and Illingworth (2000) "alpha" overlap
215  ! parameter, and assuming the two cloudy regions in the Tripleclouds
216  ! assumption have the same area
217  pure function calc_alpha_overlap_matrix_simple(nreg, op, op_inhom, &
218       &  cf_upper, cf_lower) result(overlap_matrix)
219
220    use parkind1, only : jprb
221   
222    integer, intent(in) :: nreg ! Number of regions
223
224    ! Overlap parameter for cloud boundaries and for internal
225    ! inhomogeneities
226    real(jprb), intent(in) :: op, op_inhom
227
228    ! Cloud fraction in the upper and lower layers
229    real(jprb), intent(in) :: cf_upper, cf_lower
230
231    ! Output overlap matrix
232    real(jprb) :: overlap_matrix(nreg,nreg)
233
234    ! Combined cloud cover of pair of layers
235    real(jprb) :: pair_cloud_cover
236
237    real(jprb) :: cloud_unit
238
239    pair_cloud_cover = op*max(cf_upper,cf_lower) &
240           &  + (1.0_jprb - op) &
241           &  * (cf_upper+cf_lower-cf_upper*cf_lower)
242
243    ! Clear in both layers
244    overlap_matrix(1,1) = 1.0_jprb - pair_cloud_cover
245    if (nreg == 2) then
246      ! Clear in upper layer, cloudy in lower layer
247      overlap_matrix(1,2) = pair_cloud_cover - cf_upper
248      ! Clear in lower layer, cloudy in upper layer
249      overlap_matrix(2,1) = pair_cloud_cover - cf_lower
250      ! Cloudy in both layers
251      overlap_matrix(2,2) = cf_upper + cf_lower - pair_cloud_cover
252    else
253      ! The following assumes that the two cloudy regions are of equal area.
254      ! Clear in upper layer, cloudy in lower layer
255      overlap_matrix(1,2) = 0.5_jprb * (pair_cloud_cover - cf_upper)
256      overlap_matrix(1,3) = overlap_matrix(1,2)
257      ! Clear in lower layer, cloudy in upper layer
258      overlap_matrix(2,1) = 0.5_jprb * (pair_cloud_cover - cf_lower)
259      overlap_matrix(3,1) = overlap_matrix(2,1)
260      ! Cloudy in both layers
261      cloud_unit = 0.25_jprb * (cf_upper + cf_lower - pair_cloud_cover)
262      overlap_matrix(2,2) = cloud_unit * (1.0_jprb + op_inhom)
263      overlap_matrix(2,3) = cloud_unit * (1.0_jprb - op_inhom)
264      overlap_matrix(3,3) = overlap_matrix(2,2)
265      overlap_matrix(3,2) = overlap_matrix(2,3)
266    end if
267
268  end function calc_alpha_overlap_matrix_simple
269
270
271  !---------------------------------------------------------------------
272  ! Compute the upward and downward overlap matrices u_matrix and
273  ! v_matrix, respectively, where u_matrix is defined such that
274  ! y=u_matrix*x, where x is a vector of upwelling fluxes in each
275  ! region just below an interface, and y is a vector of upwelling
276  ! fluxes in each region just above that interface. For nlev model
277  ! levels there are nlev+1 interfaces including the ground and
278  ! top-of-atmosphere, and so that is one of the dimensions of
279  ! u_matrix and v_matrix.
280  subroutine calc_overlap_matrices(nlev,nreg,istartcol,iendcol, &
281       &     region_fracs, overlap_param, u_matrix, v_matrix, decorrelation_scaling, &
282       &     cloud_fraction_threshold, cloud_cover, use_beta_overlap)
283
284    use parkind1,     only : jprb
285    use yomhook,      only : lhook, dr_hook, jphook
286
287    ! Number of levels and regions
288    integer,  intent(in) :: nlev, nreg
289
290    ! Range of columns to process (also outer dimensions of u_matrix
291    ! and v_matrix)
292    integer, intent(in) :: istartcol, iendcol
293
294    ! Area fraction of each region: region 1 is clear sky, and 2+ are
295    ! the cloudy regions (only one or two cloudy regions are
296    ! supported)
297    real(jprb), intent(in), dimension(1:nreg,nlev,istartcol:iendcol)  :: region_fracs
298
299    ! The overlap parameter: either the "alpha" of Hogan & Illingworth
300    ! (2000) or the "beta" of Shonk et al. (2010)
301    real(jprb), intent(in), dimension(:,:)  :: overlap_param  ! (ncol,nlev-1)
302
303    ! Output overlap matrices
304    real(jprb), intent(out), dimension(nreg,nreg,nlev+1,istartcol:iendcol) &
305         &  :: u_matrix, v_matrix
306
307    ! For regions 2 and above, the overlap decorrelation length for
308    ! cloud boundaries is scaled by this amount to obtain the overlap
309    ! decorrelation length for cloud inhomogeneities. Typically this
310    ! number is 0.5, but if omitted it will be assumed to be one (same
311    ! decorrelation for cloud boundaries and in-cloud inhomogeneities)
312    real(jprb), intent(in), optional :: decorrelation_scaling
313
314    ! Regions smaller than this are ignored
315    real(jprb), intent(in), optional :: cloud_fraction_threshold
316
317    ! The diagnosed cloud cover is an optional output
318    real(jprb), intent(out), optional :: cloud_cover(:)
319
320    ! Do we use Shonk et al.'s (2010) "beta" overlap parameter?
321    logical, intent(in), optional :: use_beta_overlap
322
323    ! Loop indices for column, level, region and the regions in the
324    ! upper and lower layers for an interface
325    integer  :: jcol, jlev, jupper, jlower
326
327    ! Overlap matrix (non-directional)
328    real(jprb) :: overlap_matrix(nreg,nreg)
329
330    ! Fraction of the gridbox occupied by each region in the upper and
331    ! lower layers for an interface
332    real(jprb) :: frac_upper(nreg), frac_lower(nreg)
333
334    ! Beta overlap parameter for each region
335    real(jprb) :: op(nreg)
336
337    ! In case the user doesn't supply cloud_fraction_threshold we use
338    ! a default value
339    real(jprb) :: frac_threshold
340
341    ! The decorrelation scaling to use, in case decorrelation_scaling
342    ! was not provided
343    real(jprb) :: used_decorrelation_scaling
344
345    logical :: use_beta_overlap_param
346
347    real(jphook) :: hook_handle
348
349    if (lhook) call dr_hook('radiation_overlap:calc_overlap_matrices',0,hook_handle)
350
351    if (present(decorrelation_scaling)) then
352       used_decorrelation_scaling = decorrelation_scaling
353    else
354       used_decorrelation_scaling = 1.0_jprb
355    end if
356
357    if (present(cloud_fraction_threshold)) then
358      frac_threshold = cloud_fraction_threshold
359    else
360      frac_threshold = 1.0e-20_jprb
361    end if
362
363    if (present(use_beta_overlap)) then
364      use_beta_overlap_param = use_beta_overlap
365    else
366      use_beta_overlap_param = .false.
367    end if
368
369    ! Loop through each atmospheric column
370    do jcol = istartcol, iendcol
371      ! For this column, outer space is treated as one clear-sky
372      ! region, so the fractions are assigned as such
373      frac_upper(1) = 1.0_jprb
374      frac_upper(2:nreg) = 0.0_jprb
375
376      ! Overlap parameter is irrelevant when there is only one region
377      ! in the upper layer
378      op = 1.0_jprb
379
380      ! Loop down through the atmosphere, where jlev indexes each
381      ! half-level starting at 1 for the top-of-atmosphere, as well
382      ! as indexing each level starting at 1 for the top-most level.
383      do jlev = 1,nlev+1
384        ! Fraction of each region just below the interface
385        if (jlev > nlev) then
386          ! We are at the surface: treat as a single clear-sky
387          ! region
388          frac_lower(1) = 1.0_jprb
389          frac_lower(2:nreg) = 0.0_jprb
390        else
391          frac_lower = region_fracs(1:nreg,jlev,jcol)
392        end if
393   
394        ! Compute the overlap parameter of the interface just below
395        ! the current full level
396        if (jlev == 1 .or. jlev > nlev) then
397          ! We are at the surface or top-of-atmosphere: overlap
398          ! parameter is irrelevant
399          op = 1.0_jprb
400        else
401          ! We are not at the surface
402          op(1) = overlap_param(jcol,jlev-1)
403          ! For cloudy regions, scale the cloud-boundary overlap
404          ! parameter to obtain the cloud-inhomogeneity overlap
405          ! parameter as follows
406          if (op(1) >= 0.0_jprb) then
407            op(2:nreg) = op(1)**(1.0_jprb/used_decorrelation_scaling)
408          else
409            op(2:nreg) = op(1)
410          end if
411        end if
412     
413        if (use_beta_overlap_param) then
414          overlap_matrix = calc_beta_overlap_matrix(nreg, op, &
415               &  frac_upper, frac_lower, frac_threshold)
416        else
417          ! Simpler scheme assuming the two cloudy regions have the
418          ! same fraction
419          !overlap_matrix = calc_alpha_overlap_matrix_simple(nreg, &
420          !     &  op(1), op(2), &
421          !     &  1.0_jprb - frac_upper(1), 1.0_jprb - frac_lower(1))
422          ! More general scheme
423          overlap_matrix = calc_alpha_overlap_matrix(nreg, &
424               &  op(1), op(2), frac_upper, frac_lower)
425        end if
426
427        ! Convert to directional overlap matrices
428        do jupper = 1,nreg
429          do jlower = 1,nreg
430            if (frac_lower(jlower) >= frac_threshold) then
431              u_matrix(jupper,jlower,jlev,jcol) = overlap_matrix(jupper,jlower) &
432                   &  / frac_lower(jlower)
433            else
434              u_matrix(jupper,jlower,jlev,jcol) = 0.0_jprb
435            end if
436            if (frac_upper(jupper) >= frac_threshold) then
437              v_matrix(jlower,jupper,jlev,jcol) = overlap_matrix(jupper,jlower) &
438                   &  / frac_upper(jupper)
439            else
440              v_matrix(jlower,jupper,jlev,jcol) = 0.0_jprb
441            end if
442          end do
443        end do
444        frac_upper = frac_lower
445       
446      end do ! levels
447
448      ! Compute cloud cover from one of the directional overlap matrices
449      if (present(cloud_cover)) then
450        cloud_cover(jcol) = 1.0_jprb - product(v_matrix(1,1,:,jcol))
451      end if
452
453    end do ! columns
454
455    if (lhook) call dr_hook('radiation_overlap:calc_overlap_matrices',1,hook_handle)
456
457  end subroutine calc_overlap_matrices
458 
459end module radiation_overlap
Note: See TracBrowser for help on using the repository browser.