Home

Context Navigation

source: palm/trunk/UTIL/inifor/src/inifor_transform.f90 @ 3678

Last change on this file since 3678 was 3678, checked in by eckhard, 5 years ago
inifor: bugfix: avoid empty averaging regions, check if all input files are present
Property svn:keywords set to `Id`
File size: 53.6 KB

Line
1	!> @file src/inifor_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: inifor_transform.f90 3678 2019-01-17 14:12:17Z eckhard $
28	! Check if set of averaging columns is empty
29	!
30	!
31	! 3618 2018-12-10 13:25:22Z eckhard
32	! Prefixed all INIFOR modules with inifor_, removed unused variables
33	!
34	!
35	! 3615 2018-12-10 07:21:03Z raasch
36	! bugfix: abort replaced by inifor_abort
37	!
38	! 3614 2018-12-10 07:05:46Z raasch
39	! unused variables removed
40	!
41	! 3613 2018-12-07 18:20:37Z eckhard
42	! Use averaged heights profile for level-based averaging instead of modified
43	! COSMO heights array
44	!
45	!
46	! 3557 2018-11-22 16:01:22Z eckhard
47	! Updated documentation
48	!
49	!
50	! 3537 2018-11-20 10:53:14Z eckhard
51	! bugfix: working precision added
52	!
53	! 3447 2018-10-29 15:52:54Z eckhard
54	! Renamed source files for compatibilty with PALM build system
55	!
56	!
57	! 3395 2018-10-22 17:32:49Z eckhard
58	! Switched addressing of averaging regions from index bounds to list of columns
59	! Added routines for the computation of geostrophic winds including:
60	! - the downward extrapolation of density (linear) and
61	! - pressure (hydrostatic equation) and
62	! - rotation of geostrophic wind components to PALM frame of reference
63	!
64	! 3183 2018-07-27 14:25:55Z suehring
65	! Introduced new PALM grid stretching
66	! Removed unnecessary subroutine parameters
67	! Renamed kcur to k_intermediate
68	!
69	!
70	! 3182 2018-07-27 13:36:03Z suehring
71	! Initial revision
72	!
73	!
74	!
75	! Authors:
76	! --------
77	!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
78	!
79	! Description:
80	! ------------
81	!> The transform module provides INIFOR's low-level transformation and
82	!> interpolation routines. The rotated-pole transformation routines phirot2phi,
83	!> phi2phirot, rlarot2rla, rla2rlarot, uv2uvrot, and uvrot2uv are adapted from
84	!> int2lm's utility routines.
85	!------------------------------------------------------------------------------!
86	MODULE inifor_transform
87
88	USE inifor_control
89	USE inifor_defs, &
90	ONLY: G, TO_DEGREES, TO_RADIANS, PI, dp
91	USE inifor_types
92	USE inifor_util, &
93	ONLY: real_to_str, str
94
95	IMPLICIT NONE
96
97	CONTAINS
98
99	!------------------------------------------------------------------------------!
100	! Description:
101	! ------------
102	!> Interpolates linearly in the vertical direction in very column (i,j) of the
103	!> output array outvar(i,j,:) using values of the source array invar. In cells
104	!> that are outside the COSMO-DE domain, indicated by negative interpolation
105	!> weights, extrapolate constantly from the cell above.
106	!>
107	!> Input parameters:
108	!> -----------------
109	!> invar : Array of source data
110	!>
111	!> outgrid % kk : Array of vertical neighbour indices. kk(i,j,k,:) contain the
112	!> indices of the two vertical neighbors of PALM-4U point (i,j,k) on the
113	!> input grid corresponding to the source data invar.
114	!>
115	!> outgrid % w_verti : Array of weights for vertical linear interpolation
116	!> corresponding to neighbour points indexed by kk.
117	!>
118	!> Output papameters:
119	!> ------------------
120	!> outvar : Array of interpolated data
121	!------------------------------------------------------------------------------!
122	SUBROUTINE interpolate_1d(in_arr, out_arr, outgrid)
123	TYPE(grid_definition), INTENT(IN) :: outgrid
124	REAL(dp), INTENT(IN) :: in_arr(0:,0:,0:)
125	REAL(dp), INTENT(OUT) :: out_arr(0:,0:,:)
126
127	INTEGER :: i, j, k, l, nz
128
129	nz = UBOUND(out_arr, 3)
130
131	DO j = LBOUND(out_arr, 2), UBOUND(out_arr, 2)
132	DO i = LBOUND(out_arr, 1), UBOUND(out_arr, 1)
133	DO k = nz, LBOUND(out_arr, 3), -1
134
135	!
136	!-- TODO: Remove IF clause and extrapolate based on a critical vertical
137	!-- TODO: index marking the lower bound of COSMO-DE data coverage.
138	!-- Check for negative interpolation weights indicating grid points
139	!-- below COSMO-DE domain and extrapolate from the top in such cells.
140	IF (outgrid % w_verti(i,j,k,1) < -1.0_dp .AND. k < nz) THEN
141	out_arr(i,j,k) = out_arr(i,j,k+1)
142	ELSE
143	out_arr(i,j,k) = 0.0_dp
144	DO l = 1, 2
145	out_arr(i,j,k) = out_arr(i,j,k) + &
146	outgrid % w_verti(i,j,k,l) * &
147	in_arr(i,j,outgrid % kk(i,j,k, l) )
148	END DO
149	END IF
150	END DO
151	END DO
152	END DO
153	END SUBROUTINE interpolate_1d
154
155
156	!------------------------------------------------------------------------------!
157	! Description:
158	! ------------
159	!> Interpolates bi-linearly in horizontal planes on every k level of the output
160	!> array outvar(:,:,k) using values of the source array invar(:,:,:). The source
161	!> (invar) and interpolation array (outvar) need to have matching dimensions.
162	!>
163	!> Input parameters:
164	!> -----------------
165	!> invar : Array of source data
166	!>
167	!> outgrid % ii, % jj : Array of neighbour indices in x and y direction.
168	!> ii(i,j,k,:), and jj(i,j,k,:) contain the four horizontal neighbour points
169	!> of PALM-4U point (i,j,k) on the input grid corresponding to the source
170	!> data invar. (The outgrid carries the relationship with the ingrid in the
171	! form of the interpoaltion weights.)
172	!>
173	!> outgrid % w_horiz: Array of weights for horizontal bi-linear interpolation
174	!> corresponding to neighbour points indexed by ii and jj.
175	!>
176	!> Output papameters:
177	!> ------------------
178	!> outvar : Array of interpolated data
179	!------------------------------------------------------------------------------!
180	SUBROUTINE interpolate_2d(invar, outvar, outgrid, ncvar)
181	!
182	!-- I index 0-based for the indices of the outvar to be consistent with the
183	!-- outgrid indices and interpolation weights.
184	TYPE(grid_definition), INTENT(IN) :: outgrid
185	REAL(dp), INTENT(IN) :: invar(0:,0:,0:)
186	REAL(dp), INTENT(OUT) :: outvar(0:,0:,0:)
187	TYPE(nc_var), INTENT(IN), OPTIONAL :: ncvar
188
189	INTEGER :: i, j, k, l
190
191	!
192	!-- TODO: check if input dimensions are consistent, i.e. ranges are correct
193	IF ( UBOUND(outvar, 3) .GT. UBOUND(invar, 3) ) THEN
194	message = "Output array for '" // TRIM(ncvar % name) // "' has ' more levels (" // &
195	TRIM(str(UBOUND(outvar, 3))) // ") than input variable ("//&
196	TRIM(str(UBOUND(invar, 3))) // ")."
197	CALL inifor_abort('interpolate_2d', message)
198	END IF
199
200	DO k = 0, UBOUND(outvar, 3)
201	DO j = 0, UBOUND(outvar, 2)
202	DO i = 0, UBOUND(outvar, 1)
203	outvar(i,j,k) = 0.0_dp
204	DO l = 1, 4
205
206	outvar(i,j,k) = outvar(i,j,k) + &
207	outgrid % w_horiz(i,j,l) * invar( outgrid % ii(i,j,l), &
208	outgrid % jj(i,j,l), &
209	k )
210	END DO
211	END DO
212	END DO
213	END DO
214
215	END SUBROUTINE interpolate_2d
216
217
218	!------------------------------------------------------------------------------!
219	! Description:
220	! ------------
221	!> Compute the horizontal average of the in_arr(:,:,:) and store it in
222	!> out_arr(:)
223	!------------------------------------------------------------------------------!
224	SUBROUTINE average_2d(in_arr, out_arr, ii, jj)
225	REAL(dp), INTENT(IN) :: in_arr(0:,0:,0:)
226	REAL(dp), INTENT(OUT) :: out_arr(0:)
227	INTEGER, INTENT(IN), DIMENSION(:) :: ii, jj
228
229	INTEGER :: i, j, k, l
230	REAL(dp) :: ni
231
232	IF (SIZE(ii) /= SIZE(jj)) THEN
233	message = "Length of 'ii' and 'jj' index lists do not match." // &
234	NEW_LINE(' ') // "ii has " // str(SIZE(ii)) // " elements, " // &
235	NEW_LINE(' ') // "jj has " // str(SIZE(jj)) // "."
236	CALL inifor_abort('average_2d', message)
237	END IF
238
239	IF (SIZE(ii) == 0) THEN
240	message = "No columns to average over; " // &
241	"size of index lists 'ii' and 'jj' is zero."
242	CALL inifor_abort('average_2d', message)
243	END IF
244
245	DO k = 0, UBOUND(out_arr, 1)
246
247	out_arr(k) = 0.0_dp
248	DO l = 1, UBOUND(ii, 1)
249	i = ii(l)
250	j = jj(l)
251	out_arr(k) = out_arr(k) + in_arr(i, j, k)
252	END DO
253
254	END DO
255
256	ni = 1.0_dp / SIZE(ii)
257	out_arr(:) = out_arr(:) * ni
258
259	END SUBROUTINE average_2d
260
261
262	!------------------------------------------------------------------------------!
263	! Description:
264	! ------------
265	!> Three-dimensional interpolation driver. Interpolates from the source_array to
266	!> the given palm_grid.
267	!>
268	!> The routine separates horizontal and vertical
269	!> interpolation. In the horizontal interpolation step, the source_array data is
270	!> interpolated along COSMO grid levels onto the intermediate grid (vertically
271	!> as coarse as COSMO, horizontally as fine as PALM).
272	!------------------------------------------------------------------------------!
273	SUBROUTINE interpolate_3d(source_array, palm_array, palm_intermediate, palm_grid)
274	TYPE(grid_definition), INTENT(IN) :: palm_intermediate, palm_grid
275	REAL(dp), DIMENSION(:,:,:), INTENT(IN) :: source_array
276	REAL(dp), DIMENSION(:,:,:), INTENT(OUT) :: palm_array
277	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: intermediate_array
278	INTEGER :: nx, ny, nlev
279
280	nx = palm_intermediate % nx
281	ny = palm_intermediate % ny
282	nlev = palm_intermediate % nz
283
284	!
285	!-- Interpolate from COSMO to intermediate grid. Allocating with one
286	!-- less point in the vertical, since scalars like T have 50 instead of 51
287	!-- points in COSMO.
288	ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nlev-1)) !
289
290	CALL interpolate_2d(source_array, intermediate_array, palm_intermediate)
291
292	!
293	!-- Interpolate from intermediate grid to palm_grid grid, includes
294	!-- extrapolation for cells below COSMO domain.
295	CALL interpolate_1d(intermediate_array, palm_array, palm_grid)
296
297	DEALLOCATE(intermediate_array)
298
299	END SUBROUTINE interpolate_3d
300
301
302	!------------------------------------------------------------------------------!
303	! Description:
304	! ------------
305	!> Average data horizontally from the source_array over the region given by the
306	!> averaging grid 'avg_grid' and store the result in 'profile_array'.
307	!------------------------------------------------------------------------------!
308	SUBROUTINE average_profile(source_array, profile_array, avg_grid)
309	TYPE(grid_definition), INTENT(IN) :: avg_grid
310	REAL(dp), DIMENSION(:,:,:), INTENT(IN) :: source_array
311	REAL(dp), DIMENSION(:), INTENT(OUT) :: profile_array
312
313	INTEGER :: i_source, j_source, k_profile, k_source, l, m
314
315	REAL :: ni_columns
316
317	profile_array(:) = 0.0_dp
318
319	DO l = 1, avg_grid % n_columns
320	i_source = avg_grid % iii(l)
321	j_source = avg_grid % jjj(l)
322
323	!
324	!-- Loop over PALM levels
325	DO k_profile = avg_grid % k_min, UBOUND(profile_array, 1)
326
327	!
328	!-- Loop over vertical interpolation neighbours
329	DO m = 1, 2
330
331	k_source = avg_grid % kkk(l, k_profile, m)
332
333	profile_array(k_profile) = profile_array(k_profile) &
334	+ avg_grid % w(l, k_profile, m) &
335	* source_array(i_source, j_source, k_source)
336	!
337	!-- Loop over vertical interpolation neighbours m
338	END DO
339
340	!
341	!-- Loop over PALM levels k_profile
342	END DO
343
344	!
345	!-- Loop over horizontal neighbours l
346	END DO
347
348	ni_columns = 1.0_dp / avg_grid % n_columns
349	profile_array(:) = profile_array(:) * ni_columns
350
351	!
352	!-- Constant extrapolation to the bottom
353	profile_array(1:avg_grid % k_min-1) = profile_array(avg_grid % k_min)
354
355	END SUBROUTINE average_profile
356
357
358	!------------------------------------------------------------------------------!
359	! Description:
360	! ------------
361	!> Extrapolates density linearly from the level 'k_min' downwards.
362	!------------------------------------------------------------------------------!
363	SUBROUTINE extrapolate_density(rho, avg_grid)
364	REAL(dp), DIMENSION(:), INTENT(INOUT) :: rho
365	TYPE(grid_definition), INTENT(IN) :: avg_grid
366
367	REAL(dp) :: drhodz, dz, zk, rhok
368	INTEGER :: k_min
369
370	k_min = avg_grid % k_min
371	zk = avg_grid % z(k_min)
372	rhok = rho(k_min)
373	dz = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
374	drhodz = (rho(k_min + 1) - rho(k_min)) / dz
375
376	rho(1:k_min-1) = rhok + drhodz * (avg_grid % z(1:k_min-1) - zk)
377
378	END SUBROUTINE extrapolate_density
379
380
381	!------------------------------------------------------------------------------!
382	! Description:
383	! ------------
384	!> Driver for extrapolating pressure from PALM level k_min downwards
385	!------------------------------------------------------------------------------!
386	SUBROUTINE extrapolate_pressure(p, rho, avg_grid)
387	REAL(dp), DIMENSION(:), INTENT(IN) :: rho
388	REAL(dp), DIMENSION(:), INTENT(INOUT) :: p
389	TYPE(grid_definition), INTENT(IN) :: avg_grid
390
391	REAL(dp) :: drhodz, dz, zk, rhok
392	INTEGER :: k, k_min
393
394	k_min = avg_grid % k_min
395	zk = avg_grid % z(k_min)
396	rhok = rho(k_min)
397	dz = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
398	drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz
399
400	DO k = 1, k_min-1
401	p(k) = constant_density_pressure(p(k_min), zk, rhok, drhodz, &
402	avg_grid % z(k), G)
403	END DO
404
405	END SUBROUTINE extrapolate_pressure
406
407
408	!------------------------------------------------------------------------------!
409	! Description:
410	! ------------
411	!> Takes the averaged pressure profile <p> and sets the lowest entry to the
412	!> extrapolated pressure at the surface.
413	!------------------------------------------------------------------------------!
414	SUBROUTINE get_surface_pressure(p, rho, avg_grid)
415	REAL(dp), DIMENSION(:), INTENT(IN) :: rho
416	REAL(dp), DIMENSION(:), INTENT(INOUT) :: p
417	TYPE(grid_definition), INTENT(IN) :: avg_grid
418
419	REAL(dp) :: drhodz, dz, zk, rhok
420	INTEGER :: k_min
421
422	k_min = avg_grid % k_min
423	zk = avg_grid % z(k_min)
424	rhok = rho(k_min)
425	dz = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
426	drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz
427
428	p(1) = constant_density_pressure(p(k_min), zk, rhok, drhodz, &
429	0.0_dp, G)
430
431	END SUBROUTINE get_surface_pressure
432
433
434	FUNCTION constant_density_pressure(pk, zk, rhok, drhodz, z, g) RESULT(p)
435
436	REAL(dp), INTENT(IN) :: pk, zk, rhok, drhodz, g, z
437	REAL(dp) :: p
438
439	p = pk + ( zk - z ) * g * ( rhok + 0.5drhodz (zk - z) )
440
441	END FUNCTION constant_density_pressure
442
443	!-----------------------------------------------------------------------------!
444	! Description:
445	! -----------
446	!> This routine computes profiles of the two geostrophic wind components ug and
447	!> vg.
448	!-----------------------------------------------------------------------------!
449	SUBROUTINE geostrophic_winds(p_north, p_south, p_east, p_west, rho, f3, &
450	Lx, Ly, phi_n, lam_n, phi_g, lam_g, ug, vg)
451
452	REAL(dp), DIMENSION(:), INTENT(IN) :: p_north, p_south, p_east, p_west, &
453	rho
454	REAL(dp), INTENT(IN) :: f3, Lx, Ly, phi_n, lam_n, phi_g, lam_g
455	REAL(dp), DIMENSION(:), INTENT(OUT) :: ug, vg
456	REAL(dp) :: facx, facy
457
458	facx = 1.0_dp / (Lx * f3)
459	facy = 1.0_dp / (Ly * f3)
460	ug(:) = - facy / rho(:) * (p_north(:) - p_south(:))
461	vg(:) = facx / rho(:) * (p_east(:) - p_west(:))
462
463	CALL rotate_vector_field( &
464	ug, vg, angle = meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g)&
465	)
466
467	END SUBROUTINE geostrophic_winds
468
469
470	!-----------------------------------------------------------------------------!
471	! Description:
472	! -----------
473	!> This routine computes the inverse Plate Carree projection, i.e. in projects
474	!> Cartesian coordinates (x,y) onto a sphere. It returns the latitude and
475	!> lngitude of a geographical system centered at x0 and y0.
476	!-----------------------------------------------------------------------------!
477	SUBROUTINE inv_plate_carree(x, y, x0, y0, r, lat, lon)
478	REAL(dp), INTENT(IN) :: x(:), y(:), x0, y0, r
479	REAL(dp), INTENT(OUT) :: lat(:), lon(:)
480
481	REAL(dp) :: ri
482
483	!
484	!-- TODO check dimensions of lat/lon and y/x match
485
486	ri = 1.0_dp / r
487
488	lat(:) = (y(:) - y0) * ri
489	lon(:) = (x(:) - x0) * ri
490	END SUBROUTINE
491
492
493	!-----------------------------------------------------------------------------!
494	! Description:
495	! ------------
496	!> Computes the reverse Plate-Carree projection of a x or y position on a
497	!> Cartesian grid.
498	!>
499	!> Input parameters:
500	!> -----------------
501	!> xy : x or y coordinate of the Cartasian grid point [m].
502	!>
503	!> xy0 : x or y coordinate that coincides with the origin of the underlying
504	!> sperical system (crossing point of the equator and prime meridian) [m].
505	!>
506	!> r : Radius of the of the underlying sphere, e.g. EARTH_RADIUS [m].
507	!>
508	!> Returns:
509	!> --------
510	!> project : Longitude (in case xy = x) or latitude (xy = y) of the given input
511	!> coordinate xy.
512	!------------------------------------------------------------------------------!
513	ELEMENTAL REAL(dp) FUNCTION project(xy, xy0, r)
514	REAL(dp), INTENT(IN) :: xy, xy0, r
515	REAL(dp) :: ri
516
517	!
518	!-- If this elemental function is called with a large array as xy, it is
519	!-- computationally more efficient to precompute the inverse radius and
520	!-- then muliply.
521	ri = 1.0_dp / r
522
523	project = (xy - xy0) * ri
524
525	END FUNCTION project
526
527
528	!------------------------------------------------------------------------------!
529	! Description:
530	! ------------
531	!> For a rotated-pole system with the origin at geographical latitude 'phi_c',
532	!> compute the geographical latitude of its rotated north pole.
533	!------------------------------------------------------------------------------!
534	REAL(dp) FUNCTION phic_to_phin(phi_c)
535	REAL(dp), INTENT(IN) :: phi_c
536
537	phic_to_phin = 0.5_dp * PI - ABS(phi_c)
538
539	END FUNCTION phic_to_phin
540
541
542	!------------------------------------------------------------------------------!
543	! Description:
544	! ------------
545	!> For a rotated-pole system with the origin at geographical latitude 'phi_c'
546	!> and longitude 'lam_c', compute the geographical longitude of its rotated
547	!> north pole.
548	!------------------------------------------------------------------------------!
549	REAL(dp) FUNCTION lamc_to_lamn(phi_c, lam_c)
550	REAL(dp), INTENT(IN) :: phi_c, lam_c
551
552	lamc_to_lamn = lam_c
553	IF (phi_c > 0.0_dp) THEN
554	lamc_to_lamn = lam_c - SIGN(PI, lam_c)
555	END IF
556
557	END FUNCTION lamc_to_lamn
558
559
560	!------------------------------------------------------------------------------!
561	! Description:
562	! ------------
563	!> Set gamma according to whether PALM domain is in the northern or southern
564	!> hemisphere of the COSMO rotated-pole system. Gamma assumes either the
565	!> value 0 or PI and is needed to work around around a bug in the
566	!> rotated-pole coordinate transformations.
567	!------------------------------------------------------------------------------!
568	REAL(dp) FUNCTION gamma_from_hemisphere(phi_cg, phi_ref)
569	REAL(dp), INTENT(IN) :: phi_cg
570	REAL(dp), INTENT(IN) :: phi_ref
571
572	LOGICAL :: palm_origin_is_south_of_cosmo_origin
573
574	palm_origin_is_south_of_cosmo_origin = (phi_cg < phi_ref)
575
576	IF (palm_origin_is_south_of_cosmo_origin) THEN
577	gamma_from_hemisphere = PI
578	ELSE
579	gamma_from_hemisphere = 0.0_dp
580	END IF
581	END FUNCTION gamma_from_hemisphere
582
583
584	!------------------------------------------------------------------------------!
585	! Description:
586	! ------------
587	!> Computes the geographical coordinates corresponding to the given rotated-pole
588	!> coordinates.
589	!>
590	!> In INIFOR, this routine is used to convert coordinates between two
591	!> rotated-pole systems: COSMO-DE's rotated-pole system, and one centred at the
592	!> PALM-4U domain centre. In this case, the PALM-4U system is thought of as the
593	!> rotated-pole system and the routine is used to rotate back to COSMO-DE's
594	!> system which is thought of as the geographical one.
595	!>
596	!> Input parameters:
597	!> -----------------
598	!> phir(:), lamr(: ): latitudes and longitudes of the rotated-pole grid
599	!>
600	!> phip, lamp: latitude and longitude of the rotated north pole
601	!>
602	!> gam: "angle between the north poles. If [gam] is not present, the other
603	!> system is the real geographical system." (original phiro2rot
604	!> description)
605	!>
606	!> Output parameters:
607	!> ------------------
608	!> phi(:,:), lam(:,:): geographical latitudes and logitudes
609	!------------------------------------------------------------------------------!
610	SUBROUTINE rotate_to_cosmo(phir, lamr, phip, lamp, phi, lam, gam)
611	REAL(dp), INTENT(IN) :: phir(0:), lamr(0:), phip, lamp, gam
612	REAL(dp), INTENT(OUT) :: phi(0:,0:), lam(0:,0:)
613
614	INTEGER :: i, j
615
616	IF ( SIZE(phi, 1) .NE. SIZE(lam, 1) .OR. &
617	SIZE(phi, 2) .NE. SIZE(lam, 2) ) THEN
618	PRINT *, "inifor: rotate_to_cosmo: Dimensions of phi and lambda do not match. Dimensions are:"
619	PRINT *, "inifor: rotate_to_cosmo: phi: ", SIZE(phi, 1), SIZE(phi, 2)
620	PRINT *, "inifor: rotate_to_cosmo: lam: ", SIZE(lam, 1), SIZE(lam, 2)
621	STOP
622	END IF
623
624	IF ( SIZE(phir) .NE. SIZE(phi, 2) .OR. &
625	SIZE(lamr) .NE. SIZE(phi, 1) ) THEN
626	PRINT *, "inifor: rotate_to_cosmo: Dimensions of phir and lamr do not match. Dimensions are:"
627	PRINT *, "inifor: rotate_to_cosmo: phir: ", SIZE(phir), SIZE(phi, 2)
628	PRINT *, "inifor: rotate_to_cosmo: lamr: ", SIZE(lamr), SIZE(phi, 1)
629	STOP
630	END IF
631
632	DO j = 0, UBOUND(phir, 1)
633	DO i = 0, UBOUND(lamr, 1)
634
635	phi(i,j) = phirot2phi(phir(j) * TO_DEGREES, &
636	lamr(i) * TO_DEGREES, &
637	phip * TO_DEGREES, &
638	gam * TO_DEGREES) * TO_RADIANS
639
640	lam(i,j) = rlarot2rla(phir(j) * TO_DEGREES, &
641	lamr(i) * TO_DEGREES, &
642	phip * TO_DEGREES, &
643	lamp * TO_DEGREES, &
644	gam * TO_DEGREES) * TO_RADIANS
645
646	END DO
647	END DO
648
649	END SUBROUTINE rotate_to_cosmo
650
651
652	!------------------------------------------------------------------------------!
653	! Description:
654	! ------------
655	!> Rotate the given vector field (x(:), y(:)) by the given 'angle'.
656	!------------------------------------------------------------------------------!
657	SUBROUTINE rotate_vector_field(x, y, angle)
658	REAL(dp), DIMENSION(:), INTENT(INOUT) :: x, y !< x and y coodrinate in arbitrary units
659	REAL(dp), INTENT(IN) :: angle !< rotation angle [deg]
660
661	INTEGER :: i
662	REAL(dp) :: sine, cosine, v_rot(2), rotation(2,2)
663
664	sine = SIN(angle * TO_RADIANS)
665	cosine = COS(angle * TO_RADIANS)
666	!
667	!-- RESAHPE() fills columns first, so the rotation matrix becomes
668	!--
669	!-- rotation = [ cosine -sine ]
670	!-- [ sine cosine ]
671	rotation = RESHAPE( (/cosine, sine, -sine, cosine/), (/2, 2/) )
672
673	DO i = LBOUND(x, 1), UBOUND(x, 1)
674
675	v_rot(:) = MATMUL(rotation, (/x(i), y(i)/))
676
677	x(i) = v_rot(1)
678	y(i) = v_rot(2)
679
680	END DO
681
682	END SUBROUTINE rotate_vector_field
683
684
685
686	!------------------------------------------------------------------------------!
687	! Description:
688	! ------------
689	!> This routine computes the local meridian convergence between a rotated-pole
690	!> and a geographical system using the Eq. (6) given in the DWD manual
691	!>
692	!> Baldauf et al. (2018), Beschreibung des operationelle KuÌrzestfrist-
693	!> vorhersagemodells COSMO-D2 und COSMO-D2-EPS und seiner Ausgabe in die
694	!> Datenbanken des DWD.
695	!> https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/cosmo_d2/cosmo_d2_dbbeschr_aktuell.pdf?__blob=publicationFile&v=2
696	!------------------------------------------------------------------------------!
697	FUNCTION meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g) &
698	RESULT(delta)
699
700	REAL(dp), INTENT(IN) :: phi_n, lam_n, phi_g, lam_g
701	REAL(dp) :: delta
702
703	delta = atan2( COS(phi_n) * SIN(lam_n - lam_g), &
704	COS(phi_g) * SIN(phi_n) - &
705	SIN(phi_g) * COS(phi_n) * COS(lam_n - lam_g) )
706
707	END FUNCTION meridian_convergence_rotated
708
709	!------------------------------------------------------------------------------!
710	! Description:
711	! ------------
712	!> Compute indices of PALM-4U grid point neighbours in the target
713	!> system (COSMO-DE) by rounding up and down. (i,j) are the indices of
714	!> the PALM-4U grid and (ii(i,j,1-4), jj(i,j,1-4)) contain the indices
715	!> of the its four neigbouring points in the COSMO-DE grid.
716	!>
717	!>
718	!> COSMO-DE grid
719	!> -------------
720	!> jj, lat
721	!> ^ j
722	!> \| \ i
723	!> jj(i,j,2/3) + ... 2 ---\--------/------ 3
724	!> \| \| ^ \ / \|
725	!> \| \| \|wp \ / \|
726	!> \| \| v \ / \|
727	!> latpos + ............ o/ (i,j) \|
728	!> \| \| : \|
729	!> \| \| :<----wl---->\|
730	!> jj(i,j,1/4) + ... 1 -------:----------- 4
731	!> \| : : :
732	!> \| : : :
733	!> \| : lonpos :
734	!> L-----+--------+------------+------> ii, lon
735	!> ii(i,j,1/2) ii(i,j,3/4)
736	!>
737	!>
738	!> Input parameters:
739	!> -----------------
740	!> source_lat, source_lon : (rotated-pole) coordinates of the source grid (e.g.
741	!> COSMO-DE)
742	!>
743	!> source_dxi, source_dyi : inverse grid spacings of the source grid.
744	!>
745	!> target_lat, target_lon : (rotated-pole) coordinates of the target grid (e.g.
746	!> COSMO-DE)
747	!>
748	!> Output parameters:
749	!> ------------------
750	!> palm_ii, palm_jj : x and y index arrays of horizontal neighbour columns
751	!>
752	!------------------------------------------------------------------------------!
753	SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon, &
754	palm_clat, palm_clon, &
755	palm_ii, palm_jj)
756
757	REAL(dp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon
758	REAL(dp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon
759	REAL(dp) :: cosmo_dxi, cosmo_dyi
760	INTEGER, DIMENSION(0:,0:,1:), INTENT(OUT) :: palm_ii, palm_jj
761
762	REAL(dp) :: lonpos, latpos, lon0, lat0
763	INTEGER :: i, j
764
765	lon0 = cosmo_lon(0)
766	lat0 = cosmo_lat(0)
767	cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
768	cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
769
770	DO j = 0, UBOUND(palm_clon, 2)!palm_grid % ny
771	DO i = 0, UBOUND(palm_clon, 1)!palm_grid % nx
772	!
773	!-- Compute the floating point index corrseponding to PALM-4U grid point
774	!-- location along target grid (COSMO-DE) axes.
775	lonpos = (palm_clon(i,j) - lon0) * cosmo_dxi
776	latpos = (palm_clat(i,j) - lat0) * cosmo_dyi
777
778	IF (lonpos < 0.0 .OR. latpos < 0.0) THEN
779	PRINT *, " Error while finding neighbours: lonpos or latpos out of bounds!"
780	PRINT *, " (i,j) = (", i, ",",j,")"
781	PRINT , " lonpos ", lonposTO_DEGREES, ", latpos ", latpos*TO_DEGREES
782	PRINT , " lon0 ", lon0 TO_DEGREES, ", lat0 ", lat0*TO_DEGREES
783	PRINT , " PALM lon ", palm_clon(i,j)TO_DEGREES, ", PALM lat ",palm_clat(i,j)*TO_DEGREES
784	STOP
785	END IF
786
787	palm_ii(i,j,1) = FLOOR(lonpos)
788	palm_ii(i,j,2) = FLOOR(lonpos)
789	palm_ii(i,j,3) = CEILING(lonpos)
790	palm_ii(i,j,4) = CEILING(lonpos)
791
792	palm_jj(i,j,1) = FLOOR(latpos)
793	palm_jj(i,j,2) = CEILING(latpos)
794	palm_jj(i,j,3) = CEILING(latpos)
795	palm_jj(i,j,4) = FLOOR(latpos)
796	END DO
797	END DO
798
799	END SUBROUTINE find_horizontal_neighbours
800
801
802	!------------------------------------------------------------------------------!
803	! Description:
804	! ------------
805	!> Computes linear vertical interpolation neighbour indices and weights for each
806	!> column of the given palm grid.
807	!------------------------------------------------------------------------------!
808	SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid, &
809	palm_intermediate )
810	TYPE(grid_definition), INTENT(INOUT) :: palm_grid
811	TYPE(grid_definition), INTENT(IN) :: palm_intermediate
812
813	INTEGER :: i, j, k, nx, ny, nz, nlev, k_intermediate
814	LOGICAL :: point_is_below_grid, point_is_above_grid, &
815	point_is_in_current_cell
816	REAL(dp) :: current_height, column_base, column_top, h_top, h_bottom, &
817	weight
818
819	nx = palm_grid % nx
820	ny = palm_grid % ny
821	nz = palm_grid % nz
822	nlev = palm_intermediate % nz
823
824	!
825	!-- in each column of the fine grid, find vertical neighbours of every cell
826	DO j = 0, ny
827	DO i = 0, nx
828
829	k_intermediate = 0
830
831	column_base = palm_intermediate % h(i,j,0)
832	column_top = palm_intermediate % h(i,j,nlev)
833
834	!
835	!-- scan through palm_grid column and set neighbour indices in
836	!-- case current_height is either below column_base, in the current
837	!-- cell, or above column_top. Keep increasing current cell index until
838	!-- the current cell overlaps with the current_height.
839	DO k = 1, nz
840
841	!
842	!-- Memorize the top and bottom boundaries of the coarse cell and the
843	!-- current height within it
844	current_height = palm_grid % z(k) + palm_grid % z0
845	h_top = palm_intermediate % h(i,j,k_intermediate+1)
846	h_bottom = palm_intermediate % h(i,j,k_intermediate)
847
848	point_is_above_grid = (current_height > column_top) !22000m, very unlikely
849	point_is_below_grid = (current_height < column_base)
850
851	point_is_in_current_cell = ( &
852	current_height >= h_bottom .AND. &
853	current_height < h_top &
854	)
855
856	!
857	!-- set default weights
858	palm_grid % w_verti(i,j,k,1:2) = 0.0_dp
859
860	IF (point_is_above_grid) THEN
861
862	palm_grid % kk(i,j,k,1:2) = nlev
863	palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp
864
865	message = "PALM-4U grid extends above COSMO-DE model top."
866	CALL inifor_abort('find_vertical_neighbours_and_weights', message)
867
868	ELSE IF (point_is_below_grid) THEN
869
870	palm_grid % kk(i,j,k,1:2) = 0
871	palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp
872
873	ELSE
874	!
875	!-- cycle through intermediate levels until current
876	!-- intermediate-grid cell overlaps with current_height
877	DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
878	k_intermediate = k_intermediate + 1
879
880	h_top = palm_intermediate % h(i,j,k_intermediate+1)
881	h_bottom = palm_intermediate % h(i,j,k_intermediate)
882	point_is_in_current_cell = ( &
883	current_height >= h_bottom .AND. &
884	current_height < h_top &
885	)
886	END DO
887
888	IF (k_intermediate > nlev-1) THEN
889	message = "Index " // TRIM(str(k_intermediate)) // &
890	" is above intermediate grid range."
891	CALL inifor_abort('find_vertical_neighbours', message)
892	END IF
893
894	palm_grid % kk(i,j,k,1) = k_intermediate
895	palm_grid % kk(i,j,k,2) = k_intermediate + 1
896
897	!
898	!-- compute vertical weights
899	weight = (h_top - current_height) / (h_top - h_bottom)
900	palm_grid % w_verti(i,j,k,1) = weight
901	palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight
902	END IF
903
904	END DO
905
906	END DO
907	END DO
908
909	END SUBROUTINE find_vertical_neighbours_and_weights_interp
910
911
912	!------------------------------------------------------------------------------!
913	! Description:
914	! ------------
915	!> Computes linear vertical interpolation neighbour indices and weights for each
916	!> column of the given averaging grid.
917	!>
918	!> The difference to the _interp variant of this routine lies in how columns
919	!> are adressed. While the _interp variant loops over all PALM grid columns
920	!> given by combinations of all index indices (i,j), this routine loops over a
921	!> subset of COSMO columns, which are sequantlially stored in the index lists
922	!> iii(:) and jjj(:).
923	!------------------------------------------------------------------------------!
924	SUBROUTINE find_vertical_neighbours_and_weights_average( &
925	avg_grid, level_based_averaging &
926	)
927
928	TYPE(grid_definition), INTENT(INOUT), TARGET :: avg_grid
929	LOGICAL :: level_based_averaging
930
931	INTEGER :: i, j, k_palm, k_intermediate, l, nlev
932	LOGICAL :: point_is_below_grid, point_is_above_grid, &
933	point_is_in_current_cell
934	REAL(dp) :: current_height, column_base, column_top, h_top, &
935	h_bottom, weight
936	REAL(dp), POINTER :: cosmo_h(:,:,:)
937
938
939	avg_grid % k_min = LBOUND(avg_grid % z, 1)
940
941	nlev = SIZE(avg_grid % cosmo_h, 3)
942
943	IF (level_based_averaging) THEN
944	cosmo_h => avg_grid % h
945	ELSE
946	cosmo_h => avg_grid % cosmo_h
947	END IF
948
949	!
950	!-- in each column of the fine grid, find vertical neighbours of every cell
951	DO l = 1, avg_grid % n_columns
952
953	IF (level_based_averaging) THEN
954	i = 1
955	j = 1
956	ELSE
957	i = avg_grid % iii(l)
958	j = avg_grid % jjj(l)
959	END IF
960
961	column_base = cosmo_h(i,j,1)
962	column_top = cosmo_h(i,j,nlev)
963
964	!
965	!-- scan through avg_grid column until and set neighbour indices in
966	!-- case current_height is either below column_base, in the current
967	!-- cell, or above column_top. Keep increasing current cell index until
968	!-- the current cell overlaps with the current_height.
969	k_intermediate = 1 !avg_grid % cosmo_h is indezed 1-based.
970	DO k_palm = 1, avg_grid % nz
971
972	!
973	!-- Memorize the top and bottom boundaries of the coarse cell and the
974	!-- current height within it
975	current_height = avg_grid % z(k_palm) + avg_grid % z0
976	h_top = cosmo_h(i,j,k_intermediate+1)
977	h_bottom = cosmo_h(i,j,k_intermediate)
978
979	!
980	!-- COSMO column top is located at 22000m, point_is_above_grid is very
981	!-- unlikely.
982	point_is_above_grid = (current_height > column_top)
983	point_is_below_grid = (current_height < column_base)
984
985	point_is_in_current_cell = ( &
986	current_height >= h_bottom .AND. &
987	current_height < h_top &
988	)
989
990	!
991	!-- set default weights
992	avg_grid % w(l,k_palm,1:2) = 0.0_dp
993
994	IF (point_is_above_grid) THEN
995
996	avg_grid % kkk(l,k_palm,1:2) = nlev
997	avg_grid % w(l,k_palm,1:2) = - 2.0_dp
998
999	message = "PALM-4U grid extends above COSMO-DE model top."
1000	CALL inifor_abort('find_vertical_neighbours_and_weights_average', message)
1001
1002	ELSE IF (point_is_below_grid) THEN
1003
1004	avg_grid % kkk(l,k_palm,1:2) = 0
1005	avg_grid % w(l,k_palm,1:2) = - 2.0_dp
1006	avg_grid % k_min = MAX(k_palm + 1, avg_grid % k_min)
1007	ELSE
1008	!
1009	!-- cycle through intermediate levels until current
1010	!-- intermediate-grid cell overlaps with current_height
1011	DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
1012	k_intermediate = k_intermediate + 1
1013
1014	h_top = cosmo_h(i,j,k_intermediate+1)
1015	h_bottom = cosmo_h(i,j,k_intermediate)
1016	point_is_in_current_cell = ( &
1017	current_height >= h_bottom .AND. &
1018	current_height < h_top &
1019	)
1020	END DO
1021
1022	!
1023	!-- k_intermediate = 48 indicates the last section (indices 48 and 49), i.e.
1024	!-- k_intermediate = 49 is not the beginning of a valid cell.
1025	IF (k_intermediate > nlev-1) THEN
1026	message = "Index " // TRIM(str(k_intermediate)) // &
1027	" is above intermediate grid range."
1028	CALL inifor_abort('find_vertical_neighbours', message)
1029	END IF
1030
1031	avg_grid % kkk(l,k_palm,1) = k_intermediate
1032	avg_grid % kkk(l,k_palm,2) = k_intermediate + 1
1033
1034	!
1035	!-- compute vertical weights
1036	weight = (h_top - current_height) / (h_top - h_bottom)
1037	avg_grid % w(l,k_palm,1) = weight
1038	avg_grid % w(l,k_palm,2) = 1.0_dp - weight
1039	END IF
1040
1041	!
1042	!-- Loop over PALM levels k
1043	END DO
1044
1045	!
1046	!-- Loop over averaging columns l
1047	END DO
1048
1049	END SUBROUTINE find_vertical_neighbours_and_weights_average
1050
1051	!------------------------------------------------------------------------------!
1052	! Description:
1053	! ------------
1054	!> Compute the four weights for horizontal bilinear interpolation given the
1055	!> coordinates clon(i,j) clat(i,j) of the PALM-4U grid in the COSMO-DE
1056	!> rotated-pole grid and the neightbour indices ii(i,j,1-4) and jj(i,j,1-4).
1057	!>
1058	!> Input parameters:
1059	!> -----------------
1060	!> palm_grid % clon : longitudes of PALM-4U scalars (cell centres) in COSMO-DE's rotated-pole grid [rad]
1061	!>
1062	!> palm_grid % clat : latitudes of PALM-4U cell centres in COSMO-DE's rotated-pole grid [rad]
1063	!>
1064	!> cosmo_grid % lon : rotated-pole longitudes of scalars (cell centres) of the COSMO-DE grid [rad]
1065	!>
1066	!> cosmo_grid % lat : rotated-pole latitudes of scalars (cell centers) of the COSMO-DE grid [rad]
1067	!>
1068	!> cosmo_grid % dxi : inverse grid spacing in the first dimension [m^-1]
1069	!>
1070	!> cosmo_grid % dyi : inverse grid spacing in the second dimension [m^-1]
1071	!>
1072	!> Output parameters:
1073	!> ------------------
1074	!> palm_grid % w_horiz(:,:,1-4) : weights for bilinear horizontal interpolation
1075	!
1076	! COSMO-DE grid
1077	! -------------
1078	! jj, lat
1079	! ^ j
1080	! \| \ i
1081	! jj(i,j,2/3) + ... 2 ---\--------/------ 3
1082	! \| \| ^ \ / \|
1083	! \| \| \|wp \ / \|
1084	! \| \| v \ / \|
1085	! latpos + ............ o/ (i,j) \|
1086	! \| \| : \|
1087	! \| \| :<----wl---->\|
1088	! jj(i,j,1/4) + ... 1 -------:----------- 4
1089	! \| : : :
1090	! \| : : :
1091	! \| : lonpos :
1092	! L-----+--------+------------+------> ii, lon
1093	! ii(i,j,1/2) ii(i,j,3/4)
1094	!
1095	!------------------------------------------------------------------------------!
1096	SUBROUTINE compute_horizontal_interp_weights(cosmo_lat, cosmo_lon, &
1097	palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)
1098
1099	REAL(dp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon
1100	REAL(dp) :: cosmo_dxi, cosmo_dyi
1101	REAL(dp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon
1102	INTEGER, DIMENSION(0:,0:,1:), INTENT(IN) :: palm_ii, palm_jj
1103
1104	REAL(dp), DIMENSION(0:,0:,1:), INTENT(OUT) :: palm_w_horiz
1105
1106	REAL(dp) :: wl, wp
1107	INTEGER :: i, j
1108
1109	cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
1110	cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
1111
1112	DO j = 0, UBOUND(palm_clon, 2)
1113	DO i = 0, UBOUND(palm_clon, 1)
1114
1115	!
1116	!-- weight in lambda direction
1117	wl = ( cosmo_lon(palm_ii(i,j,4)) - palm_clon(i,j) ) * cosmo_dxi
1118
1119	!
1120	!-- weight in phi direction
1121	wp = ( cosmo_lat(palm_jj(i,j,2)) - palm_clat(i,j) ) * cosmo_dyi
1122
1123	IF (wl > 1.0_dp .OR. wl < 0.0_dp) THEN
1124	message = "Horizontal weight wl = " // TRIM(real_to_str(wl)) // &
1125	" is out bounds."
1126	CALL inifor_abort('compute_horizontal_interp_weights', message)
1127	END IF
1128	IF (wp > 1.0_dp .OR. wp < 0.0_dp) THEN
1129	message = "Horizontal weight wp = " // TRIM(real_to_str(wp)) // &
1130	" is out bounds."
1131	CALL inifor_abort('compute_horizontal_interp_weights', message)
1132	END IF
1133
1134	palm_w_horiz(i,j,1) = wl * wp
1135	palm_w_horiz(i,j,2) = wl * (1.0_dp - wp)
1136	palm_w_horiz(i,j,3) = (1.0_dp - wl) * (1.0_dp - wp)
1137	palm_w_horiz(i,j,4) = 1.0_dp - SUM( palm_w_horiz(i,j,1:3) )
1138
1139	END DO
1140	END DO
1141
1142	END SUBROUTINE compute_horizontal_interp_weights
1143
1144
1145	!------------------------------------------------------------------------------!
1146	! Description:
1147	! ------------
1148	!> Interpolates u and v components of velocities located at cell faces to the
1149	!> cell centres by averaging neighbouring values.
1150	!>
1151	!> This routine is designed to be used with COSMO-DE arrays where there are the
1152	!> same number of grid points for scalars (centres) and velocities (faces). In
1153	!> COSMO-DE the velocity points are staggared one half grid spaceing up-grid
1154	!> which means the first centre point has to be omitted and is set to zero.
1155	!------------------------------------------------------------------------------!
1156	SUBROUTINE centre_velocities(u_face, v_face, u_centre, v_centre)
1157	REAL(dp), DIMENSION(0:,0:,0:), INTENT(IN) :: u_face, v_face
1158	REAL(dp), DIMENSION(0:,0:,0:), INTENT(OUT) :: u_centre, v_centre
1159	INTEGER :: nx, ny
1160
1161	nx = UBOUND(u_face, 1)
1162	ny = UBOUND(u_face, 2)
1163
1164	u_centre(0,:,:) = 0.0_dp
1165	u_centre(1:,:,:) = 0.5_dp * ( u_face(0:nx-1,:,:) + u_face(1:,:,:) )
1166
1167	v_centre(:,0,:) = 0.0_dp
1168	v_centre(:,1:,:) = 0.5_dp * ( v_face(:,0:ny-1,:) + v_face(:,1:,:) )
1169	END SUBROUTINE centre_velocities
1170
1171
1172	!------------------------------------------------------------------------------!
1173	! Description:
1174	! ------------
1175	!> Compute the geographical latitude of a point given in rotated-pole cordinates
1176	!------------------------------------------------------------------------------!
1177	FUNCTION phirot2phi (phirot, rlarot, polphi, polgam)
1178
1179	REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole
1180	REAL(dp), INTENT (IN) :: phirot !< latitude in the rotated system
1181	REAL(dp), INTENT (IN) :: rlarot !< longitude in the rotated system
1182	REAL(dp), INTENT (IN) :: polgam !< angle between the north poles of the systems
1183
1184	REAL(dp) :: phirot2phi !< latitude in the geographical system
1185
1186	REAL(dp) :: zsinpol, zcospol, zphis, zrlas, zarg, zgam
1187
1188	zsinpol = SIN(polphi * TO_RADIANS)
1189	zcospol = COS(polphi * TO_RADIANS)
1190	zphis = phirot * TO_RADIANS
1191
1192	IF (rlarot > 180.0_dp) THEN
1193	zrlas = rlarot - 360.0_dp
1194	ELSE
1195	zrlas = rlarot
1196	END IF
1197	zrlas = zrlas * TO_RADIANS
1198
1199	IF (polgam /= 0.0_dp) THEN
1200	zgam = polgam * TO_RADIANS
1201	zarg = zsinpol * SIN (zphis) + &
1202	zcospol * COS(zphis) * ( COS(zrlas) * COS(zgam) - &
1203	SIN(zgam) * SIN(zrlas) )
1204	ELSE
1205	zarg = zcospol * COS (zphis) * COS (zrlas) + zsinpol * SIN (zphis)
1206	END IF
1207
1208	phirot2phi = ASIN (zarg) * TO_DEGREES
1209
1210	END FUNCTION phirot2phi
1211
1212
1213	!------------------------------------------------------------------------------!
1214	! Description:
1215	! ------------
1216	!> Compute the geographical latitude of a point given in rotated-pole cordinates
1217	!------------------------------------------------------------------------------!
1218	FUNCTION phi2phirot (phi, rla, polphi, pollam)
1219
1220	REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole
1221	REAL(dp), INTENT (IN) :: pollam !< longitude of the rotated north pole
1222	REAL(dp), INTENT (IN) :: phi !< latitude in the geographical system
1223	REAL(dp), INTENT (IN) :: rla !< longitude in the geographical system
1224
1225	REAL(dp) :: phi2phirot !< longitude in the rotated system
1226
1227	REAL(dp) :: zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1
1228
1229	zsinpol = SIN(polphi * TO_RADIANS)
1230	zcospol = COS(polphi * TO_RADIANS)
1231	zlampol = pollam * TO_RADIANS
1232	zphi = phi * TO_RADIANS
1233
1234	IF (rla > 180.0_dp) THEN
1235	zrla1 = rla - 360.0_dp
1236	ELSE
1237	zrla1 = rla
1238	END IF
1239	zrla = zrla1 * TO_RADIANS
1240
1241	zarg1 = SIN(zphi) * zsinpol
1242	zarg2 = COS(zphi) * zcospol * COS(zrla - zlampol)
1243
1244	phi2phirot = ASIN(zarg1 + zarg2) * TO_DEGREES
1245
1246	END FUNCTION phi2phirot
1247
1248
1249	!------------------------------------------------------------------------------!
1250	! Description:
1251	! ------------
1252	!> Compute the geographical longitude of a point given in rotated-pole cordinates
1253	!------------------------------------------------------------------------------!
1254	FUNCTION rlarot2rla(phirot, rlarot, polphi, pollam, polgam)
1255
1256	REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole
1257	REAL(dp), INTENT (IN) :: pollam !< longitude of the rotated north pole
1258	REAL(dp), INTENT (IN) :: phirot !< latitude in the rotated system
1259	REAL(dp), INTENT (IN) :: rlarot !< longitude in the rotated system
1260	REAL(dp), INTENT (IN) :: polgam !< angle between the north poles of the systems
1261
1262	REAL(dp) :: rlarot2rla !< latitude in the geographical system
1263
1264	REAL(dp) :: zsinpol, zcospol, zlampol, zphis, zrlas, zarg1, zarg2, zgam
1265
1266	zsinpol = SIN(TO_RADIANS * polphi)
1267	zcospol = COS(TO_RADIANS * polphi)
1268	zlampol = TO_RADIANS * pollam
1269	zphis = TO_RADIANS * phirot
1270
1271	IF (rlarot > 180.0_dp) THEN
1272	zrlas = rlarot - 360.0_dp
1273	ELSE
1274	zrlas = rlarot
1275	END IF
1276	zrlas = TO_RADIANS * zrlas
1277
1278	IF (polgam /= 0.0_dp) THEN
1279	zgam = TO_RADIANS * polgam
1280	zarg1 = SIN(zlampol) * (zcospol * SIN(zphis) - zsinpolCOS(zphis) &
1281	(COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) - &
1282	COS(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) + &
1283	COS(zrlas) * SIN(zgam) )
1284
1285	zarg2 = COS (zlampol) * (zcospol * SIN(zphis) - zsinpolCOS(zphis) &
1286	(COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) + &
1287	SIN(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) + &
1288	COS(zrlas) * SIN(zgam) )
1289	ELSE
1290	zarg1 = SIN (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis) + &
1291	zcospol * SIN(zphis)) - &
1292	COS (zlampol) * SIN(zrlas) * COS(zphis)
1293	zarg2 = COS (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis) + &
1294	zcospol * SIN(zphis)) + &
1295	SIN (zlampol) * SIN(zrlas) * COS(zphis)
1296	END IF
1297
1298	IF (zarg2 == 0.0_dp) zarg2 = 1.0E-20_dp
1299
1300	rlarot2rla = ATAN2(zarg1,zarg2) * TO_DEGREES
1301
1302	END FUNCTION rlarot2rla
1303
1304
1305	!------------------------------------------------------------------------------!
1306	! Description:
1307	! ------------
1308	!> Compute the rotated-pole longitude of a point given in geographical cordinates
1309	!------------------------------------------------------------------------------!
1310	FUNCTION rla2rlarot ( phi, rla, polphi, pollam, polgam )
1311
1312	REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole
1313	REAL(dp), INTENT (IN) :: pollam !< longitude of the rotated north pole
1314	REAL(dp), INTENT (IN) :: phi !< latitude in geographical system
1315	REAL(dp), INTENT (IN) :: rla !< longitude in geographical system
1316	REAL(dp), INTENT (IN) :: polgam !< angle between the north poles of the systems
1317
1318	REAL (KIND=dp) :: rla2rlarot !< latitude in the the rotated system
1319
1320	REAL (KIND=dp) :: zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1
1321
1322	zsinpol = SIN(polphi * TO_RADIANS)
1323	zcospol = COS(polphi * TO_RADIANS)
1324	zlampol = pollam * TO_RADIANS
1325	zphi = phi * TO_RADIANS
1326
1327	IF (rla > 180.0_dp) THEN
1328	zrla1 = rla - 360.0_dp
1329	ELSE
1330	zrla1 = rla
1331	END IF
1332	zrla = zrla1 * TO_RADIANS
1333
1334	zarg1 = - SIN (zrla-zlampol) * COS(zphi)
1335	zarg2 = - zsinpol * COS(zphi) * COS(zrla-zlampol) + zcospol * SIN(zphi)
1336
1337	IF (zarg2 == 0.0_dp) zarg2 = 1.0E-20_dp
1338
1339	rla2rlarot = ATAN2 (zarg1,zarg2) * TO_DEGREES
1340
1341	IF (polgam /= 0.0_dp ) THEN
1342	rla2rlarot = polgam + rla2rlarot
1343	IF (rla2rlarot > 180._dp) rla2rlarot = rla2rlarot - 360.0_dp
1344	END IF
1345
1346	END FUNCTION rla2rlarot
1347
1348
1349	!------------------------------------------------------------------------------!
1350	! Description:
1351	! ------------
1352	!> Rotate the given velocity vector (u,v) from the geographical to the
1353	!> rotated-pole system
1354	!------------------------------------------------------------------------------!
1355	SUBROUTINE uv2uvrot(u, v, rlat, rlon, pollat, pollon, urot, vrot)
1356
1357	REAL(dp), INTENT (IN) :: u, v !< wind components in the true geographical system
1358	REAL(dp), INTENT (IN) :: rlat, rlon !< coordinates in the true geographical system
1359	REAL(dp), INTENT (IN) :: pollat, pollon !< latitude and longitude of the north pole of the rotated grid
1360
1361	REAL(dp), INTENT (OUT) :: urot, vrot !< wind components in the rotated grid
1362
1363	REAL (dp) :: zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm
1364
1365	zsinpol = SIN(pollat * TO_RADIANS)
1366	zcospol = COS(pollat * TO_RADIANS)
1367	zlonp = (pollon-rlon) * TO_RADIANS
1368	zlat = rlat * TO_RADIANS
1369
1370	zarg1 = zcospol * SIN(zlonp)
1371	zarg2 = zsinpol * COS(zlat) - zcospol * SIN(zlat) * COS(zlonp)
1372	znorm = 1.0_dp / SQRT(zarg1zarg1 + zarg2zarg2)
1373
1374	urot = u * zarg2 * znorm - v * zarg1 * znorm
1375	vrot = u * zarg1 * znorm + v * zarg2 * znorm
1376
1377	END SUBROUTINE uv2uvrot
1378
1379
1380	!------------------------------------------------------------------------------!
1381	! Description:
1382	! ------------
1383	!> Rotate the given velocity vector (urot, vrot) from the rotated-pole to the
1384	!> geographical system
1385	!------------------------------------------------------------------------------!
1386	SUBROUTINE uvrot2uv (urot, vrot, rlat, rlon, pollat, pollon, u, v)
1387
1388	REAL(dp), INTENT(IN) :: urot, vrot !< wind components in the rotated grid
1389	REAL(dp), INTENT(IN) :: rlat, rlon !< latitude and longitude in the true geographical system
1390	REAL(dp), INTENT(IN) :: pollat, pollon !< latitude and longitude of the north pole of the rotated grid
1391
1392	REAL(dp), INTENT(OUT) :: u, v !< wind components in the true geographical system
1393
1394	REAL(dp) :: zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm
1395
1396	zsinpol = SIN(pollat * TO_RADIANS)
1397	zcospol = COS(pollat * TO_RADIANS)
1398	zlonp = (pollon-rlon) * TO_RADIANS
1399	zlat = rlat * TO_RADIANS
1400
1401	zarg1 = zcospol * SIN(zlonp)
1402	zarg2 = zsinpol * COS(zlat) - zcospol * SIN(zlat) * COS(zlonp)
1403	znorm = 1.0_dp / SQRT(zarg1zarg1 + zarg2zarg2)
1404
1405	u = urot * zarg2 * znorm + vrot * zarg1 * znorm
1406	v = - urot * zarg1 * znorm + vrot * zarg2 * znorm
1407
1408	END SUBROUTINE uvrot2uv
1409
1410	END MODULE inifor_transform
1411

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |