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  !$acc routine vector
82  implicit none
83 
84  integer, intent(in) :: KA, KS, KE
85  real(RP), intent(in) :: temp_sfc
86  real(RP), intent(in) :: pres_sfc
87  real(RP), intent(in) :: z (KA)
88  real(RP), intent(out) :: pott(KA)
89 
90  real(RP) :: temp_isa(nref)
91  real(RP) :: pres_isa(nref)
92  real(RP) :: temp(KA)
93  real(RP) :: pres(KA)
94 
95  real(RP) :: gmr ! grav / Rdry
96  real(RP) :: RovCP ! CPdry / Rdry
97  integer :: k, n
98  !---------------------------------------------------------------------------
99 
100  gmr = grav / rdry
101  rovcp = rdry / cpdry
102 
103  !--- ISA profile
104  temp_isa(1) = temp_sfc
105  pres_isa(1) = pres_sfc
106 
107  !$acc loop seq
108  do n = 2, nref
109  temp_isa(n) = temp_isa(n-1) + gamma(n-1) * ( z_isa(n)-z_isa(n-1) )
110 
111  if ( gamma(n-1) == 0.0_rp ) then
112  pres_isa(n) = pres_isa(n-1) * exp( -gmr / temp_isa(n) * ( z_isa(n)-z_isa(n-1) ) )
113  else
114  pres_isa(n) = pres_isa(n-1) * ( temp_isa(n)/temp_isa(n-1) ) ** ( -gmr/gamma(n-1) )
115  endif
116  enddo
117 
118 #ifndef _OPENACC
119  log_newline
120  log_info("ATMOS_PROFILE_isa_1D",*) '###### ICAO International Standard Atmosphere ######'
121  log_info_cont(*) ' height: lapse rate: pressure: temperature'
122  do n = 1, nref
123  log_info_cont('(4F13.5)') z_isa(n), gamma(n), pres_isa(n), temp_isa(n)
124  enddo
125  log_info_cont(*) '####################################################'
126 #endif
127 
128  !--- make reference state
129  do k = ks, ke
130  if ( z(k) <= z_isa(1) ) then
131 
132  temp(k) = temp_isa(1) + gamma(1) * ( z(k)-z_isa(1) )
133  pres(k) = pres_isa(1) * ( temp(k)/temp_isa(1) ) ** ( -gmr/gamma(1) )
134 
135  elseif( z(k) > z_isa(nref) ) then
136 
137  temp(k) = temp_isa(nref)
138  pres(k) = pres_isa(nref) * exp( -gmr/temp_isa(nref) * ( z(k)-z_isa(nref) ) )
139 
140  else
141  do n = 2, nref
142  if ( z(k) > z_isa(n-1) .AND. z(k) <= z_isa(n) ) then
143 
144  temp(k) = temp_isa(n-1) + gamma(n-1) * ( z(k)-z_isa(n-1) )
145  if ( gamma(n-1) == 0.0_rp ) then
146  pres(k) = pres_isa(n-1) * exp( -gmr/temp_isa(n-1) * ( z(k)-z_isa(n-1) ) )
147  else
148  pres(k) = pres_isa(n-1) * ( temp(k)/temp_isa(n-1) ) ** ( -gmr/gamma(n-1) )
149  endif
150 
151  endif
152  enddo
153  endif
154 
155  pott(k) = temp(k) * ( p00/pres(k) )**rovcp
156  enddo
157 
158  return
159  end subroutine atmos_profile_isa_1d
160 
161  !-----------------------------------------------------------------------------
163  subroutine atmos_profile_isa_3d( &
164  KA, KS, KE, &
165  IA, IS, IE, &
166  JA, JS, JE, &
167  temp_sfc, &
168  pres_sfc, &
169  z, &
170  pott )
171  implicit none
172 
173  integer, intent(in) :: KA, KS, KE
174  integer, intent(in) :: IA, IS, IE
175  integer, intent(in) :: JA, JS, JE
176  real(RP), intent(in) :: temp_sfc(IA,JA)
177  real(RP), intent(in) :: pres_sfc(IA,JA)
178  real(RP), intent(in) :: z (KA,IA,JA)
179  real(RP), intent(out) :: pott (KA,IA,JA)
180 
181  real(RP) :: temp_isa(nref,IA,JA)
182  real(RP) :: pres_isa(nref,IA,JA)
183  real(RP) :: temp(KA)
184  real(RP) :: pres(KA)
185 
186  real(RP) :: gmr ! grav / Rdry
187  real(RP) :: RovCP ! CPdry / Rdry
188  integer :: k, i, j, n
189  !---------------------------------------------------------------------------
190 
191  !$acc data copyin(temp_sfc, pres_sfc, z) copyout(pott) create(temp_isa, pres_isa)
192 
193  gmr = grav / rdry
194  rovcp = rdry / cpdry
195 
196  !--- ISA profile
197  !$acc kernels
198  do j = js, je
199  do i = is, ie
200  temp_isa(1,i,j) = temp_sfc(i,j)
201  pres_isa(1,i,j) = pres_sfc(i,j)
202  enddo
203  enddo
204  !$acc end kernels
205 
206  !$acc kernels
207  do j = js, je
208  do i = is, ie
209  do n = 2, nref
210  temp_isa(n,i,j) = temp_isa(n-1,i,j) + gamma(n-1) * ( z_isa(n)-z_isa(n-1) )
211 
212  if ( gamma(n-1) == 0.0_rp ) then
213  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) ) )
214  else
215  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) )
216  endif
217  enddo
218  enddo
219  enddo
220  !$acc end kernels
221 
222  log_newline
223  log_info("ATMOS_PROFILE_isa_3D",*) 'ICAO International Standard Atmosphere'
224  log_info_cont(*) '####################################################'
225  log_info_cont(*) ' height: lapse rate: pressure: temperature'
226  do n = 1, nref
227  log_info_cont('(4F13.5)') z_isa(n), gamma(n), pres_isa(n,is,js), temp_isa(n,is,js)
228  enddo
229  log_info_cont(*) '####################################################'
230 
231  !--- make reference state
232  !$acc kernels
233  !$acc loop private(temp, pres)
234  do j = js, je
235  !$acc loop private(temp, pres)
236  do i = is, ie
237  do k = ks, ke
238  if ( z(k,i,j) <= z_isa(1) ) then
239 
240  temp(k) = temp_isa(1,i,j) + gamma(1) * ( z(k,i,j)-z_isa(1) )
241  pres(k) = pres_isa(1,i,j) * ( temp(k)/temp_isa(1,i,j) ) ** ( -gmr/gamma(1) )
242 
243  elseif ( z(k,i,j) > z_isa(nref) ) then
244 
245  temp(k) = temp_isa(nref,i,j)
246  pres(k) = pres_isa(nref,i,j) * exp( -gmr/temp_isa(nref,i,j) * ( z(k,i,j)-z_isa(nref) ) )
247 
248  else
249  do n = 2, nref
250  if ( z(k,i,j) > z_isa(n-1) .AND. z(k,i,j) <= z_isa(n) ) then
251 
252  temp(k) = temp_isa(n-1,i,j) + gamma(n-1) * ( z(k,i,j)-z_isa(n-1) )
253  if ( gamma(n-1) == 0.0_rp ) then
254  pres(k) = pres_isa(n-1,i,j) * exp( -gmr/temp_isa(n-1,i,j) * ( z(k,i,j)-z_isa(n-1) ) )
255  else
256  pres(k) = pres_isa(n-1,i,j) * ( temp(k)/temp_isa(n-1,i,j) ) ** ( -gmr/gamma(n-1) )
257  endif
258 
259  endif
260  enddo
261  endif
262 
263  pott(k,i,j) = temp(k) * ( p00/pres(k) )**rovcp
264  enddo
265  enddo
266  enddo
267  !$acc end kernels
268 
269  !$acc end data
270 
271  return
272  end subroutine atmos_profile_isa_3d
273 
274 end module scale_atmos_profile
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:49
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:60
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:59
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:97
Generated by   doxygen 1.8.17