Home

Context Navigation

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

Last change on this file since 2803 was 2718, checked in by maronga, 7 years ago
deleting of deprecated files; headers updated where needed
Property svn:keywords set to `Id`
File size: 35.4 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 2718 2018-01-02 08:49:38Z thiele $
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 ", lonposTO_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) - zsinpolCOS(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) - zsinpolCOS(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(zarg1zarg1 + zarg2zarg2)
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(zarg1zarg1 + zarg2zarg2)
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.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |