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

Last change on this file since 3182 was 3182, checked in by suehring, 3 years ago

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

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