SCALE-RM
scale_interp_vert.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 #include "scalelib.h"
12  !-----------------------------------------------------------------------------
13  !
14  !++ used modules
15  !
16  use scale_precision
17  use scale_io
18  use scale_prof
20  !-----------------------------------------------------------------------------
21  implicit none
22  private
23  !-----------------------------------------------------------------------------
24  !
25  !++ Public procedure
26  !
27  public :: interp_vert_setcoef
28  public :: interp_vert_xi2z
29  public :: interp_vert_xih2zh
30  public :: interp_vert_z2xi
31  public :: interp_vert_zh2xih
32 
33  public :: interp_vert_alloc_pres
34  public :: interp_vert_setcoef_pres
35  public :: interp_vert_xi2p
36  public :: interp_vert_xih2p
37 
38  !-----------------------------------------------------------------------------
39  !
40  !++ Public parameters & variables
41  !
42  logical, public :: interp_available = .false.
43 
44  !-----------------------------------------------------------------------------
45  !
46  !++ Private procedure
47  !
48  !-----------------------------------------------------------------------------
49  !
50  !++ Private parameters & variables
51  !
52  ! full level
53  integer, private, allocatable :: interp_xi2z_idx (:,:,:,:)
54  real(RP), private, allocatable :: interp_xi2z_coef (:,:,:,:)
55  integer, private, allocatable :: interp_z2xi_idx (:,:,:,:)
56  real(RP), private, allocatable :: interp_z2xi_coef (:,:,:,:)
57 
58  integer, private, allocatable :: interp_xi2p_idx (:,:,:,:)
59  real(RP), private, allocatable :: interp_xi2p_coef (:,:,:,:)
60 
61  ! half level
62  integer, private, allocatable :: interp_xih2zh_idx (:,:,:,:)
63  real(RP), private, allocatable :: interp_xih2zh_coef(:,:,:,:)
64  integer, private, allocatable :: interp_zh2xih_idx (:,:,:,:)
65  real(RP), private, allocatable :: interp_zh2xih_coef(:,:,:,:)
66 
67  integer, private, allocatable :: interp_xih2p_idx (:,:,:,:)
68  real(RP), private, allocatable :: interp_xih2p_coef (:,:,:,:)
69 
70  !-----------------------------------------------------------------------------
71 contains
72  !-----------------------------------------------------------------------------
74  subroutine interp_vert_setcoef( &
75  KA, KS, KE, &
76  IA, IS, IE, &
77  JA, JS, JE, &
78  TOPO_exist, &
79  Xi, Xih, &
80  Z, Zh )
81  implicit none
82 
83  integer, intent(in) :: KA, KS, KE
84  integer, intent(in) :: IA, IS, IE
85  integer, intent(in) :: JA, JS, JE
86  logical, intent(in) :: TOPO_exist
87  real(RP), intent(in) :: Xi ( ka)
88  real(RP), intent(in) :: Xih(0:ka)
89  real(RP), intent(in) :: Z ( ka,ia,ja)
90  real(RP), intent(in) :: Zh (0:ka,ia,ja)
91 
92  integer :: k, i, j, kk, kp
93  !---------------------------------------------------------------------------
94 
95  log_newline
96  log_info("INTERP_VERT_setcoef",*) 'Setup'
97  log_info("INTERP_VERT_setcoef",*) 'No namelists.'
98 
99  interp_available = topo_exist
100 
101  log_newline
102  log_info("INTERP_VERT_setcoef",*) 'Topography exists & interpolation has meaning? : ', interp_available
103 
104  ! full level
105 
106  allocate( interp_xi2z_idx(ka,ia,ja,2) )
107  allocate( interp_xi2z_coef(ka,ia,ja,3) )
108  allocate( interp_z2xi_idx(ka,ia,ja,2) )
109  allocate( interp_z2xi_coef(ka,ia,ja,3) )
110 
111  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
112  !$omp private(i,j,k,kk,kp) &
113  !$omp shared(JA,IA,KS,KE,Xi,Zh,Z,INTERP_xi2z_idx,INTERP_xi2z_coef)
114  do j = 1, ja
115  do i = 1, ia
116  do k = ks, ke
117  if ( xi(k) <= zh(ks-1,i,j) ) then
118 
119  interp_xi2z_idx(k,i,j,1) = ks ! dummmy
120  interp_xi2z_idx(k,i,j,2) = ks ! dummmy
121  interp_xi2z_coef(k,i,j,1) = 0.0_rp
122  interp_xi2z_coef(k,i,j,2) = 0.0_rp
123  interp_xi2z_coef(k,i,j,3) = 1.0_rp ! set UNDEF
124 
125  elseif( xi(k) <= z(ks,i,j) ) then
126 
127  interp_xi2z_idx(k,i,j,1) = ks ! dummmy
128  interp_xi2z_idx(k,i,j,2) = ks
129  interp_xi2z_coef(k,i,j,1) = 0.0_rp
130  interp_xi2z_coef(k,i,j,2) = 1.0_rp
131  interp_xi2z_coef(k,i,j,3) = 0.0_rp
132 
133  elseif( xi(k) > zh(ke,i,j) ) then
134 
135  interp_xi2z_idx(k,i,j,1) = ke ! dummmy
136  interp_xi2z_idx(k,i,j,2) = ke ! dummmy
137  interp_xi2z_coef(k,i,j,1) = 0.0_rp
138  interp_xi2z_coef(k,i,j,2) = 0.0_rp
139  interp_xi2z_coef(k,i,j,3) = 1.0_rp ! set UNDEF
140 
141  elseif( xi(k) > z(ke,i,j) ) then
142 
143  interp_xi2z_idx(k,i,j,1) = ke
144  interp_xi2z_idx(k,i,j,2) = ke ! dummmy
145  interp_xi2z_coef(k,i,j,1) = 1.0_rp
146  interp_xi2z_coef(k,i,j,2) = 0.0_rp
147  interp_xi2z_coef(k,i,j,3) = 0.0_rp
148 
149  else
150 
151  do kk = ks+1, ke
152  kp = kk
153  if( xi(k) <= z(kk,i,j) ) exit
154  enddo
155 
156  interp_xi2z_idx(k,i,j,1) = kp - 1
157  interp_xi2z_idx(k,i,j,2) = kp
158  interp_xi2z_coef(k,i,j,1) = ( z(kp,i,j) - xi(k) ) &
159  / ( z(kp,i,j) - z(kp-1,i,j) )
160  interp_xi2z_coef(k,i,j,2) = ( xi(k) - z(kp-1,i,j) ) &
161  / ( z(kp,i,j) - z(kp-1,i,j) )
162  interp_xi2z_coef(k,i,j,3) = 0.0_rp
163 
164  endif
165  enddo
166  enddo
167  enddo
168 
169  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
170  !$omp private(i,j,k,kk,kp) &
171  !$omp shared(JA,IA,KS,KE,Z,Xih,Xi,INTERP_z2xi_idx,INTERP_z2xi_coef)
172  do j = 1, ja
173  do i = 1, ia
174  do k = ks, ke
175  if ( z(k,i,j) < xih(ks-1) ) then
176 
177  interp_z2xi_idx(k,i,j,1) = ks ! dummmy
178  interp_z2xi_idx(k,i,j,2) = ks
179  interp_z2xi_coef(k,i,j,1) = 0.0_rp
180  interp_z2xi_coef(k,i,j,2) = 1.0_rp
181 ! INTERP_z2xi_coef(k,i,j,3) = 1.0_RP ! set UNDEF
182 
183  elseif( z(k,i,j) <= xi(ks) ) then
184 
185  interp_z2xi_idx(k,i,j,1) = ks ! dummmy
186  interp_z2xi_idx(k,i,j,2) = ks
187  interp_z2xi_coef(k,i,j,1) = 0.0_rp
188  interp_z2xi_coef(k,i,j,2) = 1.0_rp
189  interp_z2xi_coef(k,i,j,3) = 0.0_rp
190 
191  elseif( z(k,i,j) > xih(ke) ) then
192 
193  interp_z2xi_idx(k,i,j,1) = ke ! dummmy
194  interp_z2xi_idx(k,i,j,2) = ke ! dummmy
195  interp_z2xi_coef(k,i,j,1) = 0.0_rp
196  interp_z2xi_coef(k,i,j,2) = 0.0_rp
197  interp_z2xi_coef(k,i,j,3) = 1.0_rp ! set UNDEF
198 
199  elseif( z(k,i,j) > xi(ke) ) then
200 
201  interp_z2xi_idx(k,i,j,1) = ke
202  interp_z2xi_idx(k,i,j,2) = ke ! dummmy
203  interp_z2xi_coef(k,i,j,1) = 1.0_rp
204  interp_z2xi_coef(k,i,j,2) = 0.0_rp
205  interp_z2xi_coef(k,i,j,3) = 0.0_rp
206 
207  else
208 
209  do kk = ks+1, ke
210  kp = kk
211  if( z(k,i,j) <= xi(kk) ) exit
212  enddo
213 
214  interp_z2xi_idx(k,i,j,1) = kp - 1
215  interp_z2xi_idx(k,i,j,2) = kp
216  interp_z2xi_coef(k,i,j,1) = ( xi(kp) - z(k,i,j) ) &
217  / ( xi(kp) - xi(kp-1) )
218  interp_z2xi_coef(k,i,j,2) = ( z(k,i,j) - xi(kp-1) ) &
219  / ( xi(kp) - xi(kp-1) )
220  interp_z2xi_coef(k,i,j,3) = 0.0_rp
221 
222  endif
223  enddo
224  enddo
225  enddo
226 
227  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
228  !$omp private(i,j) &
229  !$omp shared(JA,IA,KS,KE,KA,INTERP_xi2z_idx,INTERP_xi2z_coef,INTERP_z2xi_idx,INTERP_z2xi_coef)
230  do j = 1, ja
231  do i = 1, ia
232  interp_xi2z_idx( 1:ks-1,i,j,1) = ks ! dummmy
233  interp_xi2z_idx( 1:ks-1,i,j,2) = ks ! dummmy
234  interp_xi2z_coef( 1:ks-1,i,j,1) = 0.0_rp
235  interp_xi2z_coef( 1:ks-1,i,j,2) = 0.0_rp
236  interp_xi2z_coef( 1:ks-1,i,j,3) = 1.0_rp ! set UNDEF
237 
238  interp_xi2z_idx(ke+1:ka,i,j,1) = ke ! dummmy
239  interp_xi2z_idx(ke+1:ka,i,j,2) = ke ! dummmy
240  interp_xi2z_coef(ke+1:ka,i,j,1) = 0.0_rp
241  interp_xi2z_coef(ke+1:ka,i,j,2) = 0.0_rp
242  interp_xi2z_coef(ke+1:ka,i,j,3) = 1.0_rp ! set UNDEF
243 
244  interp_z2xi_idx( 1:ks-1,i,j,1) = ks ! dummmy
245  interp_z2xi_idx( 1:ks-1,i,j,2) = ks ! dummmy
246  interp_z2xi_coef( 1:ks-1,i,j,1) = 0.0_rp
247  interp_z2xi_coef( 1:ks-1,i,j,2) = 0.0_rp
248  interp_z2xi_coef( 1:ks-1,i,j,3) = 1.0_rp ! set UNDEF
249 
250  interp_z2xi_idx(ke+1:ka,i,j,1) = ke ! dummmy
251  interp_z2xi_idx(ke+1:ka,i,j,2) = ke ! dummmy
252  interp_z2xi_coef(ke+1:ka,i,j,1) = 0.0_rp
253  interp_z2xi_coef(ke+1:ka,i,j,2) = 0.0_rp
254  interp_z2xi_coef(ke+1:ka,i,j,3) = 1.0_rp ! set UNDEF
255  enddo
256  enddo
257 
258 
259  ! half level
260 
261  allocate( interp_xih2zh_idx(0:ka,ia,ja,2) )
262  allocate( interp_xih2zh_coef(0:ka,ia,ja,3) )
263  allocate( interp_zh2xih_idx(0:ka,ia,ja,2) )
264  allocate( interp_zh2xih_coef(0:ka,ia,ja,3) )
265 
266  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
267  !$omp private(i,j,k,kk,kp) &
268  !$omp shared(JA,IA,KS,KE,Xih,Zh,INTERP_xih2zh_idx,INTERP_xih2zh_coef)
269  do j = 1, ja
270  do i = 1, ia
271  do k = ks-1, ke
272  if ( xih(k) < zh(ks-1,i,j) ) then
273 
274  interp_xih2zh_idx(k,i,j,1) = ks ! dummmy
275  interp_xih2zh_idx(k,i,j,2) = ks ! dummmy
276  interp_xih2zh_coef(k,i,j,1) = 0.0_rp
277  interp_xih2zh_coef(k,i,j,2) = 0.0_rp
278  interp_xih2zh_coef(k,i,j,3) = 1.0_rp ! set UNDEF
279 
280  elseif( xih(k) > zh(ke,i,j) ) then
281 
282  interp_xih2zh_idx(k,i,j,1) = ke ! dummmy
283  interp_xih2zh_idx(k,i,j,2) = ke ! dummmy
284  interp_xih2zh_coef(k,i,j,1) = 0.0_rp
285  interp_xih2zh_coef(k,i,j,2) = 0.0_rp
286  interp_xih2zh_coef(k,i,j,3) = 1.0_rp ! set UNDEF
287 
288  else
289 
290  do kk = ks, ke
291  kp = kk
292  if( xih(k) <= zh(kk,i,j) ) exit
293  enddo
294 
295  interp_xih2zh_idx(k,i,j,1) = kp - 1
296  interp_xih2zh_idx(k,i,j,2) = kp
297  interp_xih2zh_coef(k,i,j,1) = ( zh(kp,i,j) - xih(k) ) &
298  / ( zh(kp,i,j) - zh(kp-1,i,j) )
299  interp_xih2zh_coef(k,i,j,2) = ( xih(k) - zh(kp-1,i,j) ) &
300  / ( zh(kp,i,j) - zh(kp-1,i,j) )
301  interp_xih2zh_coef(k,i,j,3) = 0.0_rp
302 
303  endif
304  enddo
305  enddo
306  enddo
307 
308  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
309  !$omp private(i,j,k,kk,kp) &
310  !$omp shared(JA,IA,KS,KE,Xih,Zh,INTERP_zh2xih_idx,INTERP_zh2xih_coef)
311  do j = 1, ja
312  do i = 1, ia
313  do k = ks-1, ke
314  if ( zh(k,i,j) <= xih(ks-1) ) then
315 
316  interp_zh2xih_idx(k,i,j,1) = ks ! dummmy
317  interp_zh2xih_idx(k,i,j,2) = ks ! dummmy
318  interp_zh2xih_coef(k,i,j,1) = 0.0_rp
319  interp_zh2xih_coef(k,i,j,2) = 0.0_rp
320  interp_zh2xih_coef(k,i,j,3) = 1.0_rp ! set UNDEF
321 
322  elseif( zh(k,i,j) > xih(ke) ) then
323 
324  interp_zh2xih_idx(k,i,j,1) = ke ! dummmy
325  interp_zh2xih_idx(k,i,j,2) = ke ! dummmy
326  interp_zh2xih_coef(k,i,j,1) = 0.0_rp
327  interp_zh2xih_coef(k,i,j,2) = 0.0_rp
328  interp_zh2xih_coef(k,i,j,3) = 1.0_rp ! set UNDEF
329 
330  else
331 
332  do kk = ks, ke
333  kp = kk
334  if( zh(k,i,j) <= xih(kk) ) exit
335  enddo
336 
337  interp_zh2xih_idx(k,i,j,1) = kp - 1
338  interp_zh2xih_idx(k,i,j,2) = kp
339  interp_zh2xih_coef(k,i,j,1) = ( xih(kp) - zh(k,i,j) ) &
340  / ( xih(kp) - xih(kp-1) )
341  interp_zh2xih_coef(k,i,j,2) = ( zh(k,i,j) - xih(kp-1) ) &
342  / ( xih(kp) - xih(kp-1) )
343  interp_zh2xih_coef(k,i,j,3) = 0.0_rp
344 
345  endif
346  enddo
347  enddo
348  enddo
349 
350  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
351  !$omp private(i,j) &
352  !$omp shared(JA,IA,KS,KE,KA,INTERP_xih2zh_idx,INTERP_xih2zh_coef,INTERP_zh2xih_idx,INTERP_zh2xih_coef)
353  do j = 1, ja
354  do i = 1, ia
355  interp_xih2zh_idx( 1:ks-2,i,j,1) = ks ! dummmy
356  interp_xih2zh_idx( 1:ks-2,i,j,2) = ks ! dummmy
357  interp_xih2zh_coef( 1:ks-2,i,j,1) = 0.0_rp
358  interp_xih2zh_coef( 1:ks-2,i,j,2) = 0.0_rp
359  interp_xih2zh_coef( 1:ks-2,i,j,3) = 1.0_rp ! set UNDEF
360 
361  interp_xih2zh_idx(ke+1:ka,i,j,1) = ke ! dummmy
362  interp_xih2zh_idx(ke+1:ka,i,j,2) = ke ! dummmy
363  interp_xih2zh_coef(ke+1:ka,i,j,1) = 0.0_rp
364  interp_xih2zh_coef(ke+1:ka,i,j,2) = 0.0_rp
365  interp_xih2zh_coef(ke+1:ka,i,j,3) = 1.0_rp ! set UNDEF
366 
367  interp_zh2xih_idx( 1:ks-2,i,j,1) = ks ! dummmy
368  interp_zh2xih_idx( 1:ks-2,i,j,2) = ks ! dummmy
369  interp_zh2xih_coef( 1:ks-2,i,j,1) = 0.0_rp
370  interp_zh2xih_coef( 1:ks-2,i,j,2) = 0.0_rp
371  interp_zh2xih_coef( 1:ks-2,i,j,3) = 1.0_rp ! set UNDEF
372 
373  interp_zh2xih_idx(ke+1:ka,i,j,1) = ke ! dummmy
374  interp_zh2xih_idx(ke+1:ka,i,j,2) = ke ! dummmy
375  interp_zh2xih_coef(ke+1:ka,i,j,1) = 0.0_rp
376  interp_zh2xih_coef(ke+1:ka,i,j,2) = 0.0_rp
377  interp_zh2xih_coef(ke+1:ka,i,j,3) = 1.0_rp ! set UNDEF
378  enddo
379  enddo
380 
381  return
382  end subroutine interp_vert_setcoef
383 
384  !-----------------------------------------------------------------------------
385  subroutine interp_vert_xi2z( &
386  KA, KS, KE, &
387  IA, IS, IE, &
388  JA, JS, JE, &
389  Xi, &
390  Z, &
391  var, &
392  var_Z )
393  use scale_const, only: &
395  use scale_matrix, only: &
396  matrix_solver_tridiagonal
397  implicit none
398 
399  integer, intent(in) :: KA, KS, KE
400  integer, intent(in) :: IA, IS, IE
401  integer, intent(in) :: JA, JS, JE
402  real(RP), intent(in) :: Xi (ka)
403  real(RP), intent(in) :: Z (ka,ia,ja)
404  real(RP), intent(in) :: var (ka,ia,ja)
405  real(RP), intent(out) :: var_Z(ka,ia,ja)
406 
407  real(RP) :: FDZ(ka)
408  real(RP) :: MD (ka)
409  real(RP) :: U (ka)
410  real(RP) :: V (ka)
411  real(RP) :: c1, c2, c3, d
412  integer :: kmax
413 
414  integer :: k, i, j, kk
415  !---------------------------------------------------------------------------
416 
417  kmax = ke - ks + 1
418 
419  if ( kmax == 2 ) then
420 
421  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
422  !$omp private(k,i,j) &
423  !$omp shared(KA,IS,IE,JS,JE) &
424  !$omp shared(var,var_Z) &
425  !$omp shared(INTERP_xi2z_idx,INTERP_xi2z_coef,CONST_UNDEF)
426  do j = js, je
427  do i = is, ie
428  do k = 1, ka
429  var_z(k,i,j) = interp_xi2z_coef(k,i,j,1) * var(interp_xi2z_idx(k,i,j,1),i,j) &
430  + interp_xi2z_coef(k,i,j,2) * var(interp_xi2z_idx(k,i,j,2),i,j) &
431  + interp_xi2z_coef(k,i,j,3) * const_undef
432  enddo
433  enddo
434  enddo
435 
436  else ! cubic spline
437 
438  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
439  !$omp private(k,i,j) &
440  !$omp private(FDZ,MD,U,V,kk,c1,c2,c3,d) &
441  !$omp shared(KA,KS,KE,kmax,IS,IE,JS,JE) &
442  !$omp shared(var,var_Z) &
443  !$omp shared(INTERP_xi2z_idx,INTERP_xi2z_coef,Xi,Z,CONST_UNDEF)
444  do j = js, je
445  do i = is, ie
446 
447  do k = ks, ke-1
448  fdz(k) = z(k+1,i,j) - z(k,i,j)
449  end do
450 
451  md(ks+1) = 2.0 * ( fdz(ks) + fdz(ks+1) ) + fdz(ks)
452  do k = ks+2, ke-2
453  md(k) = 2.0 * ( fdz(k-1) + fdz(k) )
454  end do
455  md(ke-1) = 2.0 * ( fdz(ke-2) + fdz(ke-1) ) + fdz(ke-1)
456 
457  do k = ks+1, ke-1
458  v(k) = ( var(k+1,i,j) - var(k ,i,j) ) / fdz(k ) &
459  - ( var(k ,i,j) - var(k-1,i,j) ) / fdz(k-1)
460  end do
461 
462  call matrix_solver_tridiagonal( ka, ks+1, ke-1, & ! [IN]
463  fdz(:), & ! [IN]
464  md(:), & ! [IN]
465  fdz(:), & ! [IN]
466  v(:), & ! [IN]
467  u(:) ) ! [OUT]
468 
469  u(ks) = u(ks+1)
470  u(ke) = u(ke-1)
471 
472  do k = 1, ka
473  kk = min(interp_xi2z_idx(k,i,j,1),ke-1)
474 
475  c3 = ( u(kk+1) - u(kk) ) / fdz(kk)
476  c2 = 3.0_rp * u(kk)
477  c1 = ( var(kk+1,i,j) - var(kk,i,j) ) / fdz(kk) - ( u(kk) * 2.0_rp + u(kk+1) ) * fdz(kk)
478  d = xi(k) - z(kk,i,j)
479 
480  var_z(k,i,j) = ( interp_xi2z_coef(k,i,j,3) ) * const_undef &
481  + ( 1.0_rp-interp_xi2z_coef(k,i,j,3) ) * ( ( ( c3 * d + c2 ) * d + c1 ) * d + var(kk,i,j) )
482  end do
483 
484  end do
485  end do
486 
487  end if
488 
489  return
490  end subroutine interp_vert_xi2z
491 
492  !-----------------------------------------------------------------------------
493  subroutine interp_vert_z2xi( &
494  KA, KS, KE, &
495  IA, IS, IE, &
496  JA, JS, JE, &
497  Z, &
498  Xi, &
499  var, &
500  var_Xi )
501  use scale_const, only: &
503  use scale_matrix, only: &
504  matrix_solver_tridiagonal
505  implicit none
506 
507  integer, intent(in) :: KA, KS, KE
508  integer, intent(in) :: IA, IS, IE
509  integer, intent(in) :: JA, JS, JE
510  real(RP), intent(in) :: Z (ka,ia,ja)
511  real(RP), intent(in) :: Xi (ka)
512  real(RP), intent(in) :: var (ka,ia,ja)
513  real(RP), intent(out) :: var_Xi(ka,ia,ja)
514 
515  real(RP) :: FDZ(ka)
516  real(RP) :: MD (ka)
517  real(RP) :: U (ka)
518  real(RP) :: V (ka)
519  real(RP) :: c1, c2, c3, d
520  integer :: kmax
521 
522  integer :: k, i, j, kk
523  !---------------------------------------------------------------------------
524 
525  kmax = ke - ks + 1
526 
527  if ( kmax <= 2 ) then
528 
529  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
530  !$omp private(k,i,j) &
531  !$omp shared(KA,IS,IE,JS,JE) &
532  !$omp shared(var,var_Xi) &
533  !$omp shared(INTERP_z2xi_idx,INTERP_z2xi_coef,CONST_UNDEF)
534  do j = js, je
535  do i = is, ie
536  do k = 1, ka
537  var_xi(k,i,j) = interp_z2xi_coef(k,i,j,1) * var(interp_z2xi_idx(k,i,j,1),i,j) &
538  + interp_z2xi_coef(k,i,j,2) * var(interp_z2xi_idx(k,i,j,2),i,j) &
539  + interp_z2xi_coef(k,i,j,3) * const_undef
540  enddo
541  enddo
542  enddo
543 
544  else ! cubic spline
545 
546  do k = ks, ke-1
547  fdz(k) = xi(k+1) - xi(k)
548  end do
549 
550  md(ks+1) = 2.0 * ( fdz(ks) + fdz(ks+1) ) + fdz(ks)
551  do k = ks+2, ke-2
552  md(k) = 2.0 * ( fdz(k-1) + fdz(k) )
553  end do
554  md(ke-1) = 2.0 * ( fdz(ke-2) + fdz(ke-1) ) + fdz(ke-1)
555 
556  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
557  !$omp private(k,i,j) &
558  !$omp private(U,V,kk,c1,c2,c3,d) &
559  !$omp shared(KA,KS,KE,kmax,IS,IE,JS,JE) &
560  !$omp shared(var,var_Xi,MD,FDZ) &
561  !$omp shared(INTERP_z2xi_idx,INTERP_z2xi_coef,Xi,Z,CONST_UNDEF)
562  do j = js, je
563  do i = is, ie
564 
565  do k = ks+1, ke-1
566  v(k) = ( var(k+1,i,j) - var(k ,i,j) ) / fdz(k ) &
567  - ( var(k ,i,j) - var(k-1,i,j) ) / fdz(k-1)
568  end do
569 
570  call matrix_solver_tridiagonal( ka, ks+1, ke-1, & ! [IN]
571  fdz(:), & ! [IN]
572  md(:), & ! [IN]
573  fdz(:), & ! [IN]
574  v(:), & ! [IN]
575  u(:) ) ! [OUT]
576 
577  u(ks) = u(ks+1)
578  u(ke) = u(ke-1)
579 
580  do k = 1, ka
581  kk = min(interp_z2xi_idx(k,i,j,1),ke-1)
582 
583  c3 = ( u(kk+1) - u(kk) ) / fdz(kk)
584  c2 = 3.0_rp * u(kk)
585  c1 = ( var(kk+1,i,j) - var(kk,i,j) ) / fdz(kk) - ( u(kk) * 2.0_rp + u(kk+1) ) * fdz(kk)
586  d = z(k,i,j) - xi(kk)
587 
588  var_xi(k,i,j) = ( interp_z2xi_coef(k,i,j,3) ) * const_undef &
589  + ( 1.0_rp-interp_z2xi_coef(k,i,j,3) ) * ( ( ( c3 * d + c2 ) * d + c1 ) * d + var(kk,i,j) )
590  end do
591 
592  end do
593  end do
594 
595  end if
596 
597  return
598  end subroutine interp_vert_z2xi
599 
600  !-----------------------------------------------------------------------------
601  subroutine interp_vert_xih2zh( &
602  KA, KS, KE, &
603  IA, IS, IE, &
604  JA, JS, JE, &
605  Xih, &
606  Zh, &
607  var, &
608  var_Z )
609  use scale_const, only: &
611  use scale_matrix, only: &
612  matrix_solver_tridiagonal
613  implicit none
614 
615  integer, intent(in) :: KA, KS, KE
616  integer, intent(in) :: IA, IS, IE
617  integer, intent(in) :: JA, JS, JE
618  real(RP), intent(in) :: Xih (0:ka)
619  real(RP), intent(in) :: Zh (0:ka,ia,ja)
620  real(RP), intent(in) :: var (ka,ia,ja)
621  real(RP), intent(out) :: var_Z(ka,ia,ja)
622 
623  real(RP) :: CDZ(ka)
624  real(RP) :: MD (ka)
625  real(RP) :: U (ka)
626  real(RP) :: V (ka)
627  real(RP) :: c1, c2, c3, d
628  integer :: kmax
629 
630  integer :: k, i, j, kk
631  !---------------------------------------------------------------------------
632 
633  kmax = ke - ks + 1
634 
635  if ( kmax == 2 ) then
636 
637  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
638  !$omp private(k,i,j) &
639  !$omp shared(KA,IS,IE,JS,JE) &
640  !$omp shared(var,var_Z) &
641  !$omp shared(INTERP_xih2zh_idx,INTERP_xih2zh_coef,CONST_UNDEF)
642  do j = js, je
643  do i = is, ie
644  do k = 1, ka
645  var_z(k,i,j) = interp_xih2zh_coef(k,i,j,1) * var(interp_xih2zh_idx(k,i,j,1),i,j) &
646  + interp_xih2zh_coef(k,i,j,2) * var(interp_xih2zh_idx(k,i,j,2),i,j) &
647  + interp_xih2zh_coef(k,i,j,3) * const_undef
648  enddo
649  enddo
650  enddo
651 
652  else ! cubic spline
653 
654  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
655  !$omp private(k,i,j) &
656  !$omp private(CDZ,MD,U,V,kk,c1,c2,c3,d) &
657  !$omp shared(KA,KS,KE,kmax,IS,IE,JS,JE) &
658  !$omp shared(var,var_Z) &
659  !$omp shared(INTERP_xih2zh_idx,INTERP_xih2zh_coef,Xih,Zh,CONST_UNDEF)
660  do j = js, je
661  do i = is, ie
662 
663  do k = ks, ke
664  cdz(k) = zh(k,i,j) - zh(k-1,i,j)
665  end do
666 
667  md(ks) = 2.0 * ( cdz(ks) + cdz(ks+1) ) + cdz(ks)
668  do k = ks+1, ke-2
669  md(k) = 2.0 * ( cdz(k) + cdz(k+1) )
670  end do
671  md(ke-1) = 2.0 * ( cdz(ke-1) + cdz(ke) ) + cdz(ke)
672 
673  do k = ks, ke-1
674  v(k) = ( var(k+1,i,j) - var(k ,i,j) ) / cdz(k+1) &
675  - ( var(k ,i,j) - var(k-1,i,j) ) / cdz(k )
676  end do
677 
678  call matrix_solver_tridiagonal( kmax-1, 1, kmax-1, & ! [IN]
679  cdz(ks+1:ke ), & ! [IN]
680  md(ks :ke-1), & ! [IN]
681  cdz(ks+1:ke ), & ! [IN]
682  v(ks :ke-1), & ! [IN]
683  u(ks :ke-1) ) ! [OUT]
684 
685  u(ks-1) = u(ks)
686  u(ke) = u(ke-1)
687 
688  do k = 1, ka
689  kk = min(interp_xih2zh_idx(k,i,j,1),ke-1)
690 
691  c3 = ( u(kk+1) - u(kk) ) / cdz(kk+1)
692  c2 = 3.0_rp * u(kk)
693  c1 = ( var(kk+1,i,j) - var(kk,i,j) ) / cdz(kk+1) - ( u(kk) * 2.0_rp + u(kk+1) ) * cdz(kk+1)
694  d = xih(k) - zh(kk,i,j)
695 
696  var_z(k,i,j) = ( interp_xih2zh_coef(k,i,j,3) ) * const_undef &
697  + ( 1.0_rp-interp_xih2zh_coef(k,i,j,3) ) * ( ( ( c3 * d + c2 ) * d + c1 ) * d + var(kk,i,j) )
698  end do
699 
700  end do
701  end do
702 
703  end if
704 
705  return
706  end subroutine interp_vert_xih2zh
707 
708  !-----------------------------------------------------------------------------
709  subroutine interp_vert_zh2xih( &
710  KA, KS, KE, &
711  IA, IS, IE, &
712  JA, JS, JE, &
713  Zh, &
714  Xih, &
715  var, &
716  var_Xi )
717  use scale_const, only: &
719  use scale_matrix, only: &
720  matrix_solver_tridiagonal
721  implicit none
722 
723  integer, intent(in) :: KA, KS, KE
724  integer, intent(in) :: IA, IS, IE
725  integer, intent(in) :: JA, JS, JE
726  real(RP), intent(in) :: Zh (0:ka,ia,ja)
727  real(RP), intent(in) :: Xih (0:ka)
728  real(RP), intent(in) :: var (ka,ia,ja)
729  real(RP), intent(out) :: var_Xi(ka,ia,ja)
730 
731  real(RP) :: CDZ(ka)
732  real(RP) :: MD (ka)
733  real(RP) :: U (ka)
734  real(RP) :: V (ka)
735  real(RP) :: c1, c2, c3, d
736  integer :: kmax
737 
738  integer :: k, i, j, kk
739  !---------------------------------------------------------------------------
740 
741  kmax = ke - ks + 1
742 
743  if ( kmax == 2 ) then
744 
745  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
746  !$omp private(k,i,j) &
747  !$omp shared(KA,IS,IE,JS,JE) &
748  !$omp shared(var,var_Xi) &
749  !$omp shared(INTERP_zh2xih_idx,INTERP_zh2xih_coef,CONST_UNDEF)
750  do j = js, je
751  do i = is, ie
752  do k = 1, ka
753  var_xi(k,i,j) = interp_zh2xih_coef(k,i,j,1) * var(interp_zh2xih_idx(k,i,j,1),i,j) &
754  + interp_zh2xih_coef(k,i,j,2) * var(interp_zh2xih_idx(k,i,j,2),i,j) &
755  + interp_zh2xih_coef(k,i,j,3) * const_undef
756  enddo
757  enddo
758  enddo
759 
760  else ! cubic spline
761 
762  do k = ks, ke
763  cdz(k) = xih(k) - xih(k-1)
764  end do
765 
766  md(ks) = 2.0 * ( cdz(ks) + cdz(ks+1) ) + cdz(ks)
767  do k = ks+1, ke-2
768  md(k) = 2.0 * ( cdz(k) + cdz(k+1) )
769  end do
770  md(ke-1) = 2.0 * ( cdz(ke-1) + cdz(ke) ) + cdz(ke)
771 
772  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
773  !$omp private(k,i,j) &
774  !$omp private(U,V,kk,c1,c2,c3,d) &
775  !$omp shared(KA,KS,KE,kmax,IS,IE,JS,JE) &
776  !$omp shared(var,var_Xi,MD,CDZ) &
777  !$omp shared(INTERP_zh2xih_idx,INTERP_zh2xih_coef,Xih,Zh,CONST_UNDEF)
778  do j = js, je
779  do i = is, ie
780 
781  do k = ks, ke-1
782  v(k) = ( var(k+1,i,j) - var(k ,i,j) ) / cdz(k+1) &
783  - ( var(k ,i,j) - var(k-1,i,j) ) / cdz(k )
784  end do
785 
786  call matrix_solver_tridiagonal( kmax-1, 1, kmax-1, & ! [IN]
787  cdz(ks+1:ke ), & ! [IN]
788  md(ks :ke-1), & ! [IN]
789  cdz(ks+1:ke ), & ! [IN]
790  v(ks :ke-1), & ! [IN]
791  u(ks :ke-1) ) ! [OUT]
792 
793  u(ks-1) = u(ks)
794  u(ke) = u(ke-1)
795 
796  do k = 1, ka
797  kk = min(interp_zh2xih_idx(k,i,j,1),ke-1)
798 
799  c3 = ( u(kk+1) - u(kk) ) / cdz(kk+1)
800  c2 = 3.0_rp * u(kk)
801  c1 = ( var(kk+1,i,j) - var(kk,i,j) ) / cdz(kk+1) - ( u(kk) * 2.0_rp + u(kk+1) ) * cdz(kk+1)
802  d = zh(k,i,j) - xih(kk)
803  var_xi(k,i,j) = interp_zh2xih_coef(k,i,j,3) * const_undef &
804  + ( 1.0_rp - interp_zh2xih_coef(k,i,j,3) ) &
805  * ( ( ( c3 * d + c2 ) * d + c1 ) * d + var(kk,i,j) )
806  end do
807 
808  end do
809  end do
810 
811  end if
812 
813  return
814  end subroutine interp_vert_zh2xih
815 
816  !-----------------------------------------------------------------------------
818  subroutine interp_vert_alloc_pres( &
819  Kpres, IA, JA )
820  implicit none
821 
822  integer, intent(in) :: Kpres
823  integer, intent(in) :: IA
824  integer, intent(in) :: JA
825  !---------------------------------------------------------------------------
826 
827  allocate( interp_xi2p_idx(kpres,ia,ja,2) )
828  allocate( interp_xi2p_coef(kpres,ia,ja,3) )
829  allocate( interp_xih2p_idx(kpres,ia,ja,2) )
830  allocate( interp_xih2p_coef(kpres,ia,ja,3) )
831 
832  return
833  end subroutine interp_vert_alloc_pres
834 
835  !-----------------------------------------------------------------------------
836  subroutine interp_vert_setcoef_pres( &
837  Kpres, &
838  KA, KS, KE, &
839  IA, IS, IE, &
840  JA, JS, JE, &
841  PRES, &
842  PRESh, &
843  SFC_PRES, &
844  Paxis )
845  implicit none
846 
847  integer, intent(in) :: Kpres
848  integer, intent(in) :: KA, KS, KE
849  integer, intent(in) :: IA, IS, IE
850  integer, intent(in) :: JA, JS, JE
851  real(RP), intent(in) :: PRES ( ka,ia,ja) ! pressure in Xi coordinate [Pa]
852  real(RP), intent(in) :: PRESh (0:ka,ia,ja) ! pressure in Xi coordinate [Pa], layer interface
853  real(RP), intent(in) :: SFC_PRES( ia,ja) ! surface pressure [Pa]
854  real(RP), intent(in) :: Paxis (kpres) ! pressure level to output [Pa]
855 
856  real(RP) :: LnPRES ( ka,ia,ja) ! (log) pressure in Xi coordinate [Pa]
857  real(RP) :: LnPRESh (0:ka,ia,ja) ! (log) pressure in Xi coordinate [Pa], layer interface
858  real(RP) :: LnSFC_PRES( ia,ja) ! (log) surface pressure [Pa]
859  real(RP) :: LnPaxis (kpres) ! (log) pressure level to output [Pa]
860 
861  integer :: k, i, j, kk, kp
862  !---------------------------------------------------------------------------
863 
864  do j = js, je
865  do i = is, ie
866  do k = ks, ke
867  lnpres(k,i,j) = log( pres(k,i,j) )
868  enddo
869  enddo
870  enddo
871 
872  do j = js, je
873  do i = is, ie
874  do k = ks-1, ke
875  lnpresh(k,i,j) = log( presh(k,i,j) )
876  enddo
877  enddo
878  enddo
879 
880  do j = js, je
881  do i = is, ie
882  lnsfc_pres(i,j) = log( sfc_pres(i,j) )
883  enddo
884  enddo
885 
886  lnpaxis(:) = log( paxis(:) )
887 
888  ! full level
889 
890  do j = js, je
891  do i = is, ie
892  do k = 1, kpres
893  if ( lnpaxis(k) >= lnsfc_pres(i,j) ) then
894 
895  interp_xi2p_idx(k,i,j,1) = ks ! dummmy
896  interp_xi2p_idx(k,i,j,2) = ks ! dummmy
897  interp_xi2p_coef(k,i,j,1) = 0.0_rp
898  interp_xi2p_coef(k,i,j,2) = 0.0_rp
899  interp_xi2p_coef(k,i,j,3) = 1.0_rp ! set UNDEF
900 
901  elseif( lnpaxis(k) >= lnpres(ks,i,j) ) then
902 
903  interp_xi2p_idx(k,i,j,1) = ks ! dummmy
904  interp_xi2p_idx(k,i,j,2) = ks
905  interp_xi2p_coef(k,i,j,1) = 0.0_rp
906  interp_xi2p_coef(k,i,j,2) = 1.0_rp
907  interp_xi2p_coef(k,i,j,3) = 0.0_rp
908 
909  elseif( lnpaxis(k) < lnpres(ke,i,j) ) then
910 
911  interp_xi2p_idx(k,i,j,1) = ke ! dummmy
912  interp_xi2p_idx(k,i,j,2) = ke ! dummmy
913  interp_xi2p_coef(k,i,j,1) = 0.0_rp
914  interp_xi2p_coef(k,i,j,2) = 0.0_rp
915  interp_xi2p_coef(k,i,j,3) = 1.0_rp ! set UNDEF
916 
917  else
918 
919  do kk = ks+1, ke
920  kp = kk
921  if( lnpaxis(k) >= lnpres(kk,i,j) ) exit
922  enddo
923 
924  interp_xi2p_idx(k,i,j,1) = kp - 1
925  interp_xi2p_idx(k,i,j,2) = kp
926  interp_xi2p_coef(k,i,j,1) = ( lnpaxis(k) - lnpres(kp,i,j) ) &
927  / ( lnpres(kp-1,i,j) - lnpres(kp,i,j) )
928  interp_xi2p_coef(k,i,j,2) = ( lnpres(kp-1,i,j) - lnpaxis(k) ) &
929  / ( lnpres(kp-1,i,j) - lnpres(kp,i,j) )
930  interp_xi2p_coef(k,i,j,3) = 0.0_rp
931 
932  endif
933  enddo
934  enddo
935  enddo
936 
937  ! half level
938 
939  do j = js, je
940  do i = is, ie
941  do k = 1, kpres
942  if ( lnpaxis(k) > lnsfc_pres(i,j) ) then
943 
944  interp_xih2p_idx(k,i,j,1) = ks ! dummmy
945  interp_xih2p_idx(k,i,j,2) = ks ! dummmy
946  interp_xih2p_coef(k,i,j,1) = 0.0_rp
947  interp_xih2p_coef(k,i,j,2) = 0.0_rp
948  interp_xih2p_coef(k,i,j,3) = 1.0_rp ! set UNDEF
949 
950  elseif( lnpaxis(k) < lnpresh(ke,i,j) ) then
951 
952  interp_xih2p_idx(k,i,j,1) = ke ! dummmy
953  interp_xih2p_idx(k,i,j,2) = ke ! dummmy
954  interp_xih2p_coef(k,i,j,1) = 0.0_rp
955  interp_xih2p_coef(k,i,j,2) = 0.0_rp
956  interp_xih2p_coef(k,i,j,3) = 1.0_rp ! set UNDEF
957 
958  else
959 
960  do kk = ks, ke
961  kp = kk
962  if( lnpaxis(k) >= lnpresh(kk,i,j) ) exit
963  enddo
964 
965  interp_xih2p_idx(k,i,j,1) = kp - 1
966  interp_xih2p_idx(k,i,j,2) = kp
967  interp_xih2p_coef(k,i,j,1) = ( lnpaxis(k) - lnpresh(kp,i,j) ) &
968  / ( lnpresh(kp-1,i,j) - lnpresh(kp,i,j) )
969  interp_xih2p_coef(k,i,j,2) = ( lnpresh(kp-1,i,j) - lnpaxis(k) ) &
970  / ( lnpresh(kp-1,i,j) - lnpresh(kp,i,j) )
971  interp_xih2p_coef(k,i,j,3) = 0.0_rp
972 
973  endif
974  enddo
975  enddo
976  enddo
977 
978  return
979  end subroutine interp_vert_setcoef_pres
980 
981  !-----------------------------------------------------------------------------
982  subroutine interp_vert_xi2p( &
983  Kpres, &
984  KA, &
985  IA, IS, IE, &
986  JA, JS, JE, &
987  var, &
988  var_P )
989  use scale_const, only: &
991  implicit none
992 
993  integer, intent(in) :: Kpres
994  integer, intent(in) :: KA
995  integer, intent(in) :: IA, IS, IE
996  integer, intent(in) :: JA, JS, JE
997  real(RP), intent(in) :: var (ka ,ia,ja)
998  real(RP), intent(out) :: var_P(kpres,ia,ja)
999 
1000  integer :: k, i, j
1001  !---------------------------------------------------------------------------
1002 
1003  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
1004  !$omp private(k,i,j) &
1005  !$omp shared(Kpres,IS,IE,JS,JE) &
1006  !$omp shared(var,var_P) &
1007  !$omp shared(INTERP_xi2p_idx,INTERP_xi2p_coef,CONST_UNDEF)
1008  do j = js, je
1009  do i = is, ie
1010  do k = 1, kpres
1011  var_p(k,i,j) = interp_xi2p_coef(k,i,j,1) * var(interp_xi2p_idx(k,i,j,1),i,j) &
1012  + interp_xi2p_coef(k,i,j,2) * var(interp_xi2p_idx(k,i,j,2),i,j) &
1013  + interp_xi2p_coef(k,i,j,3) * const_undef
1014  enddo
1015  enddo
1016  enddo
1017 
1018  return
1019  end subroutine interp_vert_xi2p
1020 
1021  !-----------------------------------------------------------------------------
1022  subroutine interp_vert_xih2p( &
1023  Kpres, &
1024  KA, &
1025  IA, IS, IE, &
1026  JA, JS, JE, &
1027  var, &
1028  var_P )
1029  use scale_const, only: &
1030  const_undef
1031  implicit none
1032 
1033  integer, intent(in) :: Kpres
1034  integer, intent(in) :: KA
1035  integer, intent(in) :: IA, IS, IE
1036  integer, intent(in) :: JA, JS, JE
1037  real(RP), intent(in) :: var (ka ,ia,ja)
1038  real(RP), intent(out) :: var_P(kpres,ia,ja)
1039 
1040  integer :: k, i, j
1041  !---------------------------------------------------------------------------
1042 
1043  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
1044  !$omp private(k,i,j) &
1045  !$omp shared(Kpres,IS,IE,JS,JE) &
1046  !$omp shared(var,var_P) &
1047  !$omp shared(INTERP_xih2p_idx,INTERP_xih2p_coef,CONST_UNDEF)
1048  do j = js, je
1049  do i = is, ie
1050  do k = 1, kpres
1051  var_p(k,i,j) = interp_xih2p_coef(k,i,j,1) * var(interp_xih2p_idx(k,i,j,1),i,j) &
1052  + interp_xih2p_coef(k,i,j,2) * var(interp_xih2p_idx(k,i,j,2),i,j) &
1053  + interp_xih2p_coef(k,i,j,3) * const_undef
1054  enddo
1055  enddo
1056  enddo
1057 
1058  return
1059  end subroutine interp_vert_xih2p
1060 
1061 end module scale_interp_vert
subroutine, public interp_vert_xi2p(Kpres, KA, IA, IS, IE, JA, JS, JE, var, var_P)
subroutine, public interp_vert_z2xi(KA, KS, KE, IA, IS, IE, JA, JS, JE, Z, Xi, var, var_Xi)
logical, public interp_available
topography exists & vertical interpolation has meaning?
module INTERPOLATION vertical
subroutine, public interp_vert_xi2z(KA, KS, KE, IA, IS, IE, JA, JS, JE, Xi, Z, var, var_Z)
real(rp), public const_undef
Definition: scale_const.F90:41
module atmosphere / grid / cartesC index
module MATRIX
subroutine, public interp_vert_alloc_pres(Kpres, IA, JA)
Setup.
subroutine, public interp_vert_xih2zh(KA, KS, KE, IA, IS, IE, JA, JS, JE, Xih, Zh, var, var_Z)
subroutine, public interp_vert_setcoef_pres(Kpres, KA, KS, KE, IA, IS, IE, JA, JS, JE, PRES, PRESh, SFC_PRES, Paxis)
module CONSTANT
Definition: scale_const.F90:11
subroutine, public interp_vert_zh2xih(KA, KS, KE, IA, IS, IE, JA, JS, JE, Zh, Xih, var, var_Xi)
module profiler
Definition: scale_prof.F90:11
subroutine, public interp_vert_xih2p(Kpres, KA, IA, IS, IE, JA, JS, JE, var, var_P)
module PRECISION
subroutine, public interp_vert_setcoef(KA, KS, KE, IA, IS, IE, JA, JS, JE, TOPO_exist, Xi, Xih, Z, Zh)
Setup.
module STDIO
Definition: scale_io.F90:10