Home

Context Navigation

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

Last change on this file since 4834 was 4808, checked in by eckhard, 4 years ago

Small error handling improvements

Added offending height values to model-top error message
Temporarily disabling dataset validation, enabling compatibility for HHL fields without time dimension

Property svn:keywords set to Id

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |