Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |