Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |