Home

Context Navigation

transform.f90 @ 3398

Last change on this file since 3398 was 3395, checked in by eckhard, 6 years ago
inifor: Added computation of geostrophic winds from COSMO input
Property svn:keywords set to `Id`
File size: 46.3 KB

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

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

Context Navigation

source: palm/trunk/UTIL/inifor/src/transform.f90 @ 3398

Download in other formats: