Home

Context Navigation

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

Last change on this file since 4641 was 4553, checked in by eckhard, 4 years ago
Fixed domain extend check, readablity and documentation improvements
Property svn:keywords set to `Id`
File size: 59.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |