source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_cloud.F90 @ 4443

Last change on this file since 4443 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: 25.3 KB
Line 
1! radiation_cloud.F90 - Derived type to store cloud/precip properties
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!   2019-01-14  R. Hogan  Added inv_inhom_effective_size variable
17!   2019-01-14  R. Hogan  Added out_of_physical_bounds routine
18
19module radiation_cloud
20
21  use parkind1, only : jprb
22
23  implicit none
24  public
25
26  !---------------------------------------------------------------------
27  ! The intention is that all variables describing clouds and
28  ! radiatively-active precipitation are contained in this derived
29  ! type, and if cloud variables are to be added in future, they can
30  ! be added to this type without requiring extra variables to be
31  ! passed between subroutines elsewhere in the program.
32  type cloud_type
33    ! For maximum flexibility, an arbitrary number "ntype" of
34    ! cloud types could be stored, as follows:
35    !     integer :: ntype     ! number of cloud types
36    !     integer :: nfraction ! number of cloud fractions
37    !     real(jprb), allocatable, dimension(:,:,:) :: &
38    !          mixing_ratio, & ! (ncol,nwetlev,ntype) mass mixing ratio (kg/kg)
39    !          particle_size,& ! (ncol,nwetlev,ntype) effective radius/size (m)
40    !          fraction        ! (ncol,nwetlev,nfraction) areal (i.e. cloud) fraction
41    ! However, for practical purposes at the moment we consider two
42    ! cloud types, liquid cloud droplets and ice cloud
43    ! particles.  The following variables are dimensioned (ncol,nlev)
44    real(jprb), allocatable, dimension(:,:) :: &
45         &  q_liq,  q_ice,  & ! mass mixing ratio (kg/kg)
46         &  re_liq, re_ice, & ! effective radius (m)
47         &  fraction          ! (0-1) Assume liq & ice completely mixed
48
49    ! The fractional standard deviation of cloud optical depth in the
50    ! cloudy part of the gridbox.  In the Tripleclouds representation
51    ! of cloud inhomogeneity, this is implemented by splitting the
52    ! cloudy part of the gridbox into two equal-area regions, one
53    ! with the cloud optical depth scaled by 1+fractional_std and the
54    ! other scaled by 1-fractional_std. This variable is dimensioned
55    ! (ncol,nlev)
56    real(jprb), allocatable, dimension(:,:) :: fractional_std
57
58    ! The inverse of the effective horizontal size of the clouds in
59    ! the gridbox, used to compute the cloud edge length per unit
60    ! gridbox area for use in representing 3D effects. This variable
61    ! is dimensioned (ncol,nlev).
62    real(jprb), allocatable, dimension(:,:) :: inv_cloud_effective_size ! m-1
63
64    ! Similarly for the in-cloud heterogeneities, used to compute the
65    ! edge length between the optically thin and thick cloudy regions
66    ! of the gridbox.
67    real(jprb), allocatable, dimension(:,:) :: inv_inhom_effective_size ! m-1
68
69    ! The following variable describes the overlap of cloud boundaries
70    ! in adjacent layers, with dimensions (ncol,nlev-1): 1 corresponds
71    ! to maximum overlap and 0 to random overlap. Depending on the
72    ! ecRad configuration, it may be the "alpha" overlap parameter of
73    ! Hogan and Illingworth (2000) or the "beta" overlap parameter of
74    ! Shonk et al. (2010).
75    real(jprb), allocatable, dimension(:,:) :: overlap_param
76
77  contains
78    procedure :: allocate   => allocate_cloud_arrays
79    procedure :: deallocate => deallocate_cloud_arrays
80    procedure :: set_overlap_param
81    procedure :: set_overlap_param_approx
82    procedure :: create_fractional_std
83    procedure :: create_inv_cloud_effective_size
84    procedure :: create_inv_cloud_effective_size_eta
85    procedure :: param_cloud_effective_separation_eta
86    procedure :: crop_cloud_fraction
87    procedure :: out_of_physical_bounds
88
89  end type cloud_type
90
91contains
92
93  !---------------------------------------------------------------------
94  ! Allocate arrays for describing clouds and precipitation, although
95  ! in the offline code these are allocated when they are read from
96  ! the NetCDF file
97  subroutine allocate_cloud_arrays(this, ncol, nlev, use_inhom_effective_size)
98
99    use yomhook,     only : lhook, dr_hook
100
101    class(cloud_type), intent(inout) :: this
102    integer, intent(in)              :: ncol  ! Number of columns
103    integer, intent(in)              :: nlev  ! Number of levels
104    logical, intent(in), optional    :: use_inhom_effective_size
105
106    real(jprb)                       :: hook_handle
107
108    if (lhook) call dr_hook('radiation_cloud:allocate',0,hook_handle)
109
110    allocate(this%q_liq(ncol,nlev))
111    allocate(this%re_liq(ncol,nlev))
112    allocate(this%q_ice(ncol,nlev))
113    allocate(this%re_ice(ncol,nlev))
114    allocate(this%fraction(ncol,nlev))
115    allocate(this%overlap_param(ncol,nlev-1))
116    allocate(this%fractional_std(ncol,nlev))
117    allocate(this%inv_cloud_effective_size(ncol,nlev))
118
119    if (present(use_inhom_effective_size)) then
120      if (use_inhom_effective_size) then
121        allocate(this%inv_inhom_effective_size(ncol,nlev))
122      end if
123    end if
124
125    if (lhook) call dr_hook('radiation_cloud:allocate',1,hook_handle)
126
127  end subroutine allocate_cloud_arrays
128
129
130  !---------------------------------------------------------------------
131  ! Deallocate arrays
132  subroutine deallocate_cloud_arrays(this)
133
134    use yomhook,     only : lhook, dr_hook
135
136    class(cloud_type), intent(inout) :: this
137
138    real(jprb)                       :: hook_handle
139
140    if (lhook) call dr_hook('radiation_cloud:deallocate',0,hook_handle)
141
142    if (allocated(this%q_liq))    deallocate(this%q_liq)
143    if (allocated(this%re_liq))   deallocate(this%re_liq)
144    if (allocated(this%q_ice))    deallocate(this%q_ice)
145    if (allocated(this%re_ice))   deallocate(this%re_ice)
146    if (allocated(this%fraction)) deallocate(this%fraction)
147    if (allocated(this%overlap_param))  deallocate(this%overlap_param)
148    if (allocated(this%fractional_std)) deallocate(this%fractional_std)
149    if (allocated(this%inv_cloud_effective_size)) &
150         &  deallocate(this%inv_cloud_effective_size)
151    if (allocated(this%inv_inhom_effective_size)) &
152         &  deallocate(this%inv_inhom_effective_size)
153
154    if (lhook) call dr_hook('radiation_cloud:deallocate',1,hook_handle)
155
156  end subroutine deallocate_cloud_arrays
157
158
159  !---------------------------------------------------------------------
160  ! Compute and store the overlap parameter from the provided overlap
161  ! decorrelation length (in metres).  If istartcol and/or iendcol are
162  ! provided then only columns in this range are computed.  If the
163  ! overlap_param array has not been allocated then it will be
164  ! allocated to be of the correct size relative to the pressure
165  ! field.
166  subroutine set_overlap_param(this, thermodynamics, decorrelation_length, &
167       &  istartcol, iendcol)
168
169    use yomhook,                  only : lhook, dr_hook
170    use radiation_thermodynamics, only : thermodynamics_type
171    use radiation_constants,      only : GasConstantDryAir, AccelDueToGravity
172
173    class(cloud_type),         intent(inout) :: this
174    type(thermodynamics_type), intent(in)    :: thermodynamics
175    real(jprb),                intent(in)    :: decorrelation_length ! m
176    integer,         optional, intent(in)    :: istartcol, iendcol
177
178    ! Ratio of gas constant for dry air to acceleration due to gravity
179    real(jprb), parameter :: R_over_g = GasConstantDryAir / AccelDueToGravity
180
181    ! Process only columns i1 to i2, which will be istartcol to
182    ! iendcol if they were provided
183    integer :: i1, i2
184
185    integer :: ncol, nlev
186
187    integer :: jlev
188
189    real(jprb)        :: hook_handle
190
191    if (lhook) call dr_hook('radiation_cloud:set_overlap_param',0,hook_handle)
192
193    ! Pressure at half-levels, pressure_hl, is defined at nlev+1
194    ! points
195    ncol = size(thermodynamics%pressure_hl,dim=1)
196    nlev = size(thermodynamics%pressure_hl,dim=2)-1
197
198    if (present(istartcol)) then
199      i1 = istartcol
200    else
201      i1 = 1
202    end if
203
204    if (present(iendcol)) then
205      i2 = iendcol
206    else
207      i2 = ncol
208    end if
209
210    if (.not. allocated(this%overlap_param)) then
211      ! If pressure is of size (ncol,nlev+1) then overlap_param is of
212      ! size (ncol,nlev-1), since overlap parameter is only defined here
213      ! for interfaces between model layers, not for the interface to
214      ! space or the surface
215      allocate(this%overlap_param(ncol, nlev-1))
216    end if
217
218    if (thermodynamics%pressure_hl(i1,2) > thermodynamics%pressure_hl(i1,1)) then
219      ! Pressure is increasing with index (order of layers is
220      ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we
221      ! don't take the logarithm of the first pressure in each column.
222      this%overlap_param(i1:i2,1) = exp(-(R_over_g/decorrelation_length) &
223           &                            * thermodynamics%temperature_hl(i1:i2,2) &
224           &                            *log(thermodynamics%pressure_hl(i1:i2,3) &
225           &                                /thermodynamics%pressure_hl(i1:i2,2)))
226
227      do jlev = 2,nlev-1
228        this%overlap_param(i1:i2,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length) &
229             &                            * thermodynamics%temperature_hl(i1:i2,jlev+1) &
230             &                            *log(thermodynamics%pressure_hl(i1:i2,jlev+2) &
231             &                                /thermodynamics%pressure_hl(i1:i2,jlev)))
232      end do
233
234    else
235       ! Pressure is decreasing with index (order of layers is surface
236       ! to top-of-atmosphere).  In case pressure_hl(:,nlev+1)=0, we
237       ! don't take the logarithm of the last pressure in each column.
238      do jlev = 1,nlev-2
239        this%overlap_param(i1:i2,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length) &
240             &                            * thermodynamics%temperature_hl(i1:i2,jlev+1) &
241             &                            *log(thermodynamics%pressure_hl(i1:i2,jlev) &
242             &                                /thermodynamics%pressure_hl(i1:i2,jlev+2)))
243      end do
244      this%overlap_param(i1:i2,nlev-1) = exp(-(R_over_g/decorrelation_length) &
245           &                            * thermodynamics%temperature_hl(i1:i2,nlev) &
246           &                            *log(thermodynamics%pressure_hl(i1:i2,nlev-1) &
247           &                                /thermodynamics%pressure_hl(i1:i2,nlev)))
248
249    end if
250
251    if (lhook) call dr_hook('radiation_cloud:set_overlap_param',1,hook_handle)
252
253  end subroutine set_overlap_param
254
255
256  !---------------------------------------------------------------------
257  ! Compute and store the overlap parameter from the provided overlap
258  ! decorrelation length (in metres).  If istartcol and/or iendcol are
259  ! provided then only columns in this range are computed.  If the
260  ! overlap_param array has not been allocated then it will be
261  ! allocated to be of the correct size relative to the pressure
262  ! field. This is the APPROXIMATE method as it assumes a fixed
263  ! atmospheric scale height, which leads to differences particularly
264  ! in low cloud.
265  subroutine set_overlap_param_approx(this, thermodynamics, decorrelation_length, &
266       &  istartcol, iendcol)
267
268    use yomhook,                  only : lhook, dr_hook
269    use radiation_thermodynamics, only : thermodynamics_type
270
271    class(cloud_type),         intent(inout) :: this
272    type(thermodynamics_type), intent(in)    :: thermodynamics
273    real(jprb),                intent(in)    :: decorrelation_length ! m
274    integer,         optional, intent(in)    :: istartcol, iendcol
275
276    ! To convert decorrelation length (m) to overlap parameter between
277    ! layers, we need an estimate for the thickness of the layer. This
278    ! is found using the pressure difference between the edges of the
279    ! layer, along with the approximate scale height of the atmosphere
280    ! (m) given here:
281    real(jprb), parameter :: scale_height = 8000.0_jprb
282
283    ! Process only columns i1 to i2, which will be istartcol to
284    ! iendcol if they were provided
285    integer :: i1, i2
286
287    integer :: ncol, nlev
288
289    real(jprb)        :: hook_handle
290
291    if (lhook) call dr_hook('radiation_cloud:set_overlap_param_approx',0,hook_handle)
292
293    ! Pressure at half-levels, pressure_hl, is defined at nlev+1
294    ! points
295    ncol = size(thermodynamics%pressure_hl,dim=1)
296    nlev = size(thermodynamics%pressure_hl,dim=2)-1
297
298    if (present(istartcol)) then
299      i1 = istartcol
300    else
301      i1 = 1
302    end if
303
304    if (present(iendcol)) then
305      i2 = iendcol
306    else
307      i2 = ncol
308    end if
309
310    if (.not. allocated(this%overlap_param)) then
311      ! If pressure is of size (ncol,nlev+1) then overlap_param is of
312      ! size (ncol,nlev-1), since overlap parameter is only defined here
313      ! for interfaces between model layers, not for the interface to
314      ! space or the surface
315      allocate(this%overlap_param(ncol, nlev-1))
316    end if
317
318    if (thermodynamics%pressure_hl(i1,2) > thermodynamics%pressure_hl(i1,1)) then
319       ! Pressure is increasing with index (order of layers is
320       ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we
321       ! don't take the logarithm of the first pressure in each
322       ! column.
323       this%overlap_param(i1:i2,:) = exp(-(scale_height/decorrelation_length) &
324            &  * ( log(thermodynamics%pressure_hl(i1:i2,3:nlev+1) &
325            &         /thermodynamics%pressure_hl(i1:i2,2:nlev  )) ) )
326    else
327       ! Pressure is decreasing with index (order of layers is surface
328       ! to top-of-atmosphere).  In case pressure_hl(:,nlev+1)=0, we
329       ! don't take the logarithm of the last pressure in each column.
330       this%overlap_param(i1:i2,:) = exp(-(scale_height/decorrelation_length) &
331            &  * ( log(thermodynamics%pressure_hl(i1:i2,1:nlev-1) &
332            &         /thermodynamics%pressure_hl(i1:i2,2:nlev  )) ) )
333    end if
334
335    if (lhook) call dr_hook('radiation_cloud:set_overlap_param_approx',1,hook_handle)
336
337  end subroutine set_overlap_param_approx
338
339
340  !---------------------------------------------------------------------
341  ! Create a matrix of constant fractional standard deviations
342  ! (dimensionless)
343  subroutine create_fractional_std(this, ncol, nlev, frac_std)
344
345    use yomhook,                  only : lhook, dr_hook
346
347    class(cloud_type), intent(inout) :: this
348    integer,           intent(in)    :: ncol, nlev
349    real(jprb),        intent(in)    :: frac_std
350
351    real(jprb)             :: hook_handle
352
353    if (lhook) call dr_hook('radiation_cloud:create_fractional_std',0,hook_handle)
354
355    if (allocated(this%fractional_std)) then
356       deallocate(this%fractional_std)
357    end if
358   
359    allocate(this%fractional_std(ncol, nlev))
360
361    this%fractional_std = frac_std
362
363    if (lhook) call dr_hook('radiation_cloud:create_fractional_std',1,hook_handle)
364
365  end subroutine create_fractional_std
366
367
368  !---------------------------------------------------------------------
369  ! Create a matrix of constant inverse cloud effective size (m-1)
370  subroutine create_inv_cloud_effective_size(this, ncol, nlev, inv_eff_size)
371
372    use yomhook,                  only : lhook, dr_hook
373
374    class(cloud_type), intent(inout) :: this
375    integer,           intent(in)    :: ncol, nlev
376    real(jprb),        intent(in)    :: inv_eff_size
377
378    real(jprb)             :: hook_handle
379
380    if (lhook) call dr_hook('radiation_cloud:create_inv_cloud_effective_size',0,hook_handle)
381
382    if (allocated(this%inv_cloud_effective_size)) then
383       deallocate(this%inv_cloud_effective_size)
384    end if
385   
386    allocate(this%inv_cloud_effective_size(ncol, nlev))
387
388    this%inv_cloud_effective_size = inv_eff_size
389
390    if (lhook) call dr_hook('radiation_cloud:create_inv_cloud_effective_size',1,hook_handle)
391
392  end subroutine create_inv_cloud_effective_size
393
394
395  !---------------------------------------------------------------------
396  ! Create a matrix of inverse cloud effective size (m-1) according to
397  ! the value of eta (=pressure divided by surface pressure)
398  subroutine create_inv_cloud_effective_size_eta(this, ncol, nlev, &
399       &  pressure_hl, inv_eff_size_low, inv_eff_size_mid, inv_eff_size_high, &
400       &  eta_low_mid, eta_mid_high, istartcol, iendcol)
401
402    use yomhook,                  only : lhook, dr_hook
403
404    class(cloud_type), intent(inout) :: this
405    integer,           intent(in)    :: ncol, nlev
406    ! Pressure on half levels (Pa)
407    real(jprb),        intent(in)    :: pressure_hl(:,:)
408    ! Inverse effective size for low, mid and high cloud (m-1)
409    real(jprb),        intent(in)    :: inv_eff_size_low
410    real(jprb),        intent(in)    :: inv_eff_size_mid
411    real(jprb),        intent(in)    :: inv_eff_size_high
412    ! Eta values at low-mid and mid-high interfaces
413    real(jprb),        intent(in)    :: eta_low_mid, eta_mid_high
414    integer, optional, intent(in)    :: istartcol, iendcol
415
416    ! Ratio of layer midpoint pressure to surface pressure
417    real(jprb) :: eta(nlev)
418
419    ! Indices of column, level and surface half-level
420    integer :: jcol, isurf
421
422    ! Local values of istartcol, iendcol
423    integer :: i1, i2
424
425    real(jprb)             :: hook_handle
426
427    if (lhook) call dr_hook('radiation_cloud:create_inv_cloud_effective_size_eta',0,hook_handle)
428
429    if (allocated(this%inv_cloud_effective_size)) then
430      deallocate(this%inv_cloud_effective_size)
431    end if
432   
433    allocate(this%inv_cloud_effective_size(ncol, nlev))
434
435    if (present(istartcol)) then
436      i1 = istartcol
437    else
438      i1 = 1
439    end if
440
441    if (present(iendcol)) then
442      i2 = iendcol
443    else
444      i2 = ncol
445    end if
446
447    ! Locate the surface half-level
448    if (pressure_hl(1,1) > pressure_hl(1,2)) then
449      isurf = 1
450    else
451      isurf = nlev+1
452    end if
453
454    do jcol = i1,i2
455      eta = (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)) &
456           &  * (0.5_jprb / pressure_hl(jcol,isurf))
457      where (eta > eta_low_mid)
458        this%inv_cloud_effective_size(jcol,:) = inv_eff_size_low
459      elsewhere (eta > eta_mid_high)
460        this%inv_cloud_effective_size(jcol,:) = inv_eff_size_mid
461      elsewhere
462        this%inv_cloud_effective_size(jcol,:) = inv_eff_size_high
463      end where
464    end do
465
466    if (lhook) call dr_hook('radiation_cloud:create_inv_cloud_effective_size_eta',1,hook_handle)
467
468  end subroutine create_inv_cloud_effective_size_eta
469
470
471  !---------------------------------------------------------------------
472  ! Create a matrix of inverse cloud and inhomogeneity effective size
473  ! (m-1) parameterized according to the value of eta (=pressure
474  ! divided by surface pressure): effective_separation =
475  ! coeff_a + coeff_b*exp(-(eta**power)). 
476  subroutine param_cloud_effective_separation_eta(this, ncol, nlev, &
477       &  pressure_hl, separation_surf, separation_toa, power, &
478       &  inhom_separation_factor, istartcol, iendcol)
479
480    use yomhook,                  only : lhook, dr_hook
481
482    class(cloud_type), intent(inout) :: this
483    integer,           intent(in)    :: ncol, nlev
484    ! Pressure on half levels (Pa)
485    real(jprb),        intent(in)    :: pressure_hl(:,:)
486    ! Separation distances at surface and top-of-atmosphere, and power
487    ! on eta
488    real(jprb),           intent(in) :: separation_surf ! m
489    real(jprb),           intent(in) :: separation_toa ! m
490    real(jprb),           intent(in) :: power
491    real(jprb), optional, intent(in) :: inhom_separation_factor
492    integer,    optional, intent(in) :: istartcol, iendcol
493
494    ! Ratio of layer midpoint pressure to surface pressure
495    real(jprb) :: eta(nlev)
496
497    ! Effective cloud separation (m)
498    real(jprb) :: eff_separation(nlev)
499
500    ! Coefficients used to compute effective separation distance
501    real(jprb) :: coeff_e, coeff_a, coeff_b, inhom_sep_factor
502
503    ! Indices of column, level and surface half-level
504    integer :: jcol, isurf
505
506    ! Local values of istartcol, iendcol
507    integer :: i1, i2
508
509    real(jprb) :: hook_handle
510
511    if (lhook) call dr_hook('radiation_cloud:param_cloud_effective_separation_eta',0,hook_handle)
512
513    if (present(inhom_separation_factor)) then
514      inhom_sep_factor = inhom_separation_factor
515    else
516      inhom_sep_factor = 1.0_jprb
517    end if
518
519    coeff_e = 1.0_jprb - exp(-1.0_jprb)
520    coeff_b = (separation_toa - separation_surf) / coeff_e
521    coeff_a = separation_toa - coeff_b
522
523    if (allocated(this%inv_cloud_effective_size)) then
524      deallocate(this%inv_cloud_effective_size)
525    end if
526     if (allocated(this%inv_inhom_effective_size)) then
527      deallocate(this%inv_inhom_effective_size)
528    end if
529   
530    allocate(this%inv_cloud_effective_size(ncol, nlev))
531    allocate(this%inv_inhom_effective_size(ncol, nlev))
532
533    if (present(istartcol)) then
534      i1 = istartcol
535    else
536      i1 = 1
537    end if
538
539    if (present(iendcol)) then
540      i2 = iendcol
541    else
542      i2 = ncol
543    end if
544
545    ! Locate the surface half-level
546    if (pressure_hl(1,1) > pressure_hl(1,2)) then
547      isurf = 1
548    else
549      isurf = nlev+1
550    end if
551
552    do jcol = i1,i2
553      eta = (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)) &
554           &  * (0.5_jprb / pressure_hl(jcol,isurf))
555      eff_separation = coeff_a + coeff_b * exp(-eta**power)
556      this%inv_cloud_effective_size(jcol,:) = 1.0_jprb / (eff_separation &
557           &  * sqrt(max(1.0e-5_jprb,this%fraction(jcol,:)*(1.0_jprb-this%fraction(jcol,:)))))
558      this%inv_inhom_effective_size(jcol,:) = 1.0_jprb / (eff_separation * inhom_sep_factor &
559           &  * sqrt(max(1.0e-5_jprb,0.5_jprb*this%fraction(jcol,:)*(1.0_jprb-0.5_jprb*this%fraction(jcol,:)))))
560    end do
561
562    if (lhook) call dr_hook('radiation_cloud:param_cloud_effective_separation_eta',1,hook_handle)
563
564  end subroutine param_cloud_effective_separation_eta
565
566
567  !---------------------------------------------------------------------
568  ! Remove "ghost" clouds: those with a cloud fraction that is too
569  ! small to treat sensibly (e.g. because it implies that the
570  ! "in-cloud" water content is too high), or with a cloud water
571  ! content that is too small.  We do this in one place to ensure that
572  ! all subsequent subroutines can assume that if cloud_fraction > 0.0
573  ! then cloud is really present and should be treated.
574  subroutine crop_cloud_fraction(this, istartcol, iendcol, &
575       &    cloud_fraction_threshold, cloud_mixing_ratio_threshold)
576   
577    use yomhook, only : lhook, dr_hook
578
579    class(cloud_type), intent(inout) :: this
580    integer,           intent(in)    :: istartcol, iendcol
581
582    integer :: nlev
583    integer :: jcol, jlev
584
585    real(jprb) :: cloud_fraction_threshold, cloud_mixing_ratio_threshold
586
587    real(jprb) :: hook_handle
588
589    if (lhook) call dr_hook('radiation_cloud:crop_cloud_fraction',0,hook_handle)
590
591    nlev = size(this%fraction,2)
592
593    do jlev = 1,nlev
594      do jcol = istartcol,iendcol
595        if (this%fraction(jcol,jlev) < cloud_fraction_threshold &
596             &  .or. this%q_liq(jcol,jlev)+this%q_ice(jcol,jlev) &
597             &        < cloud_mixing_ratio_threshold) then
598          this%fraction(jcol,jlev) = 0.0_jprb
599        end if
600      end do
601    end do
602
603    if (lhook) call dr_hook('radiation_cloud:crop_cloud_fraction',1,hook_handle)
604
605  end subroutine crop_cloud_fraction
606
607
608  !---------------------------------------------------------------------
609  ! Return .true. if variables are out of a physically sensible range,
610  ! optionally only considering columns between istartcol and iendcol
611  function out_of_physical_bounds(this, istartcol, iendcol, do_fix) result(is_bad)
612
613    use yomhook,          only : lhook, dr_hook
614    use radiation_config, only : out_of_bounds_2d
615
616    class(cloud_type), intent(inout) :: this
617    integer,  optional,intent(in) :: istartcol, iendcol
618    logical,  optional,intent(in) :: do_fix
619    logical                       :: is_bad
620
621    logical    :: do_fix_local
622
623    real(jprb) :: hook_handle
624
625    if (lhook) call dr_hook('radiation_cloud:out_of_physical_bounds',0,hook_handle)
626
627    if (present(do_fix)) then
628      do_fix_local = do_fix
629    else
630      do_fix_local = .false.
631    end if
632
633    is_bad =    out_of_bounds_2d(this%q_liq, 'q_liq', 0.0_jprb, 1.0_jprb, &
634         &                       do_fix_local, i1=istartcol, i2=iendcol) &
635         & .or. out_of_bounds_2d(this%q_ice, 'q_ice', 0.0_jprb, 1.0_jprb, &
636         &                       do_fix_local, i1=istartcol, i2=iendcol) &
637         & .or. out_of_bounds_2d(this%re_liq, 're_liq', 0.0_jprb, 0.01_jprb, &
638         &                       do_fix_local, i1=istartcol, i2=iendcol) &
639         & .or. out_of_bounds_2d(this%re_ice, 're_ice', 0.0_jprb, 0.1_jprb, &
640         &                       do_fix_local, i1=istartcol, i2=iendcol) &
641         & .or. out_of_bounds_2d(this%fraction, 'cloud%fraction', 0.0_jprb, 1.0_jprb, &
642         &                       do_fix_local, i1=istartcol, i2=iendcol) &
643         & .or. out_of_bounds_2d(this%fractional_std, 'fractional_std', 0.0_jprb, 10.0_jprb, &
644         &                       do_fix_local, i1=istartcol, i2=iendcol) &
645         & .or. out_of_bounds_2d(this%inv_cloud_effective_size, 'inv_cloud_effective_size', &
646         &                       0.0_jprb, 1.0_jprb, do_fix_local, i1=istartcol, i2=iendcol) &
647         & .or. out_of_bounds_2d(this%inv_inhom_effective_size, 'inv_inhom_effective_size', &
648         &                       0.0_jprb, 1.0_jprb, do_fix_local, i1=istartcol, i2=iendcol) &
649         & .or. out_of_bounds_2d(this%overlap_param, 'overlap_param', -0.5_jprb, 1.0_jprb, &
650         &                       do_fix_local, i1=istartcol, i2=iendcol)
651
652    if (lhook) call dr_hook('radiation_cloud:out_of_physical_bounds',1,hook_handle)
653
654  end function out_of_physical_bounds
655 
656end module radiation_cloud
Note: See TracBrowser for help on using the repository browser.