Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |