source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_regions.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: 8.7 KB
Line 
1! radiation_regions.F90 -- Properties of horizontal regions in Tripleclouds & SPARTACUS
2!
3! (C) Copyright 2016- 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-07-14  R. Hogan  Incorporate gamma distribution option
17!   2018-10-06  R. Hogan  Merged from radiation_optical_depth_scaling.h and radiation_overlap.F90
18
19module radiation_regions
20
21  implicit none
22
23  public :: calc_region_properties
24
25contains
26
27  !---------------------------------------------------------------------
28  ! Compute the optical depth scalings for the optically "thick" and
29  ! "thin" regions of a Tripleclouds representation of a sub-grid PDF
30  ! of cloud optical depth. Following Shonk and Hogan (2008), the 16th
31  ! percentile is used for the thin region, and the formulas estimate
32  ! this for both lognormal and gamma distributions. However, an
33  ! adjustment is needed for the gamma distribution at large
34  ! fractional standard deviations.
35  subroutine calc_region_properties(nlev, nreg, istartcol, iendcol, do_gamma, &
36       &  cloud_fraction, frac_std, reg_fracs, od_scaling, cloud_fraction_threshold)
37
38    use parkind1,     only : jprb
39    use yomhook,      only : lhook, dr_hook
40    use radiation_io, only : nulerr, radiation_abort
41
42    ! Minimum od_scaling in the case of a Gamma distribution
43    real(jprb), parameter :: MinGammaODScaling = 0.025_jprb
44
45    ! At large fractional standard deviations (FSDs), we cannot
46    ! capture the behaviour of a gamma distribution with two equally
47    ! weighted points; we need to weight the first ("lower") point
48    ! more.  The weight of the first point is normally 0.5, but for
49    ! FSDs between 1.5 and 3.725 it rises linearly to 0.9, and for
50    ! higher FSD it is capped at this value.  The weight of the second
51    ! point is one minus the first point.
52    real(jprb), parameter :: MinLowerFrac      = 0.5_jprb
53    real(jprb), parameter :: MaxLowerFrac      = 0.9_jprb
54    real(jprb), parameter :: FSDAtMinLowerFrac = 1.5_jprb
55    real(jprb), parameter :: FSDAtMaxLowerFrac = 3.725_jprb
56    ! Between FSDAtMinLowerFrac and FSDAtMaxLowerFrac,
57    ! LowerFrac=LowerFracFSDIntercept+FSD*LowerFracFSDGradient
58    real(jprb), parameter :: LowerFracFSDGradient &
59         &  = (MaxLowerFrac-MinLowerFrac) / (FSDAtMaxLowerFrac-FSDAtMinLowerFrac)
60    real(jprb), parameter :: LowerFracFSDIntercept &
61         &  = MinLowerFrac - FSDAtMinLowerFrac*LowerFracFSDGradient
62
63    ! Number of levels and regions
64    integer, intent(in) :: nlev, nreg
65
66    ! Range of columns to process
67    integer, intent(in) :: istartcol, iendcol
68
69    ! Do we do a lognormal or gamma distribution?
70    logical, intent(in) :: do_gamma
71
72    ! Cloud fraction, i.e. the fraction of the gridbox assigned to all
73    ! regions numbered 2 and above (region 1 is clear sky)
74    real(jprb), intent(in), dimension(:,:)  :: cloud_fraction ! (ncol,nlev)
75
76    ! Fractional standard deviation of in-cloud water content
77    real(jprb), intent(in), dimension(:,:)  :: frac_std       ! (ncol,nlev)
78
79    ! Fractional area coverage of each region
80    real(jprb), intent(out) :: reg_fracs(1:nreg,nlev,istartcol:iendcol)
81
82    ! Optical depth scaling for the cloudy regions
83    real(jprb), intent(out) :: od_scaling(2:nreg,nlev,istartcol:iendcol)
84
85    ! Regions smaller than this are ignored
86    real(jprb), intent(in), optional :: cloud_fraction_threshold
87
88    ! In case the user doesn't supply cloud_fraction_threshold we use
89    ! a default value
90    real(jprb) :: frac_threshold
91
92    ! Loop indices
93    integer :: jcol, jlev
94
95    real(jprb) :: hook_handle
96
97    if (lhook) call dr_hook('radiation_region_properties:calc_region_properties',0,hook_handle)
98
99    if (present(cloud_fraction_threshold)) then
100      frac_threshold = cloud_fraction_threshold
101    else
102      frac_threshold = 1.0e-20_jprb
103    end if
104
105    if (nreg == 2) then
106      ! Only one clear-sky and one cloudy region: cloudy region is
107      ! homogeneous
108      reg_fracs(2,1:nlev,istartcol:iendcol)  = transpose(cloud_fraction(istartcol:iendcol,1:nlev))
109      reg_fracs(1,1:nlev,istartcol:iendcol)  = 1.0_jprb - reg_fracs(2,1:nlev,istartcol:iendcol)
110      od_scaling(2,1:nlev,istartcol:iendcol) = 1.0_jprb
111
112    else if (nreg == 3) then
113      ! Two cloudy regions with optical depth scaled by 1-x and
114      ! 1+x.
115      ! Simple version which fails when fractional_std >= 1:
116      !od_scaling(2) = 1.0_jprb-cloud%fractional_std(jcol,jlev)
117      ! According to Shonk and Hogan (2008), 1-FSD should correspond to
118      ! the 16th percentile.
119      if (.not. do_gamma) then
120        ! If we treat the distribution as a lognormal such that the
121        ! equivalent Normal has a mean mu and standard deviation
122        ! sigma, then the 16th percentile of the lognormal is very
123        ! close to exp(mu-sigma).
124        do jcol = istartcol,iendcol
125          do jlev = 1,nlev
126            if (cloud_fraction(jcol,jlev) < frac_threshold) then
127              reg_fracs(1,jlev,jcol)       = 1.0_jprb
128              reg_fracs(2:3,jlev,jcol)  = 0.0_jprb
129              od_scaling(2:3,jlev,jcol) = 1.0_jprb
130            else
131              reg_fracs(1,jlev,jcol) = 1.0_jprb - cloud_fraction(jcol,jlev)
132              reg_fracs(2:3,jlev,jcol) = cloud_fraction(jcol,jlev)*0.5_jprb
133              od_scaling(2,jlev,jcol) &
134                   &  = exp(-sqrt(log(frac_std(jcol,jlev)**2+1))) &
135                   &  / sqrt(frac_std(jcol,jlev)**2+1)
136              od_scaling(3,jlev,jcol) = 2.0_jprb - od_scaling(2,jlev,jcol)
137            end if
138          end do
139        end do
140      else
141        ! If we treat the distribution as a gamma then the 16th
142        ! percentile is close to the following.  Note that because it
143        ! becomes vanishingly small for FSD >~ 2, we have a lower
144        ! limit of 1/40, and for higher FSDs reduce the fractional
145        ! cover of the denser region - see region_fractions routine
146        ! below
147        do jcol = istartcol,iendcol
148          do jlev = 1,nlev
149            if (cloud_fraction(jcol,jlev) < frac_threshold) then
150              reg_fracs(1,jlev,jcol)    = 1.0_jprb
151              reg_fracs(2:3,jlev,jcol)  = 0.0_jprb
152              od_scaling(2:3,jlev,jcol) = 1.0_jprb
153            else
154              ! Fraction of the clear-sky region
155              reg_fracs(1,jlev,jcol) = 1.0_jprb - cloud_fraction(jcol,jlev)
156!#define OLD_GAMMA_REGION_BEHAVIOUR 1
157#ifdef OLD_GAMMA_REGION_BEHAVIOUR
158              ! Use previous behaviour (ecRad version 1.1.5 and
159              ! earlier): cloudy fractions are the same and there is
160              ! no minimum optical depth scaling; this tends to lead
161              ! to an overprediction of the reflection from scenes
162              ! with a large fractional standard deviation of optical
163              ! depth.
164              ! Fraction and optical-depth scaling of the lower of the
165              ! two cloudy regions
166              reg_fracs(2,jlev,jcol) = cloud_fraction(jcol,jlev) * 0.5_jprb
167              od_scaling(2,jlev,jcol) = &
168                   &  exp(-frac_std(jcol,jlev)*(1.0_jprb + 0.5_jprb*frac_std(jcol,jlev) &
169                   &                     *(1.0_jprb+0.5_jprb*frac_std(jcol,jlev))))
170
171#else
172              ! Improved behaviour.
173              ! Fraction and optical-depth scaling of the lower of the
174              ! two cloudy regions
175              reg_fracs(2,jlev,jcol) = cloud_fraction(jcol,jlev) &
176                   &  * max(MinLowerFrac, min(MaxLowerFrac, &
177                   &  LowerFracFSDIntercept + frac_std(jcol,jlev)*LowerFracFSDGradient))
178              od_scaling(2,jlev,jcol) = MinGammaODScaling &
179                   &  + (1.0_jprb - MinGammaODScaling) &
180                   &    * exp(-frac_std(jcol,jlev)*(1.0_jprb + 0.5_jprb*frac_std(jcol,jlev) &
181                   &                     *(1.0_jprb+0.5_jprb*frac_std(jcol,jlev))))
182#endif
183              ! Fraction of the upper of the two cloudy regions
184              reg_fracs(3,jlev,jcol) = 1.0_jprb-reg_fracs(1,jlev,jcol)-reg_fracs(2,jlev,jcol)
185              ! Ensure conservation of the mean optical depth
186              od_scaling(3,jlev,jcol) = (cloud_fraction(jcol,jlev) &
187                   &  -reg_fracs(2,jlev,jcol)*od_scaling(2,jlev,jcol)) / reg_fracs(3,jlev,jcol)
188            end if
189          end do
190        end do
191      end if
192    else ! nreg > 3
193      write(nulerr,'(a)') '*** Error: only 2 or 3 regions may be specified'
194      call radiation_abort()
195    end if
196
197    if (lhook) call dr_hook('radiation_region_properties:calc_region_properties',1,hook_handle)
198
199  end subroutine calc_region_properties
200
201end module radiation_regions
202
Note: See TracBrowser for help on using the repository browser.