Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |