Home

Context Navigation

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

Last change on this file since 3614 was 3614, checked in by raasch, 6 years ago
unused variables removed, abort renamed inifor_abort to avoid intrinsic problem in Fortran
Property svn:keywords set to `Id`
File size: 53.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |