Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |