SCALE-RM
scale_atmos_profile.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 !-------------------------------------------------------------------------------
11 #include "scalelib.h"
12 module scale_atmos_profile
13  !-----------------------------------------------------------------------------
14  !
15  !++ used modules
16  !
17  use scale_precision
18  use scale_io
19  use scale_prof
20 
21  use scale_const, only: &
22  grav => const_grav, &
23  cpdry => const_cpdry, &
24  rdry => const_rdry, &
25  p00 => const_pre00
26  !-----------------------------------------------------------------------------
27  implicit none
28  private
29  !-----------------------------------------------------------------------------
30  !
31  !++ Public procedure
32  !
33  public :: atmos_profile_isa
34 
35  interface atmos_profile_isa
36  module procedure atmos_profile_isa_1d
37  module procedure atmos_profile_isa_3d
38  end interface atmos_profile_isa
39 
40  !-----------------------------------------------------------------------------
41  !
42  !++ Public parameters & variables
43  !
44  !-----------------------------------------------------------------------------
45  !
46  !++ Private procedure
47  !
48  !-----------------------------------------------------------------------------
49  !
50  !++ Private parameters & variables
51  !
52  ! ISA profile
53  integer, private, parameter :: nref = 8
54  real(RP), private, parameter :: z_isa(nref) = (/ 0.0_rp, &
55  11000.0_rp, &
56  20000.0_rp, &
57  32000.0_rp, &
58  47000.0_rp, &
59  51000.0_rp, &
60  71000.0_rp, &
61  84852.0_rp /)
62  real(RP), private, parameter :: GAMMA(nref) = (/ -6.5e-3_rp, &
63  0.0_rp, &
64  1.0e-3_rp, &
65  2.8e-3_rp, &
66  0.0e-3_rp, &
67  -2.8e-3_rp, &
68  -2.0e-3_rp, &
69  0.0_rp /)
70 
71  !-----------------------------------------------------------------------------
72 contains
73  !-----------------------------------------------------------------------------
75  subroutine atmos_profile_isa_1d( &
76  KA, KS, KE, &
77  temp_sfc, &
78  pres_sfc, &
79  z, &
80  pott )
81  implicit none
82 
83  integer, intent(in) :: KA, KS, KE
84  real(RP), intent(in) :: temp_sfc
85  real(RP), intent(in) :: pres_sfc
86  real(RP), intent(in) :: z (KA)
87  real(RP), intent(out) :: pott(KA)
88 
89  real(RP) :: temp_isa(nref)
90  real(RP) :: pres_isa(nref)
91  real(RP) :: temp(KA)
92  real(RP) :: pres(KA)
93 
94  real(RP) :: gmr ! grav / Rdry
95  real(RP) :: RovCP ! CPdry / Rdry
96  integer :: k, n
97  !---------------------------------------------------------------------------
98 
99  gmr = grav / rdry
100  rovcp = rdry / cpdry
101 
102  !--- ISA profile
103  temp_isa(1) = temp_sfc
104  pres_isa(1) = pres_sfc
105 
106  do n = 2, nref
107  temp_isa(n) = temp_isa(n-1) + gamma(n-1) * ( z_isa(n)-z_isa(n-1) )
108 
109  if ( gamma(n-1) == 0.0_rp ) then
110  pres_isa(n) = pres_isa(n-1) * exp( -gmr / temp_isa(n) * ( z_isa(n)-z_isa(n-1) ) )
111  else
112  pres_isa(n) = pres_isa(n-1) * ( temp_isa(n)/temp_isa(n-1) ) ** ( -gmr/gamma(n-1) )
113  endif
114  enddo
115 
116  log_newline
117  log_info("ATMOS_PROFILE_isa_1D",*) '###### ICAO International Standard Atmosphere ######'
118  log_info_cont(*) ' height: lapse rate: pressure: temperature'
119  do n = 1, nref
120  log_info_cont('(4F13.5)') z_isa(n), gamma(n), pres_isa(n), temp_isa(n)
121  enddo
122  log_info_cont(*) '####################################################'
123 
124  !--- make reference state
125  do k = ks, ke
126  if ( z(k) <= z_isa(1) ) then
127 
128  temp(k) = temp_isa(1) + gamma(1) * ( z(k)-z_isa(1) )
129  pres(k) = pres_isa(1) * ( temp(k)/temp_isa(1) ) ** ( -gmr/gamma(1) )
130 
131  elseif( z(k) > z_isa(nref) ) then
132 
133  temp(k) = temp_isa(nref)
134  pres(k) = pres_isa(nref) * exp( -gmr/temp_isa(nref) * ( z(k)-z_isa(nref) ) )
135 
136  else
137  do n = 2, nref
138  if ( z(k) > z_isa(n-1) .AND. z(k) <= z_isa(n) ) then
139 
140  temp(k) = temp_isa(n-1) + gamma(n-1) * ( z(k)-z_isa(n-1) )
141  if ( gamma(n-1) == 0.0_rp ) then
142  pres(k) = pres_isa(n-1) * exp( -gmr/temp_isa(n-1) * ( z(k)-z_isa(n-1) ) )
143  else
144  pres(k) = pres_isa(n-1) * ( temp(k)/temp_isa(n-1) ) ** ( -gmr/gamma(n-1) )
145  endif
146 
147  endif
148  enddo
149  endif
150 
151  pott(k) = temp(k) * ( p00/pres(k) )**rovcp
152  enddo
153 
154  return
155  end subroutine atmos_profile_isa_1d
156 
157  !-----------------------------------------------------------------------------
159  subroutine atmos_profile_isa_3d( &
160  KA, KS, KE, &
161  IA, IS, IE, &
162  JA, JS, JE, &
163  temp_sfc, &
164  pres_sfc, &
165  z, &
166  pott )
167  implicit none
168 
169  integer, intent(in) :: KA, KS, KE
170  integer, intent(in) :: IA, IS, IE
171  integer, intent(in) :: JA, JS, JE
172  real(RP), intent(in) :: temp_sfc(IA,JA)
173  real(RP), intent(in) :: pres_sfc(IA,JA)
174  real(RP), intent(in) :: z (KA,IA,JA)
175  real(RP), intent(out) :: pott (KA,IA,JA)
176 
177  real(RP) :: temp_isa(nref,IA,JA)
178  real(RP) :: pres_isa(nref,IA,JA)
179  real(RP) :: temp(KA)
180  real(RP) :: pres(KA)
181 
182  real(RP) :: gmr ! grav / Rdry
183  real(RP) :: RovCP ! CPdry / Rdry
184  integer :: k, i, j, n
185  !---------------------------------------------------------------------------
186 
187  gmr = grav / rdry
188  rovcp = rdry / cpdry
189 
190  !--- ISA profile
191  do j = js, je
192  do i = is, ie
193  temp_isa(1,i,j) = temp_sfc(i,j)
194  pres_isa(1,i,j) = pres_sfc(i,j)
195  enddo
196  enddo
197 
198  do j = js, je
199  do i = is, ie
200  do n = 2, nref
201  temp_isa(n,i,j) = temp_isa(n-1,i,j) + gamma(n-1) * ( z_isa(n)-z_isa(n-1) )
202 
203  if ( gamma(n-1) == 0.0_rp ) then
204  pres_isa(n,i,j) = pres_isa(n-1,i,j) * exp( -gmr / temp_isa(n,i,j) * ( z_isa(n)-z_isa(n-1) ) )
205  else
206  pres_isa(n,i,j) = pres_isa(n-1,i,j) * ( temp_isa(n,i,j)/temp_isa(n-1,i,j) ) ** ( -gmr/gamma(n-1) )
207  endif
208  enddo
209  enddo
210  enddo
211 
212  log_newline
213  log_info("ATMOS_PROFILE_isa_3D",*) 'ICAO International Standard Atmosphere'
214  log_info_cont(*) '####################################################'
215  log_info_cont(*) ' height: lapse rate: pressure: temperature'
216  do n = 1, nref
217  log_info_cont('(4F13.5)') z_isa(n), gamma(n), pres_isa(n,is,js), temp_isa(n,is,js)
218  enddo
219  log_info_cont(*) '####################################################'
220 
221  !--- make reference state
222  do j = js, je
223  do i = is, ie
224  do k = ks, ke
225  if ( z(k,i,j) <= z_isa(1) ) then
226 
227  temp(k) = temp_isa(1,i,j) + gamma(1) * ( z(k,i,j)-z_isa(1) )
228  pres(k) = pres_isa(1,i,j) * ( temp(k)/temp_isa(1,i,j) ) ** ( -gmr/gamma(1) )
229 
230  elseif ( z(k,i,j) > z_isa(nref) ) then
231 
232  temp(k) = temp_isa(nref,i,j)
233  pres(k) = pres_isa(nref,i,j) * exp( -gmr/temp_isa(nref,i,j) * ( z(k,i,j)-z_isa(nref) ) )
234 
235  else
236  do n = 2, nref
237  if ( z(k,i,j) > z_isa(n-1) .AND. z(k,i,j) <= z_isa(n) ) then
238 
239  temp(k) = temp_isa(n-1,i,j) + gamma(n-1) * ( z(k,i,j)-z_isa(n-1) )
240  if ( gamma(n-1) == 0.0_rp ) then
241  pres(k) = pres_isa(n-1,i,j) * exp( -gmr/temp_isa(n-1,i,j) * ( z(k,i,j)-z_isa(n-1) ) )
242  else
243  pres(k) = pres_isa(n-1,i,j) * ( temp(k)/temp_isa(n-1,i,j) ) ** ( -gmr/gamma(n-1) )
244  endif
245 
246  endif
247  enddo
248  endif
249 
250  pott(k,i,j) = temp(k) * ( p00/pres(k) )**rovcp
251  enddo
252  enddo
253  enddo
254 
255  return
256  end subroutine atmos_profile_isa_3d
257 
258 end module scale_atmos_profile
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:46
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_profile
module atmosphere / vertical profile
Definition: scale_atmos_profile.F90:12
scale_io
module STDIO
Definition: scale_io.F90:10
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:56
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:55
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:88
Generated by   doxygen 1.8.17