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