Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |