Home

Context Navigation

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

Last change on this file since 3182 was 3182, checked in by suehring, 6 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 ", lonposTO_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) - zsinpolCOS(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) - zsinpolCOS(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(zarg1zarg1 + zarg2zarg2)
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(zarg1zarg1 + zarg2zarg2)
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.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |