Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |