Home

Context Navigation

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

Last change on this file since 3615 was 3615, checked in by raasch, 5 years ago
bugfix for last commit: abort replaced by inifor_abort
Property svn:keywords set to `Id`
File size: 53.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |