Home

Context Navigation

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

Last change on this file since 3359 was 3183, checked in by suehring, 6 years ago
last commit documented
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	!
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: transform.f90 3183 2018-07-27 14:25:55Z knoop $
28	! Introduced new PALM grid stretching
29	! Removed unnecessary subroutine parameters
30	! Renamed kcur to k_intermediate
31	!
32	!
33	! 3182 2018-07-27 13:36:03Z suehring
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:)
88	REAL(dp), INTENT(OUT) :: out_arr(0:,0:,:)
89
90	INTEGER :: i, j, k, l, nz
91
92	nz = UBOUND(out_arr, 3)
93
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
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
107	out_arr(i,j,k) = out_arr(i,j,k) + &
108	outgrid % w_verti(i,j,k,l) * &
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.
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:)
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
419
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
434	!> ^ j
435	!> \| \ i
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	!------------------------------------------------------------------------------!
466	SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon, &
467	palm_clat, palm_clon, &
468	palm_ii, palm_jj)
469
470	REAL(dp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon
471	REAL(dp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon
472	REAL(dp) :: cosmo_dxi, cosmo_dyi
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)
480	cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
481	cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
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 ", lonposTO_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
518	INTEGER :: i, j, k, nx, ny, nz, nlev, k_intermediate
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
533	k_intermediate = 0
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.
542	DO k = 1, nz
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
547	h_top = palm_intermediate % h(i,j,k_intermediate+1)
548	h_bottom = palm_intermediate % h(i,j,k_intermediate)
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
566	message = "PALM-4U grid extends above COSMO-DE model top."
567	CALL abort('find_vertical_neighbours_and_weights', message)
568
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
577	DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
578	k_intermediate = k_intermediate + 1
579
580	h_top = palm_intermediate % h(i,j,k_intermediate+1)
581	h_bottom = palm_intermediate % h(i,j,k_intermediate)
582	point_is_in_current_cell = ( &
583	current_height >= h_bottom .AND. &
584	current_height < h_top &
585	)
586	END DO
587
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."
593	CALL abort('find_vertical_neighbours', message)
594	END IF
595
596	palm_grid % kk(i,j,k,1) = k_intermediate
597	palm_grid % kk(i,j,k,2) = k_intermediate + 1
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, &
657	palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)
658
659	REAL(dp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon
660	REAL(dp) :: cosmo_dxi, cosmo_dyi
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
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
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) - zsinpolCOS(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) - zsinpolCOS(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(zarg1zarg1 + zarg2zarg2)
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(zarg1zarg1 + zarg2zarg2)
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.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |