Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |