Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |