Home

Context Navigation

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

Last change on this file since 3601 was 3557, checked in by eckhard, 6 years ago
inifor: Updated documentation
Property svn:keywords set to `Id`
File size: 52.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |