Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |