Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |