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

Last change on this file since 2696 was 2696, checked in by kanani, 6 years ago

Merge of branch palm4u into trunk

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