Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |