source: palm/trunk/UTIL/inifor/src/transform.f90 @ 3250

Last change on this file since 3250 was 3183, checked in by suehring, 6 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 36.2 KB
RevLine 
[2696]1!> @file src/transform.f90
2!------------------------------------------------------------------------------!
[2718]3! This file is part of the PALM model system.
[2696]4!
[2718]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
[2696]8! version.
9!
[2718]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.
[2696]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!
[2718]17! Copyright 2017-2018 Leibniz Universitaet Hannover
18! Copyright 2017-2018 Deutscher Wetterdienst Offenbach
[2696]19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
[3183]23!
24!
25! Former revisions:
26! -----------------
27! $Id: transform.f90 3183 2018-07-27 14:25:55Z sward $
[3182]28! Introduced new PALM grid stretching
29! Removed unnecessary subroutine parameters
30! Renamed kcur to k_intermediate
[2696]31!
32!
[3183]33! 3182 2018-07-27 13:36:03Z suehring
[2696]34! Initial revision
35!
36!
37!
38! Authors:
39! --------
40! @author Eckhard Kadasch
41!
42! Description:
43! ------------
44!> The transform module provides INIFOR's low-level transformation and
45!> interpolation routines. The rotated-pole transformation routines phirot2phi,
46!> phi2phirot, rlarot2rla, rla2rlarot, uv2uvrot, and uvrot2uv are adapted from
47!> int2lm's utility routines.
48!------------------------------------------------------------------------------!
49 MODULE transform
50
51    USE control
52    USE defs,                                                                  &
53        ONLY: TO_DEGREES, TO_RADIANS, PI, dp
54    USE types
55    USE util,                                                                  &       
56        ONLY: real_to_str, str
57
58    IMPLICIT NONE
59
60 CONTAINS
61
62!------------------------------------------------------------------------------!
63! Description:
64! ------------
65!> Interpolates linearly in the vertical direction in very column (i,j) of the
66!> output array outvar(i,j,:) using values of the source array invar. In cells
67!> that are outside the COSMO-DE domain, indicated by negative interpolation
68!> weights, extrapolate constantly from the cell above.
69!>
70!> Input parameters:
71!> -----------------
72!> invar : Array of source data
73!>
74!> outgrid % kk : Array of vertical neighbour indices. kk(i,j,k,:) contain the
75!>     indices of the two vertical neighbors of PALM-4U point (i,j,k) on the
76!>     input grid corresponding to the source data invar.
77!>
78!> outgrid % w_verti : Array of weights for vertical linear interpolation
79!>     corresponding to neighbour points indexed by kk.
80!>
81!> Output papameters:
82!> ------------------
83!> outvar : Array of interpolated data
84!------------------------------------------------------------------------------!
85    SUBROUTINE interpolate_1d(in_arr, out_arr, outgrid)
86       TYPE(grid_definition), INTENT(IN) ::  outgrid
87       REAL(dp), INTENT(IN)              ::  in_arr(0:,0:,0:)
[3182]88       REAL(dp), INTENT(OUT)             ::  out_arr(0:,0:,:)
[2696]89
[3182]90       INTEGER :: i, j, k, l, nz
[2696]91
92       nz = UBOUND(out_arr, 3)
93
[3182]94       DO j = LBOUND(out_arr, 2), UBOUND(out_arr, 2)
95       DO i = LBOUND(out_arr, 1), UBOUND(out_arr, 1)
96       DO k = nz, LBOUND(out_arr, 3), -1
[2696]97
98          ! TODO: Remove IF clause and extrapolate based on a critical vertical
99          ! TODO: index marking the lower bound of COSMO-DE data coverage.
100          ! Check for negative interpolation weights indicating grid points
101          ! below COSMO-DE domain and extrapolate from the top in such cells.
102          IF (outgrid % w_verti(i,j,k,1) < -1.0_dp .AND. k < nz)  THEN
103             out_arr(i,j,k) = out_arr(i,j,k+1)
104          ELSE
105             out_arr(i,j,k) = 0.0_dp
106             DO l = 1, 2
[3182]107                out_arr(i,j,k) = out_arr(i,j,k) +                              &
108                    outgrid % w_verti(i,j,k,l) *                               &
[2696]109                    in_arr(i,j,outgrid % kk(i,j,k, l) )
110             END DO
111          END IF
112       END DO
113       END DO
114       END DO
115    END SUBROUTINE interpolate_1d
116
117
118!------------------------------------------------------------------------------!
119! Description:
120! ------------
121!> Interpolates bi-linearly in horizontal planes on every k level of the output
122!> array outvar(:,:,k) using values of the source array invar(:,:,:). The source
123!> (invar) and interpolation array (outvar) need to have matching dimensions.
124!>
125!> Input parameters:
126!> -----------------
127!> invar : Array of source data
128!>
129!> outgrid % ii, % jj : Array of neighbour indices in x and y direction.
130!>     ii(i,j,k,:), and jj(i,j,k,:) contain the four horizontal neighbour points
131!>     of PALM-4U point (i,j,k) on the input grid corresponding to the source
132!>     data invar. (The outgrid carries the relationship with the ingrid in the
133!      form of the interpoaltion weights.)
134!>
135!> outgrid % w_horiz: Array of weights for horizontal bi-linear interpolation
136!>     corresponding to neighbour points indexed by ii and jj.
137!>
138!> Output papameters:
139!> ------------------
140!> outvar : Array of interpolated data
141!------------------------------------------------------------------------------!
142    SUBROUTINE interpolate_2d(invar, outvar, outgrid, ncvar)
143    ! I index 0-based for the indices of the outvar to be consistent with the
144    ! outgrid indices and interpolation weights.
[3182]145       TYPE(grid_definition), INTENT(IN)  ::  outgrid
146       REAL(dp), INTENT(IN)               ::  invar(0:,0:,0:)
147       REAL(dp), INTENT(OUT)              ::  outvar(0:,0:,0:)
[2696]148       TYPE(nc_var), INTENT(IN), OPTIONAL ::  ncvar
149
150       INTEGER ::  i, j, k, l
151
152       ! TODO: check if input dimensions are consistent, i.e. ranges are correct
153       IF (UBOUND(outvar, 3) .GT. UBOUND(invar, 3))  THEN
154           message = "Output array for '" // TRIM(ncvar % name) // "' has ' more levels (" // &
155              TRIM(str(UBOUND(outvar, 3))) // ") than input variable ("//&
156              TRIM(str(UBOUND(invar, 3))) // ")."
157           CALL abort('interpolate_2d', message) 
158       END IF
159
160       DO k = 0, UBOUND(outvar, 3)
161       DO j = 0, UBOUND(outvar, 2)
162       DO i = 0, UBOUND(outvar, 1)
163          outvar(i,j,k) = 0.0_dp
164          DO l = 1, 4
165             
166             outvar(i,j,k) = outvar(i,j,k) +                                   &
167                outgrid % w_horiz(i,j,l) * invar( outgrid % ii(i,j,l),         & 
168                                                  outgrid % jj(i,j,l),         &
169                                                  k )
170          END DO
171       END DO
172       END DO
173       END DO
174       
175    END SUBROUTINE interpolate_2d
176
177
178    SUBROUTINE average_2d(in_arr, out_arr, imin, imax, jmin, jmax)
179       REAL(dp), INTENT(IN)  ::  in_arr(0:,0:,0:)
180       REAL(dp), INTENT(OUT) ::  out_arr(0:)
181       INTEGER, INTENT(IN)   ::  imin, imax, jmin, jmax
182
183       INTEGER  ::  i, j, k
184       REAL(dp) ::  ni
185       
186       IF (imin < 0)  CALL abort('average_2d', "imin < 0.")
187       IF (jmin < 0)  CALL abort('average_2d', "jmin < 0.")
188       IF (imax > UBOUND(in_arr, 1))  CALL abort('average_2d', "imax out of i bound.")
189       IF (jmax > UBOUND(in_arr, 2))  CALL abort('average_2d', "jmax out of j bound.")
190
191       DO k = 0, UBOUND(out_arr, 1)
192
193          out_arr(k) = 0.0_dp
194          DO j = jmin, jmax
195          DO i = imin, imax
196             out_arr(k) = out_arr(k) + in_arr(i, j, k)
197          END DO
198          END DO
199
200       END DO
201   
202       ! devide by number of grid points
203       ni = 1.0_dp / ( (imax - imin + 1) * (jmax - jmin + 1) )
204       out_arr(:) = out_arr(:) * ni
205
206    END SUBROUTINE average_2d
207
208
209    SUBROUTINE interpolate_3d(source_array, palm_array, palm_intermediate, palm_grid)
210       TYPE(grid_definition), INTENT(IN) ::  palm_intermediate, palm_grid
211       REAL(dp), DIMENSION(:,:,:), INTENT(IN)  ::  source_array
212       REAL(dp), DIMENSION(:,:,:), INTENT(OUT) ::  palm_array
213       REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::  intermediate_array
214       INTEGER ::  nx, ny, nz
215
216       nx = palm_intermediate % nx
217       ny = palm_intermediate % ny
218       nz = palm_intermediate % nz ! nlev
219
220       ! Interpolate from COSMO-DE to intermediate grid. Allocating with one
221       ! less point in the vertical, since scalars like T have 50 instead of 51
222       ! points in COSMO-DE.
223       ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nz-1)) !
224
225       CALL interpolate_2d(source_array, intermediate_array, palm_intermediate)
226
227       ! Interpolate from intermediate grid to palm_grid grid, includes
228       ! extrapolation for cells below COSMO-DE domain.
229       CALL interpolate_1d(intermediate_array, palm_array, palm_grid)
230
231       DEALLOCATE(intermediate_array)
232
233    END SUBROUTINE interpolate_3d
234
235
236    SUBROUTINE average_profile(source_array, profile_array, imin, imax, jmin, jmax,&
237                               palm_intermediate, palm_grid)
238       TYPE(grid_definition), INTENT(IN)          ::  palm_intermediate, palm_grid
239       REAL(dp), DIMENSION(:,:,:), INTENT(IN)     ::  source_array
240       INTEGER, INTENT(IN)                        ::  imin, imax, jmin, jmax
241       REAL(dp), DIMENSION(0:,0:,0:), INTENT(OUT) ::  profile_array
242       REAL(dp), DIMENSION(:,:,:), ALLOCATABLE    ::  intermediate_array
243       INTEGER                                    ::  nx, ny, nz
244
245       nx = palm_intermediate % nx
246       ny = palm_intermediate % ny
247       nz = palm_intermediate % nz
248       ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nz-1))
249       intermediate_array(:,:,:) = 0.0_dp
250
251       ! average input array to intermediate profile
252       CALL average_2d(source_array, intermediate_array(0,0,:), imin, imax, jmin, jmax)
253
254       ! vertically interpolate to ouput array
255       CALL interpolate_1d(intermediate_array, profile_array, palm_grid)
256
257       DEALLOCATE(intermediate_array)
258
259    END SUBROUTINE average_profile
260
261
262
263!-----------------------------------------------------------------------------!
264! Description:
265! -----------
266!> This routine computes the inverse Plate Carree projection, i.e. in projects
267!> Cartesian coordinates (x,y) onto a sphere. It returns the latitude and
268!> lngitude of a geographical system centered at x0 and y0.
269!-----------------------------------------------------------------------------!
270    SUBROUTINE inv_plate_carree(x, y, x0, y0, r, lat, lon)
271       REAL(dp), INTENT(IN)  ::  x(:), y(:), x0, y0, r
272       REAL(dp), INTENT(OUT) ::  lat(:), lon(:)
273       
274       REAL(dp) :: ri
275
276       ! TODO check dimensions of lat/lon and y/x match
277
278       ri = 1.0_dp / r
279       
280       lat(:) = (y(:) - y0) * ri
281       lon(:) = (x(:) - x0) * ri
282    END SUBROUTINE 
283
284
285!-----------------------------------------------------------------------------!
286! Description:
287! ------------
288!> Computes the reverse Plate-Carree projection of a x or y position on a
289!> Cartesian grid.
290!>
291!> Input parameters:
292!> -----------------
293!> xy : x or y coordinate of the Cartasian grid point [m].
294!>
295!> xy0 : x or y coordinate that coincides with the origin of the underlying
296!>     sperical system (crossing point of the equator and prime meridian) [m].
297!>
298!> r : Radius of the of the underlying sphere, e.g. EARTH_RADIUS [m].
299!>
300!> Returns:
301!> --------
302!> project : Longitude (in case xy = x) or latitude (xy = y) of the given input
303!>     coordinate xy.
304!------------------------------------------------------------------------------!
305    ELEMENTAL REAL(dp) FUNCTION project(xy, xy0, r)
306       REAL(dp), INTENT(IN)  ::  xy, xy0, r
307       REAL(dp) :: ri
308
309       ! If this elemental function is called with a large array as xy, it is
310       ! computationally more efficient to precompute the inverse radius and
311       ! then muliply.
312       ri = 1.0_dp / r
313
314       project = (xy - xy0) * ri
315
316    END FUNCTION project
317
318
319    REAL(dp) FUNCTION phic_to_phin(phi_c)
320        REAL(dp), INTENT(IN) ::  phi_c
321
322        phic_to_phin = 0.5_dp * PI - ABS(phi_c)
323
324    END FUNCTION phic_to_phin
325
326   
327    REAL(dp) FUNCTION lamc_to_lamn(phi_c, lam_c)
328       REAL(dp), INTENT(IN) ::  phi_c, lam_c
329       
330       lamc_to_lamn = lam_c
331       IF (phi_c > 0.0_dp)  THEN
332          lamc_to_lamn = lam_c - SIGN(PI, lam_c)
333       END IF
334
335    END FUNCTION lamc_to_lamn
336
337
338    REAL(dp) FUNCTION gamma_from_hemisphere(phi_cg, phi_ref)
339       REAL(dp), INTENT(IN) ::  phi_cg, phi_ref
340       LOGICAL ::  palm_centre_is_south_of_cosmo_origin
341       
342       palm_centre_is_south_of_cosmo_origin = (phi_cg < phi_ref)
343
344       IF (palm_centre_is_south_of_cosmo_origin)  THEN
345           gamma_from_hemisphere = PI
346       ELSE
347           gamma_from_hemisphere = 0.0_dp
348       END IF
349    END FUNCTION gamma_from_hemisphere
350
351
352!------------------------------------------------------------------------------!
353! Description:
354! ------------
355!> Computes the geographical coordinates corresponding to the given rotated-pole
356!> coordinates.
357!>
358!> In INIFOR, this routine is used to convert coordinates between two
359!> rotated-pole systems: COSMO-DE's rotated-pole system, and one centred at the
360!> PALM-4U domain centre. In this case, the PALM-4U system is thought of as the
361!> rotated-pole system and the routine is used to rotate back to COSMO-DE's
362!> system which is thought of as the geographical one.
363!>
364!> Input parameters:
365!> -----------------
366!> phir(:), lamr(: ): latitudes and longitudes of the rotated-pole grid
367!>
368!> phip, lamp: latitude and longitude of the rotated north pole
369!>
370!> gam: "angle between the north poles. If [gam] is not present, the other
371!>       system is the real geographical system." (original phiro2rot
372!>       description)
373!>
374!> Output parameters:
375!> ------------------
376!> phi(:,:), lam(:,:): geographical latitudes and logitudes
377!------------------------------------------------------------------------------!
378    SUBROUTINE rotate_to_cosmo(phir, lamr, phip, lamp, phi, lam, gam)
379       REAL(dp), INTENT(IN)  ::  phir(0:), lamr(0:), phip, lamp, gam
380       REAL(dp), INTENT(OUT) ::  phi(0:,0:), lam(0:,0:)
381
382       INTEGER ::  i, j
383       
384       IF ( SIZE(phi, 1) .NE. SIZE(lam, 1) .OR. &
385            SIZE(phi, 2) .NE. SIZE(lam, 2) )  THEN
386          PRINT *, "inifor: rotate_to_cosmo: Dimensions of phi and lambda do not match. Dimensions are:"
387          PRINT *, "inifor: rotate_to_cosmo: phi: ", SIZE(phi, 1), SIZE(phi, 2)
388          PRINT *, "inifor: rotate_to_cosmo: lam: ", SIZE(lam, 1), SIZE(lam, 2)
389          STOP
390       END IF
391
392       IF ( SIZE(phir) .NE. SIZE(phi, 2) .OR. &
393            SIZE(lamr) .NE. SIZE(phi, 1) )  THEN
394          PRINT *, "inifor: rotate_to_cosmo: Dimensions of phir and lamr do not match. Dimensions are:"
395          PRINT *, "inifor: rotate_to_cosmo: phir: ", SIZE(phir), SIZE(phi, 2)
396          PRINT *, "inifor: rotate_to_cosmo: lamr: ", SIZE(lamr), SIZE(phi, 1)
397          STOP
398       END IF
399       
400       DO j = 0, UBOUND(phir, 1)
401          DO i = 0, UBOUND(lamr, 1)
402
403             phi(i,j) = phirot2phi(phir(j) * TO_DEGREES,                       &
404                                   lamr(i) * TO_DEGREES,                       &
405                                   phip * TO_DEGREES,                          &
406                                   lamp * TO_DEGREES,                          &
407                                   gam  * TO_DEGREES) * TO_RADIANS
408
409             lam(i,j) = rlarot2rla(phir(j) * TO_DEGREES,                       &
410                                   lamr(i) * TO_DEGREES,                       &
411                                   phip * TO_DEGREES,                          &
412                                   lamp * TO_DEGREES,                          &
413                                   gam  * TO_DEGREES) * TO_RADIANS
414
415          END DO
416       END DO
417
418    END SUBROUTINE rotate_to_cosmo
[3182]419       
[2696]420
421
422!------------------------------------------------------------------------------!
423! Description:
424! ------------
425!> Compute indices of PALM-4U grid point neighbours in the target
426!> system (COSMO-DE) by rounding up and down. (i,j) are the indices of
427!> the PALM-4U grid and (ii(i,j,1-4), jj(i,j,1-4)) contain the indices
428!> of the its four neigbouring points in the COSMO-DE grid.
429!>
430!>
431!>                     COSMO-DE grid
432!>                     -------------
433!>           jj, lat
[3182]434!>              ^        j
435!>              |         \          i
[2696]436!>  jj(i,j,2/3) + ... 2 ---\--------/------ 3
437!>              |     | ^   \      /        |
438!>              |     | |wp  \    /         |
439!>              |     | v     \  /          |
440!>       latpos + ............ o/ (i,j)     |
441!>              |     |        :            |
442!>              |     |        :<----wl---->|
443!>  jj(i,j,1/4) + ... 1 -------:----------- 4
444!>              |     :        :            :
445!>              |     :        :            :
446!>              |     :      lonpos         :
447!>              L-----+--------+------------+------> ii, lon
448!>               ii(i,j,1/2)        ii(i,j,3/4)
449!>
450!>
451!> Input parameters:
452!> -----------------
453!> source_lat, source_lon : (rotated-pole) coordinates of the source grid (e.g.
454!>    COSMO-DE)
455!>
456!> source_dxi, source_dyi : inverse grid spacings of the source grid.
457!>
458!> target_lat, target_lon : (rotated-pole) coordinates of the target grid (e.g.
459!>    COSMO-DE)
460!>
461!> Output parameters:
462!> ------------------
463!> palm_ii, palm_jj : x and y index arrays of horizontal neighbour columns
464!>
465!------------------------------------------------------------------------------!
[3182]466    SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon,                &
467                                          palm_clat, palm_clon,                &
468                                          palm_ii, palm_jj)
[2696]469
470       REAL(dp), DIMENSION(0:), INTENT(IN)        ::  cosmo_lat, cosmo_lon
471       REAL(dp), DIMENSION(0:,0:), INTENT(IN)     ::  palm_clat, palm_clon
[3182]472       REAL(dp)                                   ::  cosmo_dxi, cosmo_dyi
[2696]473       INTEGER, DIMENSION(0:,0:,1:), INTENT(OUT)  ::  palm_ii, palm_jj
474
475       REAL(dp) ::  lonpos, latpos, lon0, lat0
476       INTEGER  ::  i, j
477
478       lon0 = cosmo_lon(0)
479       lat0 = cosmo_lat(0)
[3182]480       cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
481       cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
[2696]482
483       DO j = 0, UBOUND(palm_clon, 2)!palm_grid % ny
484       DO i = 0, UBOUND(palm_clon, 1)!palm_grid % nx
485          ! Compute the floating point index corrseponding to PALM-4U grid point
486          ! location along target grid (COSMO-DE) axes.
487          lonpos = (palm_clon(i,j) - lon0) * cosmo_dxi
488          latpos = (palm_clat(i,j) - lat0) * cosmo_dyi
489
490          IF (lonpos < 0.0 .OR. latpos < 0.0)  THEN
491             PRINT *, " Error while finding neighbours: lonpos or latpos out of bounds!"
492             PRINT *, "     (i,j) = (", i, ",",j,")"
493             PRINT *, "      lonpos ", lonpos*TO_DEGREES, ", latpos ", latpos*TO_DEGREES
494             PRINT *, "        lon0 ", lon0  *TO_DEGREES,   ", lat0   ", lat0*TO_DEGREES
495             PRINT *, "    PALM lon ", palm_clon(i,j)*TO_DEGREES,   ", PALM lat ",palm_clat(i,j)*TO_DEGREES
496             STOP
497          END IF
498
499          palm_ii(i,j,1) = FLOOR(lonpos)
500          palm_ii(i,j,2) = FLOOR(lonpos)
501          palm_ii(i,j,3) = CEILING(lonpos)
502          palm_ii(i,j,4) = CEILING(lonpos)
503
504          palm_jj(i,j,1) = FLOOR(latpos)
505          palm_jj(i,j,2) = CEILING(latpos)
506          palm_jj(i,j,3) = CEILING(latpos)
507          palm_jj(i,j,4) = FLOOR(latpos)
508       END DO
509       END DO
510
511    END SUBROUTINE find_horizontal_neighbours
512
513   
514    SUBROUTINE find_vertical_neighbours_and_weights(palm_grid, palm_intermediate)
515       TYPE(grid_definition), INTENT(INOUT) ::  palm_grid
516       TYPE(grid_definition), INTENT(IN)    ::  palm_intermediate
517
[3182]518       INTEGER  ::  i, j, k, nx, ny, nz, nlev, k_intermediate
[2696]519       LOGICAL  ::  point_is_below_grid, point_is_above_grid,                  &
520                    point_is_in_current_cell
521       REAL(dp) ::  current_height, column_base, column_top, h_top, h_bottom,  &
522                    weight
523
524       nx   = palm_grid % nx
525       ny   = palm_grid % ny
526       nz   = palm_grid % nz
527       nlev = palm_intermediate % nz
528
529       ! in each column of the fine grid, find vertical neighbours of every cell
530       DO i = 0, nx
531       DO j = 0, ny
532
[3182]533          k_intermediate = 0
[2696]534
535          column_base = palm_intermediate % h(i,j,0)
536          column_top  = palm_intermediate % h(i,j,nlev)
537
538          ! scan through palm_grid column until and set neighbour indices in
539          ! case current_height is either below column_base, in the current
540          ! cell, or above column_top. Keep increasing current cell index until
541          ! the current cell overlaps with the current_height.
[3182]542          DO k = 1, nz
[2696]543
544             ! Memorize the top and bottom boundaries of the coarse cell and the
545             ! current height within it
546             current_height = palm_grid % z(k) + palm_grid % z0
[3182]547             h_top    = palm_intermediate % h(i,j,k_intermediate+1)
548             h_bottom = palm_intermediate % h(i,j,k_intermediate)
[2696]549
550             point_is_above_grid = (current_height > column_top) !22000m, very unlikely
551             point_is_below_grid = (current_height < column_base)
552
553             point_is_in_current_cell = (                                      &
554                current_height >= h_bottom .AND.                               &
555                current_height <  h_top                                        &
556             )
557
558             ! set default weights
559             palm_grid % w_verti(i,j,k,1:2) = 0.0_dp
560
561             IF (point_is_above_grid)  THEN
562
563                palm_grid % kk(i,j,k,1:2) = nlev
564                palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp
565
[3182]566                message = "PALM-4U grid extends above COSMO-DE model top."
567                CALL abort('find_vertical_neighbours_and_weights', message)
568
[2696]569             ELSE IF (point_is_below_grid)  THEN
570
571                palm_grid % kk(i,j,k,1:2) = 0
572                palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp
573
574             ELSE
575                ! cycle through intermediate levels until current
576                ! intermediate-grid cell overlaps with current_height
[3182]577                DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
578                   k_intermediate = k_intermediate + 1
[2696]579
[3182]580                   h_top    = palm_intermediate % h(i,j,k_intermediate+1)
581                   h_bottom = palm_intermediate % h(i,j,k_intermediate)
[2696]582                   point_is_in_current_cell = (                                &
583                      current_height >= h_bottom .AND.                         &
584                      current_height <  h_top                                  &
585                   )
586                END DO
587
[3182]588                ! k_intermediate = 48 indicates the last section (indices 48 and 49), i.e.
589                ! k_intermediate = 49 is not the beginning of a valid cell.
590                IF (k_intermediate > nlev-1)  THEN
591                   message = "Index " // TRIM(str(k_intermediate)) //          &
592                             " is above intermediate grid range."
[2696]593                   CALL abort('find_vertical_neighbours', message)
594                END IF
595   
[3182]596                palm_grid % kk(i,j,k,1) = k_intermediate
597                palm_grid % kk(i,j,k,2) = k_intermediate + 1
[2696]598
599                ! copmute vertical weights
600                weight = (h_top - current_height) / (h_top - h_bottom)
601                palm_grid % w_verti(i,j,k,1) = weight
602                palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight
603             END IF
604
605          END DO
606
607       END DO
608       END DO
609
610    END SUBROUTINE find_vertical_neighbours_and_weights
611
612!------------------------------------------------------------------------------!
613! Description:
614! ------------
615!> Compute the four weights for horizontal bilinear interpolation given the
616!> coordinates clon(i,j) clat(i,j) of the PALM-4U grid in the COSMO-DE
617!> rotated-pole grid and the neightbour indices ii(i,j,1-4) and jj(i,j,1-4).
618!>
619!> Input parameters:
620!> -----------------
621!> palm_grid % clon : longitudes of PALM-4U scalars (cell centres) in COSMO-DE's rotated-pole grid [rad]
622!>
623!> palm_grid % clat : latitudes of PALM-4U cell centres in COSMO-DE's rotated-pole grid [rad]
624!>
625!> cosmo_grid % lon : rotated-pole longitudes of scalars (cell centres) of the COSMO-DE grid [rad]
626!>
627!> cosmo_grid % lat : rotated-pole latitudes of scalars (cell centers) of the COSMO-DE grid [rad]
628!>
629!> cosmo_grid % dxi : inverse grid spacing in the first dimension [m^-1]
630!>
631!> cosmo_grid % dyi : inverse grid spacing in the second dimension [m^-1]
632!>
633!> Output parameters:
634!> ------------------
635!> palm_grid % w_horiz(:,:,1-4) : weights for bilinear horizontal interpolation
636!
637!                               COSMO-DE grid
638!                               -------------
639!                     jj, lat
640!                        ^        j
641!                        |         \          i
642!            jj(i,j,2/3) + ... 2 ---\--------/------ 3
643!                        |     | ^   \      /        |
644!                        |     | |wp  \    /         |
645!                        |     | v     \  /          |
646!                 latpos + ............ o/ (i,j)     |
647!                        |     |        :            |
648!                        |     |        :<----wl---->|
649!            jj(i,j,1/4) + ... 1 -------:----------- 4
650!                        |     :        :            :
651!                        |     :        :            :
652!                        |     :      lonpos         :
653!                        L-----+--------+------------+------> ii, lon
654!                         ii(i,j,1/2)        ii(i,j,3/4)
655!         
656    SUBROUTINE compute_horizontal_interp_weights(cosmo_lat, cosmo_lon,         &
[3182]657       palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)
[2696]658       
659       REAL(dp), DIMENSION(0:), INTENT(IN)        ::  cosmo_lat, cosmo_lon
[3182]660       REAL(dp)                                   ::  cosmo_dxi, cosmo_dyi
[2696]661       REAL(dp), DIMENSION(0:,0:), INTENT(IN)     ::  palm_clat, palm_clon
662       INTEGER, DIMENSION(0:,0:,1:), INTENT(IN)   ::  palm_ii, palm_jj
663
664       REAL(dp), DIMENSION(0:,0:,1:), INTENT(OUT) ::  palm_w_horiz
665
666       REAL(dp) ::  wl, wp
667       INTEGER  ::  i, j
668
[3182]669       cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
670       cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
671
[2696]672       DO j = 0, UBOUND(palm_clon, 2)
673       DO i = 0, UBOUND(palm_clon, 1)
674     
675          ! weight in lambda direction
676          wl = ( cosmo_lon(palm_ii(i,j,4)) - palm_clon(i,j) ) * cosmo_dxi
677
678          ! weight in phi direction
679          wp = ( cosmo_lat(palm_jj(i,j,2)) - palm_clat(i,j) ) * cosmo_dyi
680
681          IF (wl > 1.0_dp .OR. wl < 0.0_dp)  THEN
682              message = "Horizontal weight wl = " // TRIM(real_to_str(wl)) //   &
683                        " is out bounds."
684              CALL abort('compute_horizontal_interp_weights', message)
685          END IF
686          IF (wp > 1.0_dp .OR. wp < 0.0_dp)  THEN
687              message = "Horizontal weight wp = " // TRIM(real_to_str(wp)) //   &
688                        " is out bounds."
689              CALL abort('compute_horizontal_interp_weights', message)
690          END IF
691
692          palm_w_horiz(i,j,1) = wl * wp
693          palm_w_horiz(i,j,2) = wl * (1.0_dp - wp)
694          palm_w_horiz(i,j,3) = (1.0_dp - wl) * (1.0_dp - wp)
695          palm_w_horiz(i,j,4) = 1.0_dp - SUM( palm_w_horiz(i,j,1:3) )
696
697       END DO
698       END DO
699       
700    END SUBROUTINE compute_horizontal_interp_weights
701
702
703!------------------------------------------------------------------------------!
704! Description:
705! ------------
706!> Interpolates u and v components of velocities located at cell faces to the
707!> cell centres by averaging neighbouring values.
708!>
709!> This routine is designed to be used with COSMO-DE arrays where there are the
710!> same number of grid points for scalars (centres) and velocities (faces). In
711!> COSMO-DE the velocity points are staggared one half grid spaceing up-grid
712!> which means the first centre point has to be omitted and is set to zero.
713    SUBROUTINE centre_velocities(u_face, v_face, u_centre, v_centre)
714       REAL(dp), DIMENSION(0:,0:,0:), INTENT(IN)  ::  u_face, v_face
715       REAL(dp), DIMENSION(0:,0:,0:), INTENT(OUT) ::  u_centre, v_centre
716       INTEGER ::  nx, ny
717
718       nx = UBOUND(u_face, 1)
719       ny = UBOUND(u_face, 2)
720
721       u_centre(0,:,:)  = 0.0_dp
722       u_centre(1:,:,:) = 0.5_dp * ( u_face(0:nx-1,:,:) + u_face(1:,:,:) )
723
724       v_centre(:,0,:)  = 0.0_dp
725       v_centre(:,1:,:) = 0.5_dp * ( v_face(:,0:ny-1,:) + v_face(:,1:,:) )
726    END SUBROUTINE centre_velocities
727
728
729    FUNCTION phirot2phi (phirot, rlarot, polphi, pollam, polgam)
730   
731       REAL(dp), INTENT (IN) ::  polphi      !< latitude of the rotated north pole
732       REAL(dp), INTENT (IN) ::  pollam      !< longitude of the rotated north pole
733       REAL(dp), INTENT (IN) ::  phirot      !< latitude in the rotated system
734       REAL(dp), INTENT (IN) ::  rlarot      !< longitude in the rotated system
735       REAL(dp), INTENT (IN) ::  polgam      !< angle between the north poles of the systems
736
737       REAL(dp)              ::  phirot2phi  !< latitude in the geographical system
738       
739       REAL(dp)              ::  zsinpol, zcospol, zphis, zrlas, zarg, zgam
740   
741       zsinpol = SIN(polphi * TO_RADIANS)
742       zcospol = COS(polphi * TO_RADIANS)
743       zphis   = phirot * TO_RADIANS
744
745       IF (rlarot > 180.0_dp)  THEN
746          zrlas = rlarot - 360.0_dp
747       ELSE
748          zrlas = rlarot
749       END IF
750       zrlas = zrlas * TO_RADIANS
751     
752       IF (polgam /= 0.0_dp)  THEN
753          zgam = polgam * TO_RADIANS
754          zarg = zsinpol * SIN (zphis) +                                       &
755                 zcospol * COS(zphis) * ( COS(zrlas) * COS(zgam) -             &
756                                          SIN(zgam)  * SIN(zrlas) )
757       ELSE
758          zarg = zcospol * COS (zphis) * COS (zrlas) + zsinpol * SIN (zphis)
759       END IF
760     
761       phirot2phi = ASIN (zarg) * TO_DEGREES
762   
763    END FUNCTION phirot2phi
764
765
766    FUNCTION phi2phirot (phi, rla, polphi, pollam)
767   
768       REAL(dp), INTENT (IN) ::  polphi !< latitude of the rotated north pole
769       REAL(dp), INTENT (IN) ::  pollam !< longitude of the rotated north pole
770       REAL(dp), INTENT (IN) ::  phi    !< latitude in the geographical system
771       REAL(dp), INTENT (IN) ::  rla    !< longitude in the geographical system
772       
773       REAL(dp) ::  phi2phirot          !< longitude in the rotated system
774       
775       REAL(dp) ::  zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1
776       
777       zsinpol = SIN(polphi * TO_RADIANS)
778       zcospol = COS(polphi * TO_RADIANS)
779       zlampol = pollam * TO_RADIANS
780       zphi    = phi * TO_RADIANS
781
782       IF (rla > 180.0_dp)  THEN
783          zrla1 = rla - 360.0_dp
784       ELSE
785          zrla1 = rla
786       END IF
787       zrla = zrla1 * TO_RADIANS
788       
789       zarg1 = SIN(zphi) * zsinpol
790       zarg2 = COS(zphi) * zcospol * COS(zrla - zlampol)
791       
792       phi2phirot = ASIN(zarg1 + zarg2) * TO_DEGREES
793   
794    END FUNCTION phi2phirot
795
796
797    FUNCTION rlarot2rla(phirot, rlarot, polphi, pollam, polgam)
798   
799       REAL(dp), INTENT (IN) ::  polphi !< latitude of the rotated north pole
800       REAL(dp), INTENT (IN) ::  pollam !< longitude of the rotated north pole
801       REAL(dp), INTENT (IN) ::  phirot !< latitude in the rotated system
802       REAL(dp), INTENT (IN) ::  rlarot !< longitude in the rotated system
803       REAL(dp), INTENT (IN) ::  polgam !< angle between the north poles of the systems
804       
805       REAL(dp) ::  rlarot2rla          !< latitude in the geographical system
806       
807       REAL(dp) ::  zsinpol, zcospol, zlampol, zphis, zrlas, zarg1, zarg2, zgam
808       
809       zsinpol = SIN(TO_RADIANS * polphi)
810       zcospol = COS(TO_RADIANS * polphi)
811       zlampol = TO_RADIANS * pollam
812       zphis   = TO_RADIANS * phirot
813
814       IF (rlarot > 180.0_dp)  THEN
815          zrlas = rlarot - 360.0_dp
816       ELSE
817          zrlas = rlarot
818       END IF
819       zrlas   = TO_RADIANS * zrlas
820     
821       IF (polgam /= 0.0_dp)  THEN
822          zgam  = TO_RADIANS * polgam
823          zarg1 = SIN(zlampol) * (zcospol * SIN(zphis) - zsinpol*COS(zphis) *  &
824                  (COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) -        &
825                  COS(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) +       &
826                                                COS(zrlas) * SIN(zgam) )
827       
828          zarg2 = COS (zlampol) * (zcospol * SIN(zphis) - zsinpol*COS(zphis) * &
829                  (COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) +        &
830                  SIN(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) +       &
831                                                COS(zrlas) * SIN(zgam) )
832       ELSE
833          zarg1   = SIN (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis)  +     &
834                                      zcospol *              SIN(zphis)) -     &
835                    COS (zlampol) *             SIN(zrlas) * COS(zphis)
836          zarg2   = COS (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis)  +     &
837                                      zcospol *              SIN(zphis)) +     &
838                    SIN (zlampol) *             SIN(zrlas) * COS(zphis)
839       END IF
840     
841       IF (zarg2 == 0.0_dp)  zarg2 = 1.0E-20_dp
842     
843       rlarot2rla = ATAN2(zarg1,zarg2) * TO_DEGREES
844       
845    END FUNCTION rlarot2rla
846
847
848    FUNCTION rla2rlarot ( phi, rla, polphi, pollam, polgam )
849
850       REAL(dp), INTENT (IN) ::  polphi !< latitude of the rotated north pole
851       REAL(dp), INTENT (IN) ::  pollam !< longitude of the rotated north pole
852       REAL(dp), INTENT (IN) ::  phi    !< latitude in geographical system
853       REAL(dp), INTENT (IN) ::  rla    !< longitude in geographical system
854       REAL(dp), INTENT (IN) ::  polgam !< angle between the north poles of the systems
855       
856       REAL (KIND=dp) ::  rla2rlarot    !< latitude in the the rotated system
857       
858       REAL (KIND=dp) ::  zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1
859       
860       zsinpol = SIN(polphi * TO_RADIANS)
861       zcospol = COS(polphi * TO_RADIANS)
862       zlampol = pollam * TO_RADIANS
863       zphi    = phi * TO_RADIANS
864
865       IF (rla > 180.0_dp)  THEN
866          zrla1 = rla - 360.0_dp
867       ELSE
868          zrla1 = rla
869       END IF
870       zrla = zrla1 * TO_RADIANS
871       
872       zarg1 = - SIN (zrla-zlampol) * COS(zphi)
873       zarg2 = - zsinpol * COS(zphi) * COS(zrla-zlampol) + zcospol * SIN(zphi)
874       
875       IF (zarg2 == 0.0_dp)  zarg2 = 1.0E-20_dp
876       
877       rla2rlarot = ATAN2 (zarg1,zarg2) * TO_DEGREES
878       
879       IF (polgam /= 0.0_dp )  THEN
880          rla2rlarot = polgam + rla2rlarot
881          IF (rla2rlarot > 180._dp)  rla2rlarot = rla2rlarot - 360.0_dp
882       END IF
883       
884    END FUNCTION rla2rlarot
885
886
887    SUBROUTINE uv2uvrot(u, v, rlat, rlon, pollat, pollon, urot, vrot)
888   
889       REAL(dp), INTENT (IN)  ::  u, v           !< wind components in the true geographical system
890       REAL(dp), INTENT (IN)  ::  rlat, rlon     !< coordinates in the true geographical system
891       REAL(dp), INTENT (IN)  ::  pollat, pollon !< latitude and longitude of the north pole of the rotated grid
892       
893       REAL(dp), INTENT (OUT) ::  urot, vrot     !< wind components in the rotated grid             
894       
895       REAL (dp) ::  zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm
896       
897       zsinpol = SIN(pollat * TO_RADIANS)
898       zcospol = COS(pollat * TO_RADIANS)
899       zlonp   = (pollon-rlon) * TO_RADIANS
900       zlat    = rlat * TO_RADIANS
901       
902       zarg1 = zcospol * SIN(zlonp)
903       zarg2 = zsinpol * COS(zlat) - zcospol * SIN(zlat) * COS(zlonp)
904       znorm = 1.0_dp / SQRT(zarg1*zarg1 + zarg2*zarg2)
905       
906       urot = u * zarg2 * znorm - v * zarg1 * znorm
907       vrot = u * zarg1 * znorm + v * zarg2 * znorm
908   
909    END SUBROUTINE uv2uvrot
910
911
912    SUBROUTINE uvrot2uv (urot, vrot, rlat, rlon, pollat, pollon, u, v)
913   
914       REAL(dp), INTENT(IN) ::  urot, vrot     !< wind components in the rotated grid
915       REAL(dp), INTENT(IN) ::  rlat, rlon     !< latitude and longitude in the true geographical system
916       REAL(dp), INTENT(IN) ::  pollat, pollon !< latitude and longitude of the north pole of the rotated grid
917       
918       REAL(dp), INTENT(OUT) ::  u, v          !< wind components in the true geographical system
919       
920       REAL(dp) ::  zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm
921     
922       zsinpol = SIN(pollat * TO_RADIANS)
923       zcospol = COS(pollat * TO_RADIANS)
924       zlonp   = (pollon-rlon) * TO_RADIANS
925       zlat    = rlat * TO_RADIANS
926     
927       zarg1 = zcospol * SIN(zlonp)
928       zarg2 = zsinpol * COS(zlat) - zcospol * SIN(zlat) * COS(zlonp)
929       znorm = 1.0_dp / SQRT(zarg1*zarg1 + zarg2*zarg2)
930     
931       u =   urot * zarg2 * znorm + vrot * zarg1 * znorm
932       v = - urot * zarg1 * znorm + vrot * zarg2 * znorm
933   
934    END SUBROUTINE uvrot2uv
935
936 END MODULE
937
Note: See TracBrowser for help on using the repository browser.