Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |