Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |