source: LMDZ6/branches/Ocean_skin/libf/phylmd/ecrad/radiation_thermodynamics.F90 @ 5185

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

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

  • Inclusion of the ecrad directory containing the sources of the ECRAD code
    • interface routine : radiation_scheme.F90
  • Adaptation of compilation scripts :
    • compilation under CPP key CPP_ECRAD
    • compilation with option "-rad ecard" or "-ecard true"
    • The "-rad old/rtm/ecran" build option will need to replace the "-rrtm true" and "-ecrad true" options in the future.
  • Runing LMDZ simulations with ecrad, you need :
    • logical key iflag_rrtm = 2 in physiq.def
    • namelist_ecrad (DefLists?)
    • the directory "data" containing the configuration files is temporarily placed in ../libfphylmd/ecrad/
  • Compilation and execution are tested in the 1D case. The repository under svn would allow to continue the implementation work: tests, verification of the results, ...
File size: 11.9 KB
Line 
1! radiation_thermodynamics.F90 - Derived type for pressure & temperature
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-05-11  R. Hogan  Fix startcol/endcol for get_layer_mass
17!   2019-01-14  R. Hogan  Added out_of_physical_bounds routine
18!   2019-01-14  R. Hogan  Capped h2o_sat_liq at 1
19
20module radiation_thermodynamics
21
22  use parkind1, only : jprb
23
24  implicit none
25  public
26
27  !---------------------------------------------------------------------
28  ! Derived type for storing pressure and temperature at half levels
29  type thermodynamics_type
30     real(jprb), allocatable, dimension(:,:) :: &
31          &  pressure_hl, &   ! (ncol,nlev+1) pressure (Pa)
32          &  temperature_hl   ! (ncol,nlev+1) temperature (K)
33
34     ! The following is a function of pressure and temperature: you
35     ! can calculate it according to your favourite formula, or the
36     ! calc_saturation_wrt_liquid subroutine can be used to do this
37     ! approximately
38     real(jprb), allocatable, dimension(:,:) :: &
39          &  h2o_sat_liq ! (ncol,nlev) specific humidity at liquid
40                         ! saturation (kg/kg)
41   contains
42     procedure :: allocate   => allocate_thermodynamics_arrays
43     procedure :: deallocate => deallocate_thermodynamics_arrays
44     procedure :: calc_saturation_wrt_liquid
45     procedure :: get_layer_mass
46     procedure :: get_layer_mass_column
47     procedure :: out_of_physical_bounds
48
49  end type thermodynamics_type
50
51contains
52
53
54  !---------------------------------------------------------------------
55  ! Allocate variables with specified dimensions
56  subroutine allocate_thermodynamics_arrays(this, ncol, nlev, &
57       &                                    use_h2o_sat)
58
59    use yomhook,  only : lhook, dr_hook
60
61    class(thermodynamics_type), intent(inout) :: this
62    integer, intent(in)           :: ncol  ! Number of columns
63    integer, intent(in)           :: nlev  ! Number of levels
64    logical, intent(in), optional :: use_h2o_sat ! Allocate h2o_sat_liq?
65
66    logical :: use_h2o_sat_local
67
68    real(jprb) :: hook_handle
69
70    if (lhook) call dr_hook('radiation_thermodynamics:allocate',0,hook_handle)
71
72    allocate(this%pressure_hl(ncol,nlev+1))
73    allocate(this%temperature_hl(ncol,nlev+1))
74
75    use_h2o_sat_local = .false.
76    if (present(use_h2o_sat)) then
77      use_h2o_sat_local = use_h2o_sat
78    end if
79   
80    if (use_h2o_sat_local) then
81      allocate(this%h2o_sat_liq(ncol,nlev))
82    end if   
83
84    if (lhook) call dr_hook('radiation_thermodynamics:allocate',1,hook_handle)
85
86  end subroutine allocate_thermodynamics_arrays
87
88
89  !---------------------------------------------------------------------
90  ! Deallocate variables
91  subroutine deallocate_thermodynamics_arrays(this)
92
93    use yomhook,  only : lhook, dr_hook
94
95    class(thermodynamics_type), intent(inout) :: this
96
97    real(jprb) :: hook_handle
98
99    if (lhook) call dr_hook('radiation_thermodynamics:deallocate',0,hook_handle)
100
101    if (allocated(this%pressure_hl)) then
102      deallocate(this%pressure_hl)
103    end if
104    if (allocated(this%temperature_hl)) then
105      deallocate(this%temperature_hl)
106    end if
107    if (allocated(this%h2o_sat_liq)) then
108      deallocate(this%h2o_sat_liq)
109    end if
110
111    if (lhook) call dr_hook('radiation_thermodynamics:deallocate',1,hook_handle)
112 
113  end subroutine deallocate_thermodynamics_arrays
114
115
116  !---------------------------------------------------------------------
117  ! Calculate approximate saturation with respect to liquid
118  subroutine calc_saturation_wrt_liquid(this,istartcol,iendcol)
119
120    use yomhook,  only : lhook, dr_hook
121
122    class(thermodynamics_type), intent(inout) :: this
123    integer, intent(in)                       :: istartcol, iendcol
124
125    ! Pressure and temperature at full levels
126    real(jprb) :: pressure, temperature
127
128    ! Vapour pressure (Pa)
129    real(jprb) :: e_sat
130
131    integer :: ncol, nlev ! Dimension sizes
132    integer :: jcol, jlev ! Loop indices for column and level
133
134    real(jprb) :: hook_handle
135
136    if (lhook) call dr_hook('radiation_thermodynamics:calc_saturation_wrt_liquid',0,hook_handle)
137
138    ncol = size(this%pressure_hl,1)
139    nlev = size(this%pressure_hl,2) - 1
140
141    if (.not. allocated(this%h2o_sat_liq)) then
142      allocate(this%h2o_sat_liq(ncol,nlev))
143    end if
144
145    do jlev = 1,nlev
146       do jcol = istartcol,iendcol
147          pressure = 0.5 * (this%pressure_hl(jcol,jlev)+this%pressure_hl(jcol,jlev+1))
148          temperature = 0.5 * (this%temperature_hl(jcol,jlev)+this%temperature_hl(jcol,jlev+1))
149          e_sat = 6.11e2_jprb * exp( 17.269_jprb * (temperature-273.16_jprb) / (temperature-35.86_jprb) )
150          ! This formula can go above 1 at low pressure so needs to be
151          ! capped
152          this%h2o_sat_liq(jcol,jlev) = min(1.0_jprb, 0.622_jprb * e_sat / pressure)
153       end do
154    end do
155
156    if (lhook) call dr_hook('radiation_thermodynamics:calc_saturation_wrt_liquid',1,hook_handle)
157
158  end subroutine calc_saturation_wrt_liquid
159
160
161  !---------------------------------------------------------------------
162  ! Calculate the dry mass of each layer, neglecting humidity effects.
163  ! The first version is for all columns.
164  subroutine get_layer_mass(this,istartcol,iendcol,layer_mass)
165
166    use yomhook,              only : lhook, dr_hook
167    use radiation_constants,  only : AccelDueToGravity
168
169    class(thermodynamics_type), intent(in)  :: this
170    integer,                    intent(in)  :: istartcol, iendcol
171    real(jprb),                 intent(out) :: layer_mass(:,:)
172
173    integer    :: nlev
174    real(jprb) :: inv_g
175
176    real(jprb) :: hook_handle
177
178    if (lhook) call dr_hook('radiation_thermodynamics:get_layer_mass',0,hook_handle)
179
180    nlev  = ubound(this%pressure_hl,2) - 1
181    inv_g = 1.0_jprb / AccelDueToGravity
182
183    layer_mass(istartcol:iendcol,1:nlev) &
184         &  = ( this%pressure_hl(istartcol:iendcol,2:nlev+1) &
185         &     -this%pressure_hl(istartcol:iendcol,1:nlev  )  ) &
186         &  * inv_g
187   
188    if (lhook) call dr_hook('radiation_thermodynamics:get_layer_mass',1,hook_handle)
189
190  end subroutine get_layer_mass
191
192  !---------------------------------------------------------------------
193  ! Calculate the dry mass of each layer, neglecting humidity effects.
194  ! The second version is for one column, the one numbered "icol".
195  subroutine get_layer_mass_column(this, icol, layer_mass)
196
197    use yomhook,              only : lhook, dr_hook
198    use radiation_constants,  only : AccelDueToGravity
199
200    class(thermodynamics_type), intent(in)  :: this
201    integer,                    intent(in)  :: icol
202    real(jprb),                 intent(out) :: layer_mass(:)
203
204    integer    :: nlev
205    real(jprb) :: inv_g
206
207    real(jprb) :: hook_handle
208
209    if (lhook) call dr_hook('radiation_thermodynamics:get_layer_mass_column',0,hook_handle)
210
211    nlev  = ubound(this%pressure_hl,2) - 1
212    inv_g = 1.0_jprb / AccelDueToGravity
213
214    layer_mass = ( this%pressure_hl(icol,2:nlev+1) &
215             &    -this%pressure_hl(icol,1:nlev  )  ) &
216             &   * inv_g
217   
218    if (lhook) call dr_hook('radiation_thermodynamics:get_layer_mass_column',1,hook_handle)
219
220  end subroutine get_layer_mass_column
221
222
223  !---------------------------------------------------------------------
224  ! Estimate the separation between the mid-points of model layers
225  ! given the half-level pressure and temperature.  This is not in
226  ! terms of the "thermodynamics" type as it is useful for computing
227  ! overlap decorrelation lengths and hence cloud cover outside the
228  ! radiation scheme.
229  subroutine get_layer_separation(pressure_hl, temperature_hl, layer_separation)
230
231    use yomhook,              only : lhook, dr_hook
232    use radiation_constants,  only : GasConstantDryAir, AccelDueToGravity
233
234    ! Pressure (Pa) and temperature (K) at half-levels, dimensioned
235    ! (ncol,nlev+1) where ncol is the number of columns and nlev is
236    ! the number of model levels
237    real(jprb), dimension(:,:), intent(in)  :: pressure_hl, temperature_hl
238
239    ! Layer separation in metres, dimensioned (ncol,nlev-1)
240    real(jprb), dimension(:,:), intent(out) :: layer_separation
241
242    ! Ratio of gas constant for dry air to acceleration due to gravity
243    real(jprb), parameter :: R_over_g = GasConstantDryAir / AccelDueToGravity
244
245    ! Loop indices and array bounds
246    integer    :: jlev
247    integer    :: i1, i2, nlev
248
249    real(jprb) :: hook_handle
250
251    if (lhook) call dr_hook('radiation_thermodynamics:get_layer_separation',0,hook_handle)
252
253    i1   = lbound(pressure_hl,1)
254    i2   = ubound(pressure_hl,1)
255    nlev =   size(pressure_hl,2)-1
256
257    if (pressure_hl(i1,2) > pressure_hl(i1,1)) then
258      ! Pressure is increasing with index (order of layers is
259      ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we
260      ! don't take the logarithm of the first pressure in each column.
261      layer_separation(i1:i2,1) = R_over_g * temperature_hl(i1:i2,2) &
262           &                    * log(pressure_hl(i1:i2,3)/pressure_hl(i1:i2,2))
263     
264      ! For other layers we take the separation between midpoints to
265      ! be half the separation between the half-levels at the edge of
266      ! the two adjacent layers
267      do jlev = 2,nlev-1
268        layer_separation(i1:i2,jlev) = (0.5_jprb * R_over_g) * temperature_hl(i1:i2,jlev+1) &
269             &                    * log(pressure_hl(i1:i2,jlev+2)/pressure_hl(i1:i2,jlev))
270
271      end do
272
273    else
274      ! Pressure is decreasing with index (order of layers is surface
275      ! to top-of-atmosphere).  In case pressure_hl(:,nlev+1)=0, we
276      ! don't take the logarithm of the last pressure in each column.
277
278      do jlev = 1,nlev-2
279        layer_separation(i1:i2,jlev) = (0.5_jprb * R_over_g) * temperature_hl(i1:i2,jlev+1) &
280             &                    * log(pressure_hl(i1:i2,jlev)/pressure_hl(i1:i2,jlev+2))
281
282      end do
283      layer_separation(i1:i2,nlev-1) = R_over_g * temperature_hl(i1:i2,nlev) &
284           &                    * log(pressure_hl(i1:i2,nlev-1)/pressure_hl(i1:i2,nlev))
285
286    end if
287
288    if (lhook) call dr_hook('radiation_thermodynamics:get_layer_separation',1,hook_handle)   
289
290  end subroutine get_layer_separation
291
292
293  !---------------------------------------------------------------------
294  ! Return .true. if variables are out of a physically sensible range,
295  ! optionally only considering columns between istartcol and iendcol
296  function out_of_physical_bounds(this, istartcol, iendcol, do_fix) result(is_bad)
297
298    use yomhook,          only : lhook, dr_hook
299    use radiation_config, only : out_of_bounds_2d
300
301    class(thermodynamics_type), intent(inout) :: this
302    integer,           optional,intent(in) :: istartcol, iendcol
303    logical,           optional,intent(in) :: do_fix
304    logical                                :: is_bad
305
306    logical    :: do_fix_local
307
308    real(jprb) :: hook_handle
309
310    if (lhook) call dr_hook('radiation_thermodynamics:out_of_physical_bounds',0,hook_handle)
311
312    if (present(do_fix)) then
313      do_fix_local = do_fix
314    else
315      do_fix_local = .false.
316    end if
317
318    ! Dangerous to cap pressure_hl as then the pressure difference across a layer could be zero
319    is_bad =    out_of_bounds_2d(this%pressure_hl, 'pressure_hl', 0.0_jprb, 110000.0_jprb, &
320         &                       .false., i1=istartcol, i2=iendcol) &
321         & .or. out_of_bounds_2d(this%temperature_hl, 'temperature_hl', 100.0_jprb,  400.0_jprb, &
322         &                       do_fix_local, i1=istartcol, i2=iendcol) &
323         & .or. out_of_bounds_2d(this%h2o_sat_liq, 'h2o_sat_liq', 0.0_jprb, 1.0_jprb, &
324         &                       do_fix_local, i1=istartcol, i2=iendcol)
325
326    if (lhook) call dr_hook('radiation_thermodynamics:out_of_physical_bounds',1,hook_handle)
327
328  end function out_of_physical_bounds
329 
330end module radiation_thermodynamics
Note: See TracBrowser for help on using the repository browser.