source: palm/trunk/UTIL/inifor/tests/test-interpolation.f90

Last change on this file was 4843, checked in by raasch, 21 months ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

  • Property svn:keywords set to Id
File size: 9.6 KB
Line 
1!> @file tests/test-interpolation.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2017-2021 Leibniz Universitaet Hannover
18! Copyright 2017-2021 Deutscher Wetterdienst Offenbach
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: test-interpolation.f90 4843 2021-01-15 15:22:11Z banzhafs $
28! Prefixed all INIFOR modules with inifor_
29!
30!
31! 3183 2018-07-27 14:25:55Z suehring
32! Updated test for new grid_definition
33!
34!
35! 3182 2018-07-27 13:36:03Z suehring
36! Initial revision
37!
38!
39!
40! Authors:
41! --------
42! @author Eckhard Kadasch
43!
44! Description:
45! ------------
46!> This program tests INIFOR's horizontal interpolation.
47!------------------------------------------------------------------------------!
48 PROGRAM test_interpolation
49
50    USE inifor_grid,                                                           &
51        ONLY:  grid_definition, init_grid_definition, TO_RADIANS, TO_DEGREES,  &
52               linspace, hhl
53    USE inifor_transform,                                                      &
54        ONLY:  find_horizontal_neighbours, compute_horizontal_interp_weights
55    USE test_utils
56   
57    IMPLICIT NONE
58
59!
60!------------------------------------------------------------------------------
61!- Test 1: Find neighbours
62!------------------------------------------------------------------------------
63    CHARACTER(LEN=30) ::  title = "find neighbours"
64    LOGICAL           ::  res
65
66    TYPE(grid_definition) ::  palm_grid, cosmo_grid
67    INTEGER               ::  i, j, ii_ref(0:1, 0:1, 4), jj_ref(0:1, 0:1, 4)
68    INTEGER, PARAMETER    ::  nlon=3, nlat=3, nlev=2
69    REAL                  ::  w_ref(4), lat(0:2), lon(0:2)
70
71    title = "find neighbours"
72    CALL begin_test(title, res)
73
74    ! Arange.
75    ! Make a COSMO-DE grid with just two horizotal cells/three h. points
76    PRINT *, "INIT GRID"
77
78    ! Allocate grid.hhl for use in init_grid_definition. In INIFOR, this is done
79    ! in get_netcdf_variable_2d. In this test, grid.hhl is not used and only
80    ! defined manually because it is used in init_grid_definition.
81    ALLOCATE (hhl (nlon, nlat, nlev) )
82    hhl(:,:,:) = 0.0
83    CALL init_grid_definition('cosmo-de', grid = cosmo_grid,                   &
84                              xmin = -5.0 * TO_RADIANS, xmax = 5.5 * TO_RADIANS, &
85                              ymin = -5.0 * TO_RADIANS, ymax = 6.5 * TO_RADIANS, &
86                              x0 = 0.0, y0 = 0.0, z0 = 0.0,                    &
87                              nx = nlon-1, ny = nlat-1, nz = nlev-1)
88
89    PRINT *, "GRID DONE"
90    PRINT *, "COSMO lats: ", cosmo_grid % lat * TO_DEGREES
91    PRINT *, "COSMO lons: ", cosmo_grid % lon * TO_DEGREES
92   
93    res = assert_equal( (/cosmo_grid%lat(0), cosmo_grid % lon(0),              &
94                          cosmo_grid%lat(2), cosmo_grid % lon(2),              &
95                          (cosmo_grid%lon(1) - cosmo_grid%lon(0))*TO_DEGREES,  &
96                          (cosmo_grid%lat(1) - cosmo_grid%lat(0))*TO_DEGREES/),&
97                        (/-5.0 * TO_RADIANS, -5.0 * TO_RADIANS,                &
98                           6.5 * TO_RADIANS,  5.5 * TO_RADIANS,                &
99                           5.25, 5.75 /),    &
100                        "COSMO grid coordinates" )
101    ! Define a PALM-4U grid with only one cell, i.e. four points in the
102    ! horizontal plane. The points are located at the centres of
103    ! the COSMO-DE cells.
104    CALL init_grid_definition('palm intermediate', grid = palm_grid,           &
105                              xmin = 0.0, xmax = 1.0,                          &
106                              ymin = 0.0, ymax = 1.0,                          &
107                              x0 = 0.0, y0 = 0.0, z0 = 0.0,                    &
108                              nx = 1, ny = 1, nz = 1)
109
110    palm_grid % clon(0,0) = 0.5 * cosmo_grid % lon(0)
111    palm_grid % clat(0,0) = 0.5 * cosmo_grid % lat(0)
112   
113    palm_grid % clon(0,1) = 0.5 * cosmo_grid % lon(0)
114    palm_grid % clat(0,1) = 0.5 * cosmo_grid % lat(2)
115
116    palm_grid % clon(1,1) = 0.5 * cosmo_grid % lon(2)
117    palm_grid % clat(1,1) = 0.5 * cosmo_grid % lat(2)
118
119    palm_grid % clon(1,0) = 0.5 * cosmo_grid % lon(2)
120    palm_grid % clat(1,0) = 0.5 * cosmo_grid % lat(0)
121
122    ii_ref(0,0,:) = (/0, 0, 1, 1/)
123    jj_ref(0,0,:) = (/0, 1, 1, 0/)
124
125    ii_ref(0,1,:) = (/0, 0, 1, 1/)
126    jj_ref(0,1,:) = (/1, 2, 2, 1/)
127   
128    ii_ref(1,1,:) = (/1, 1, 2, 2/)
129    jj_ref(1,1,:) = (/1, 2, 2, 1/)
130
131    ii_ref(1,0,:) = (/1, 1, 2, 2/)
132    jj_ref(1,0,:) = (/0, 1, 1, 0/)
133
134    ! Act
135    CALL find_horizontal_neighbours(cosmo_grid % lat, cosmo_grid % lon,        &
136                                    palm_grid % clat, palm_grid % clon,        &
137                                    palm_grid % ii, palm_grid % jj)
138
139    ! Assert
140    DO j = 0, 1
141    DO i = 0, 1
142       res = res .AND. ALL(palm_grid%ii(i,j,:) == ii_ref(i,j,:))
143       PRINT *, "ii     : ", palm_grid%ii(i,j,:)
144       PRINT *, "ii_ref : ", ii_ref(i,j,:), " indices match? ", res
145       res = res .AND. ALL(palm_grid%jj(i,j,:) == jj_ref(i,j,:))
146       PRINT *, "jj     : ", palm_grid%jj(i,j,:)
147       PRINT *, "jj_ref : ", jj_ref(i,j,:), " indices match? ", res
148    ENDDO
149    ENDDO
150
151    CALL end_test(title, res)
152   
153
154!
155!------------------------------------------------------------------------------
156!- Test 2: Compute weights for linear interpolation
157!------------------------------------------------------------------------------
158    title = "interpolation weights"
159    CALL begin_test(title, res)
160
161    ! Arange
162    ! defining some shorthands
163    lon(:) = cosmo_grid % lon(:)
164    lat(:) = cosmo_grid % lat(:)
165
166    ! set up PALM-4U points at 1/4 and 1/3 of the COSMO grid widths
167    palm_grid % clon(0,0) = -0.25  * (lon(1) - lon(0)) + lon(1)
168    palm_grid % clat(0,0) = -2./3. * (lat(1) - lat(0)) + lat(1)
169   
170    palm_grid % clon(0,1) = -2./3. * (lon(1) - lon(0)) + lon(1)
171    palm_grid % clat(0,1) = +0.25  * (lat(2) - lat(1)) + lat(1)
172
173    palm_grid % clon(1,1) = +0.25  * (lon(2) - lon(1)) + lon(1)
174    palm_grid % clat(1,1) = +2./3. * (lat(2) - lat(1)) + lat(1)
175
176    palm_grid % clon(1,0) = +2./3. * (lon(2) - lon(1)) + lon(1)
177    palm_grid % clat(1,0) = -0.25  * (lat(1) - lat(0)) + lat(1)
178
179    DO j = 0, 1
180    DO i = 0, 1
181       PRINT *, "PALM lon, lat: ", palm_grid % clon(i,j) * TO_DEGREES, palm_grid % clat(i,j)*TO_DEGREES
182    ENDDO
183    ENDDO
184
185    ! Act
186    CALL find_horizontal_neighbours(cosmo_grid % lat, cosmo_grid % lon,        &
187                                    palm_grid % clat, palm_grid % clon,        &
188                                    palm_grid % ii, palm_grid % jj)
189
190    CALL compute_horizontal_interp_weights(cosmo_grid % lat, cosmo_grid % lon, &
191                                           palm_grid % clat, palm_grid % clon, &
192                                           palm_grid % ii, palm_grid % jj,     &
193                                           palm_grid % w_horiz)
194
195    ! Assert
196    ! asserting that neighbours are still correct
197    DO j = 0, 1
198    DO i = 0, 1
199       res = res .AND. ALL(palm_grid%ii(i,j,:) == ii_ref(i,j,:))
200       PRINT *, "ii     : ", palm_grid%ii(i,j,:)
201       PRINT *, "ii_ref : ", ii_ref(i,j,:), " indices match? ", res
202       res = res .AND. ALL(palm_grid%jj(i,j,:) == jj_ref(i,j,:))
203       PRINT *, "jj     : ", palm_grid%jj(i,j,:)
204       PRINT *, "jj_ref : ", jj_ref(i,j,:), " indices match? ", res
205    ENDDO
206    ENDDO
207
208    ! asserting that all four weights equal, 0.5, 0.25, 1./6., and 1./12., resp.
209    w_ref = (/1./6., 1./12., 0.25, 0.5/)
210    res = res .AND. assert_equal(palm_grid % w_horiz(0, 0, :), w_ref(:), "weights at (0,0)")
211    !res = res .AND. palm_grid % w_horiz(0, 0, 1) == w_ref(1)
212    !res = res .AND. palm_grid % w_horiz(0, 0, 2) == w_ref(2)
213    !res = res .AND. palm_grid % w_horiz(0, 0, 3) == w_ref(3)
214    !res = res .AND. palm_grid % w_horiz(0, 0, 4) == w_ref(4)
215
216    w_ref = (/0.5, 1./6., 1./12., 0.25/)
217    res = res .AND. assert_equal(palm_grid % w_horiz(0, 1, :), w_ref(:), "weights at (0,1)")
218    !res = res .AND. palm_grid % w_horiz(0, 1, 1) == w_ref(4)
219    !res = res .AND. palm_grid % w_horiz(0, 1, 2) == w_ref(1)
220    !res = res .AND. palm_grid % w_horiz(0, 1, 3) == w_ref(2)
221    !res = res .AND. palm_grid % w_horiz(0, 1, 4) == w_ref(3)
222
223    w_ref = (/0.25, 0.5, 1./6., 1./12./)
224    res = res .AND. assert_equal(palm_grid % w_horiz(1, 1, :), w_ref(:), "weights at (1,1)")
225    !res = res .AND. palm_grid % w_horiz(1, 1, 1) == w_ref(3)
226    !res = res .AND. palm_grid % w_horiz(1, 1, 2) == w_ref(4)
227    !res = res .AND. palm_grid % w_horiz(1, 1, 3) == w_ref(1)
228    !res = res .AND. palm_grid % w_horiz(1, 1, 4) == w_ref(2)
229
230    w_ref = (/1./12., 0.25, 0.5, 1./6./)
231    res = res .AND. assert_equal(palm_grid % w_horiz(1, 0, :), w_ref(:), "weights at (1,0)")
232    !res = res .AND. palm_grid % w_horiz(1, 0, 1) == w_ref(2)
233    !res = res .AND. palm_grid % w_horiz(1, 0, 2) == w_ref(3)
234    !res = res .AND. palm_grid % w_horiz(1, 0, 3) == w_ref(4)
235    !res = res .AND. palm_grid % w_horiz(1, 0, 4) == w_ref(1)
236
237    CALL end_test(title, res)
238
239 END PROGRAM test_interpolation
Note: See TracBrowser for help on using the repository browser.