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

Last change on this file since 4313 was 3785, checked in by eckhard, 3 years ago

inifor: Removed unused variables, improved coding style

  • 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-2019 Leibniz Universitaet Hannover
18! Copyright 2017-2019 Deutscher Wetterdienst Offenbach
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: test-interpolation.f90 3785 2019-03-06 10:41:14Z suehring $
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.