Home

Context Navigation

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

Last change on this file since 3318 was 3183, checked in by suehring, 6 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 36.2 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 3183 2018-07-27 14:25:55Z sward $
[3182]	28	! Introduced new PALM grid stretching
	29	! Removed unnecessary subroutine parameters
	30	! Renamed kcur to k_intermediate
[2696]	31	!
	32	!
[3183]	33	! 3182 2018-07-27 13:36:03Z suehring
[2696]	34	! Initial revision
	35	!
	36	!
	37	!
	38	! Authors:
	39	! --------
	40	! @author Eckhard Kadasch
	41	!
	42	! Description:
	43	! ------------
	44	!> The transform module provides INIFOR's low-level transformation and
	45	!> interpolation routines. The rotated-pole transformation routines phirot2phi,
	46	!> phi2phirot, rlarot2rla, rla2rlarot, uv2uvrot, and uvrot2uv are adapted from
	47	!> int2lm's utility routines.
	48	!------------------------------------------------------------------------------!
	49	MODULE transform
	50
	51	USE control
	52	USE defs, &
	53	ONLY: TO_DEGREES, TO_RADIANS, PI, dp
	54	USE types
	55	USE util, &
	56	ONLY: real_to_str, str
	57
	58	IMPLICIT NONE
	59
	60	CONTAINS
	61
	62	!------------------------------------------------------------------------------!
	63	! Description:
	64	! ------------
	65	!> Interpolates linearly in the vertical direction in very column (i,j) of the
	66	!> output array outvar(i,j,:) using values of the source array invar. In cells
	67	!> that are outside the COSMO-DE domain, indicated by negative interpolation
	68	!> weights, extrapolate constantly from the cell above.
	69	!>
	70	!> Input parameters:
	71	!> -----------------
	72	!> invar : Array of source data
	73	!>
	74	!> outgrid % kk : Array of vertical neighbour indices. kk(i,j,k,:) contain the
	75	!> indices of the two vertical neighbors of PALM-4U point (i,j,k) on the
	76	!> input grid corresponding to the source data invar.
	77	!>
	78	!> outgrid % w_verti : Array of weights for vertical linear interpolation
	79	!> corresponding to neighbour points indexed by kk.
	80	!>
	81	!> Output papameters:
	82	!> ------------------
	83	!> outvar : Array of interpolated data
	84	!------------------------------------------------------------------------------!
	85	SUBROUTINE interpolate_1d(in_arr, out_arr, outgrid)
	86	TYPE(grid_definition), INTENT(IN) :: outgrid
	87	REAL(dp), INTENT(IN) :: in_arr(0:,0:,0:)
[3182]	88	REAL(dp), INTENT(OUT) :: out_arr(0:,0:,:)
[2696]	89
[3182]	90	INTEGER :: i, j, k, l, nz
[2696]	91
	92	nz = UBOUND(out_arr, 3)
	93
[3182]	94	DO j = LBOUND(out_arr, 2), UBOUND(out_arr, 2)
	95	DO i = LBOUND(out_arr, 1), UBOUND(out_arr, 1)
	96	DO k = nz, LBOUND(out_arr, 3), -1
[2696]	97
	98	! TODO: Remove IF clause and extrapolate based on a critical vertical
	99	! TODO: index marking the lower bound of COSMO-DE data coverage.
	100	! Check for negative interpolation weights indicating grid points
	101	! below COSMO-DE domain and extrapolate from the top in such cells.
	102	IF (outgrid % w_verti(i,j,k,1) < -1.0_dp .AND. k < nz) THEN
	103	out_arr(i,j,k) = out_arr(i,j,k+1)
	104	ELSE
	105	out_arr(i,j,k) = 0.0_dp
	106	DO l = 1, 2
[3182]	107	out_arr(i,j,k) = out_arr(i,j,k) + &
	108	outgrid % w_verti(i,j,k,l) * &
[2696]	109	in_arr(i,j,outgrid % kk(i,j,k, l) )
	110	END DO
	111	END IF
	112	END DO
	113	END DO
	114	END DO
	115	END SUBROUTINE interpolate_1d
	116
	117
	118	!------------------------------------------------------------------------------!
	119	! Description:
	120	! ------------
	121	!> Interpolates bi-linearly in horizontal planes on every k level of the output
	122	!> array outvar(:,:,k) using values of the source array invar(:,:,:). The source
	123	!> (invar) and interpolation array (outvar) need to have matching dimensions.
	124	!>
	125	!> Input parameters:
	126	!> -----------------
	127	!> invar : Array of source data
	128	!>
	129	!> outgrid % ii, % jj : Array of neighbour indices in x and y direction.
	130	!> ii(i,j,k,:), and jj(i,j,k,:) contain the four horizontal neighbour points
	131	!> of PALM-4U point (i,j,k) on the input grid corresponding to the source
	132	!> data invar. (The outgrid carries the relationship with the ingrid in the
	133	! form of the interpoaltion weights.)
	134	!>
	135	!> outgrid % w_horiz: Array of weights for horizontal bi-linear interpolation
	136	!> corresponding to neighbour points indexed by ii and jj.
	137	!>
	138	!> Output papameters:
	139	!> ------------------
	140	!> outvar : Array of interpolated data
	141	!------------------------------------------------------------------------------!
	142	SUBROUTINE interpolate_2d(invar, outvar, outgrid, ncvar)
	143	! I index 0-based for the indices of the outvar to be consistent with the
	144	! outgrid indices and interpolation weights.
[3182]	145	TYPE(grid_definition), INTENT(IN) :: outgrid
	146	REAL(dp), INTENT(IN) :: invar(0:,0:,0:)
	147	REAL(dp), INTENT(OUT) :: outvar(0:,0:,0:)
[2696]	148	TYPE(nc_var), INTENT(IN), OPTIONAL :: ncvar
	149
	150	INTEGER :: i, j, k, l
	151
	152	! TODO: check if input dimensions are consistent, i.e. ranges are correct
	153	IF (UBOUND(outvar, 3) .GT. UBOUND(invar, 3)) THEN
	154	message = "Output array for '" // TRIM(ncvar % name) // "' has ' more levels (" // &
	155	TRIM(str(UBOUND(outvar, 3))) // ") than input variable ("//&
	156	TRIM(str(UBOUND(invar, 3))) // ")."
	157	CALL abort('interpolate_2d', message)
	158	END IF
	159
	160	DO k = 0, UBOUND(outvar, 3)
	161	DO j = 0, UBOUND(outvar, 2)
	162	DO i = 0, UBOUND(outvar, 1)
	163	outvar(i,j,k) = 0.0_dp
	164	DO l = 1, 4
	165
	166	outvar(i,j,k) = outvar(i,j,k) + &
	167	outgrid % w_horiz(i,j,l) * invar( outgrid % ii(i,j,l), &
	168	outgrid % jj(i,j,l), &
	169	k )
	170	END DO
	171	END DO
	172	END DO
	173	END DO
	174
	175	END SUBROUTINE interpolate_2d
	176
	177
	178	SUBROUTINE average_2d(in_arr, out_arr, imin, imax, jmin, jmax)
	179	REAL(dp), INTENT(IN) :: in_arr(0:,0:,0:)
	180	REAL(dp), INTENT(OUT) :: out_arr(0:)
	181	INTEGER, INTENT(IN) :: imin, imax, jmin, jmax
	182
	183	INTEGER :: i, j, k
	184	REAL(dp) :: ni
	185
	186	IF (imin < 0) CALL abort('average_2d', "imin < 0.")
	187	IF (jmin < 0) CALL abort('average_2d', "jmin < 0.")
	188	IF (imax > UBOUND(in_arr, 1)) CALL abort('average_2d', "imax out of i bound.")
	189	IF (jmax > UBOUND(in_arr, 2)) CALL abort('average_2d', "jmax out of j bound.")
	190
	191	DO k = 0, UBOUND(out_arr, 1)
	192
	193	out_arr(k) = 0.0_dp
	194	DO j = jmin, jmax
	195	DO i = imin, imax
	196	out_arr(k) = out_arr(k) + in_arr(i, j, k)
	197	END DO
	198	END DO
	199
	200	END DO
	201
	202	! devide by number of grid points
	203	ni = 1.0_dp / ( (imax - imin + 1) * (jmax - jmin + 1) )
	204	out_arr(:) = out_arr(:) * ni
	205
	206	END SUBROUTINE average_2d
	207
	208
	209	SUBROUTINE interpolate_3d(source_array, palm_array, palm_intermediate, palm_grid)
	210	TYPE(grid_definition), INTENT(IN) :: palm_intermediate, palm_grid
	211	REAL(dp), DIMENSION(:,:,:), INTENT(IN) :: source_array
	212	REAL(dp), DIMENSION(:,:,:), INTENT(OUT) :: palm_array
	213	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: intermediate_array
	214	INTEGER :: nx, ny, nz
	215
	216	nx = palm_intermediate % nx
	217	ny = palm_intermediate % ny
	218	nz = palm_intermediate % nz ! nlev
	219
	220	! Interpolate from COSMO-DE to intermediate grid. Allocating with one
	221	! less point in the vertical, since scalars like T have 50 instead of 51
	222	! points in COSMO-DE.
	223	ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nz-1)) !
	224
	225	CALL interpolate_2d(source_array, intermediate_array, palm_intermediate)
	226
	227	! Interpolate from intermediate grid to palm_grid grid, includes
	228	! extrapolation for cells below COSMO-DE domain.
	229	CALL interpolate_1d(intermediate_array, palm_array, palm_grid)
	230
	231	DEALLOCATE(intermediate_array)
	232
	233	END SUBROUTINE interpolate_3d
	234
	235
	236	SUBROUTINE average_profile(source_array, profile_array, imin, imax, jmin, jmax,&
	237	palm_intermediate, palm_grid)
	238	TYPE(grid_definition), INTENT(IN) :: palm_intermediate, palm_grid
	239	REAL(dp), DIMENSION(:,:,:), INTENT(IN) :: source_array
	240	INTEGER, INTENT(IN) :: imin, imax, jmin, jmax
	241	REAL(dp), DIMENSION(0:,0:,0:), INTENT(OUT) :: profile_array
	242	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: intermediate_array
	243	INTEGER :: nx, ny, nz
	244
	245	nx = palm_intermediate % nx
	246	ny = palm_intermediate % ny
	247	nz = palm_intermediate % nz
	248	ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nz-1))
	249	intermediate_array(:,:,:) = 0.0_dp
	250
	251	! average input array to intermediate profile
	252	CALL average_2d(source_array, intermediate_array(0,0,:), imin, imax, jmin, jmax)
	253
	254	! vertically interpolate to ouput array
	255	CALL interpolate_1d(intermediate_array, profile_array, palm_grid)
	256
	257	DEALLOCATE(intermediate_array)
	258
	259	END SUBROUTINE average_profile
	260
	261
	262
	263	!-----------------------------------------------------------------------------!
	264	! Description:
	265	! -----------
	266	!> This routine computes the inverse Plate Carree projection, i.e. in projects
	267	!> Cartesian coordinates (x,y) onto a sphere. It returns the latitude and
	268	!> lngitude of a geographical system centered at x0 and y0.
	269	!-----------------------------------------------------------------------------!
	270	SUBROUTINE inv_plate_carree(x, y, x0, y0, r, lat, lon)
	271	REAL(dp), INTENT(IN) :: x(:), y(:), x0, y0, r
	272	REAL(dp), INTENT(OUT) :: lat(:), lon(:)
	273
	274	REAL(dp) :: ri
	275
	276	! TODO check dimensions of lat/lon and y/x match
	277
	278	ri = 1.0_dp / r
	279
	280	lat(:) = (y(:) - y0) * ri
	281	lon(:) = (x(:) - x0) * ri
	282	END SUBROUTINE
	283
	284
	285	!-----------------------------------------------------------------------------!
	286	! Description:
	287	! ------------
	288	!> Computes the reverse Plate-Carree projection of a x or y position on a
	289	!> Cartesian grid.
	290	!>
	291	!> Input parameters:
	292	!> -----------------
	293	!> xy : x or y coordinate of the Cartasian grid point [m].
	294	!>
	295	!> xy0 : x or y coordinate that coincides with the origin of the underlying
	296	!> sperical system (crossing point of the equator and prime meridian) [m].
	297	!>
	298	!> r : Radius of the of the underlying sphere, e.g. EARTH_RADIUS [m].
	299	!>
	300	!> Returns:
	301	!> --------
	302	!> project : Longitude (in case xy = x) or latitude (xy = y) of the given input
	303	!> coordinate xy.
	304	!------------------------------------------------------------------------------!
	305	ELEMENTAL REAL(dp) FUNCTION project(xy, xy0, r)
	306	REAL(dp), INTENT(IN) :: xy, xy0, r
	307	REAL(dp) :: ri
	308
	309	! If this elemental function is called with a large array as xy, it is
	310	! computationally more efficient to precompute the inverse radius and
	311	! then muliply.
	312	ri = 1.0_dp / r
	313
	314	project = (xy - xy0) * ri
	315
	316	END FUNCTION project
	317
	318
	319	REAL(dp) FUNCTION phic_to_phin(phi_c)
	320	REAL(dp), INTENT(IN) :: phi_c
	321
	322	phic_to_phin = 0.5_dp * PI - ABS(phi_c)
	323
	324	END FUNCTION phic_to_phin
	325
	326
	327	REAL(dp) FUNCTION lamc_to_lamn(phi_c, lam_c)
	328	REAL(dp), INTENT(IN) :: phi_c, lam_c
	329
	330	lamc_to_lamn = lam_c
	331	IF (phi_c > 0.0_dp) THEN
	332	lamc_to_lamn = lam_c - SIGN(PI, lam_c)
	333	END IF
	334
	335	END FUNCTION lamc_to_lamn
	336
	337
	338	REAL(dp) FUNCTION gamma_from_hemisphere(phi_cg, phi_ref)
	339	REAL(dp), INTENT(IN) :: phi_cg, phi_ref
	340	LOGICAL :: palm_centre_is_south_of_cosmo_origin
	341
	342	palm_centre_is_south_of_cosmo_origin = (phi_cg < phi_ref)
	343
	344	IF (palm_centre_is_south_of_cosmo_origin) THEN
	345	gamma_from_hemisphere = PI
	346	ELSE
	347	gamma_from_hemisphere = 0.0_dp
	348	END IF
	349	END FUNCTION gamma_from_hemisphere
	350
	351
	352	!------------------------------------------------------------------------------!
	353	! Description:
	354	! ------------
	355	!> Computes the geographical coordinates corresponding to the given rotated-pole
	356	!> coordinates.
	357	!>
	358	!> In INIFOR, this routine is used to convert coordinates between two
	359	!> rotated-pole systems: COSMO-DE's rotated-pole system, and one centred at the
	360	!> PALM-4U domain centre. In this case, the PALM-4U system is thought of as the
	361	!> rotated-pole system and the routine is used to rotate back to COSMO-DE's
	362	!> system which is thought of as the geographical one.
	363	!>
	364	!> Input parameters:
	365	!> -----------------
	366	!> phir(:), lamr(: ): latitudes and longitudes of the rotated-pole grid
	367	!>
	368	!> phip, lamp: latitude and longitude of the rotated north pole
	369	!>
	370	!> gam: "angle between the north poles. If [gam] is not present, the other
	371	!> system is the real geographical system." (original phiro2rot
	372	!> description)
	373	!>
	374	!> Output parameters:
	375	!> ------------------
	376	!> phi(:,:), lam(:,:): geographical latitudes and logitudes
	377	!------------------------------------------------------------------------------!
	378	SUBROUTINE rotate_to_cosmo(phir, lamr, phip, lamp, phi, lam, gam)
	379	REAL(dp), INTENT(IN) :: phir(0:), lamr(0:), phip, lamp, gam
	380	REAL(dp), INTENT(OUT) :: phi(0:,0:), lam(0:,0:)
	381
	382	INTEGER :: i, j
	383
	384	IF ( SIZE(phi, 1) .NE. SIZE(lam, 1) .OR. &
	385	SIZE(phi, 2) .NE. SIZE(lam, 2) ) THEN
	386	PRINT *, "inifor: rotate_to_cosmo: Dimensions of phi and lambda do not match. Dimensions are:"
	387	PRINT *, "inifor: rotate_to_cosmo: phi: ", SIZE(phi, 1), SIZE(phi, 2)
	388	PRINT *, "inifor: rotate_to_cosmo: lam: ", SIZE(lam, 1), SIZE(lam, 2)
	389	STOP
	390	END IF
	391
	392	IF ( SIZE(phir) .NE. SIZE(phi, 2) .OR. &
	393	SIZE(lamr) .NE. SIZE(phi, 1) ) THEN
	394	PRINT *, "inifor: rotate_to_cosmo: Dimensions of phir and lamr do not match. Dimensions are:"
	395	PRINT *, "inifor: rotate_to_cosmo: phir: ", SIZE(phir), SIZE(phi, 2)
	396	PRINT *, "inifor: rotate_to_cosmo: lamr: ", SIZE(lamr), SIZE(phi, 1)
	397	STOP
	398	END IF
	399
	400	DO j = 0, UBOUND(phir, 1)
	401	DO i = 0, UBOUND(lamr, 1)
	402
	403	phi(i,j) = phirot2phi(phir(j) * TO_DEGREES, &
	404	lamr(i) * TO_DEGREES, &
	405	phip * TO_DEGREES, &
	406	lamp * TO_DEGREES, &
	407	gam * TO_DEGREES) * TO_RADIANS
	408
	409	lam(i,j) = rlarot2rla(phir(j) * TO_DEGREES, &
	410	lamr(i) * TO_DEGREES, &
	411	phip * TO_DEGREES, &
	412	lamp * TO_DEGREES, &
	413	gam * TO_DEGREES) * TO_RADIANS
	414
	415	END DO
	416	END DO
	417
	418	END SUBROUTINE rotate_to_cosmo
[3182]	419
[2696]	420
	421
	422	!------------------------------------------------------------------------------!
	423	! Description:
	424	! ------------
	425	!> Compute indices of PALM-4U grid point neighbours in the target
	426	!> system (COSMO-DE) by rounding up and down. (i,j) are the indices of
	427	!> the PALM-4U grid and (ii(i,j,1-4), jj(i,j,1-4)) contain the indices
	428	!> of the its four neigbouring points in the COSMO-DE grid.
	429	!>
	430	!>
	431	!> COSMO-DE grid
	432	!> -------------
	433	!> jj, lat
[3182]	434	!> ^ j
	435	!> \| \ i
[2696]	436	!> jj(i,j,2/3) + ... 2 ---\--------/------ 3
	437	!> \| \| ^ \ / \|
	438	!> \| \| \|wp \ / \|
	439	!> \| \| v \ / \|
	440	!> latpos + ............ o/ (i,j) \|
	441	!> \| \| : \|
	442	!> \| \| :<----wl---->\|
	443	!> jj(i,j,1/4) + ... 1 -------:----------- 4
	444	!> \| : : :
	445	!> \| : : :
	446	!> \| : lonpos :
	447	!> L-----+--------+------------+------> ii, lon
	448	!> ii(i,j,1/2) ii(i,j,3/4)
	449	!>
	450	!>
	451	!> Input parameters:
	452	!> -----------------
	453	!> source_lat, source_lon : (rotated-pole) coordinates of the source grid (e.g.
	454	!> COSMO-DE)
	455	!>
	456	!> source_dxi, source_dyi : inverse grid spacings of the source grid.
	457	!>
	458	!> target_lat, target_lon : (rotated-pole) coordinates of the target grid (e.g.
	459	!> COSMO-DE)
	460	!>
	461	!> Output parameters:
	462	!> ------------------
	463	!> palm_ii, palm_jj : x and y index arrays of horizontal neighbour columns
	464	!>
	465	!------------------------------------------------------------------------------!
[3182]	466	SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon, &
	467	palm_clat, palm_clon, &
	468	palm_ii, palm_jj)
[2696]	469
	470	REAL(dp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon
	471	REAL(dp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon
[3182]	472	REAL(dp) :: cosmo_dxi, cosmo_dyi
[2696]	473	INTEGER, DIMENSION(0:,0:,1:), INTENT(OUT) :: palm_ii, palm_jj
	474
	475	REAL(dp) :: lonpos, latpos, lon0, lat0
	476	INTEGER :: i, j
	477
	478	lon0 = cosmo_lon(0)
	479	lat0 = cosmo_lat(0)
[3182]	480	cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
	481	cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
[2696]	482
	483	DO j = 0, UBOUND(palm_clon, 2)!palm_grid % ny
	484	DO i = 0, UBOUND(palm_clon, 1)!palm_grid % nx
	485	! Compute the floating point index corrseponding to PALM-4U grid point
	486	! location along target grid (COSMO-DE) axes.
	487	lonpos = (palm_clon(i,j) - lon0) * cosmo_dxi
	488	latpos = (palm_clat(i,j) - lat0) * cosmo_dyi
	489
	490	IF (lonpos < 0.0 .OR. latpos < 0.0) THEN
	491	PRINT *, " Error while finding neighbours: lonpos or latpos out of bounds!"
	492	PRINT *, " (i,j) = (", i, ",",j,")"
	493	PRINT , " lonpos ", lonposTO_DEGREES, ", latpos ", latpos*TO_DEGREES
	494	PRINT , " lon0 ", lon0 TO_DEGREES, ", lat0 ", lat0*TO_DEGREES
	495	PRINT , " PALM lon ", palm_clon(i,j)TO_DEGREES, ", PALM lat ",palm_clat(i,j)*TO_DEGREES
	496	STOP
	497	END IF
	498
	499	palm_ii(i,j,1) = FLOOR(lonpos)
	500	palm_ii(i,j,2) = FLOOR(lonpos)
	501	palm_ii(i,j,3) = CEILING(lonpos)
	502	palm_ii(i,j,4) = CEILING(lonpos)
	503
	504	palm_jj(i,j,1) = FLOOR(latpos)
	505	palm_jj(i,j,2) = CEILING(latpos)
	506	palm_jj(i,j,3) = CEILING(latpos)
	507	palm_jj(i,j,4) = FLOOR(latpos)
	508	END DO
	509	END DO
	510
	511	END SUBROUTINE find_horizontal_neighbours
	512
	513
	514	SUBROUTINE find_vertical_neighbours_and_weights(palm_grid, palm_intermediate)
	515	TYPE(grid_definition), INTENT(INOUT) :: palm_grid
	516	TYPE(grid_definition), INTENT(IN) :: palm_intermediate
	517
[3182]	518	INTEGER :: i, j, k, nx, ny, nz, nlev, k_intermediate
[2696]	519	LOGICAL :: point_is_below_grid, point_is_above_grid, &
	520	point_is_in_current_cell
	521	REAL(dp) :: current_height, column_base, column_top, h_top, h_bottom, &
	522	weight
	523
	524	nx = palm_grid % nx
	525	ny = palm_grid % ny
	526	nz = palm_grid % nz
	527	nlev = palm_intermediate % nz
	528
	529	! in each column of the fine grid, find vertical neighbours of every cell
	530	DO i = 0, nx
	531	DO j = 0, ny
	532
[3182]	533	k_intermediate = 0
[2696]	534
	535	column_base = palm_intermediate % h(i,j,0)
	536	column_top = palm_intermediate % h(i,j,nlev)
	537
	538	! scan through palm_grid column until and set neighbour indices in
	539	! case current_height is either below column_base, in the current
	540	! cell, or above column_top. Keep increasing current cell index until
	541	! the current cell overlaps with the current_height.
[3182]	542	DO k = 1, nz
[2696]	543
	544	! Memorize the top and bottom boundaries of the coarse cell and the
	545	! current height within it
	546	current_height = palm_grid % z(k) + palm_grid % z0
[3182]	547	h_top = palm_intermediate % h(i,j,k_intermediate+1)
	548	h_bottom = palm_intermediate % h(i,j,k_intermediate)
[2696]	549
	550	point_is_above_grid = (current_height > column_top) !22000m, very unlikely
	551	point_is_below_grid = (current_height < column_base)
	552
	553	point_is_in_current_cell = ( &
	554	current_height >= h_bottom .AND. &
	555	current_height < h_top &
	556	)
	557
	558	! set default weights
	559	palm_grid % w_verti(i,j,k,1:2) = 0.0_dp
	560
	561	IF (point_is_above_grid) THEN
	562
	563	palm_grid % kk(i,j,k,1:2) = nlev
	564	palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp
	565
[3182]	566	message = "PALM-4U grid extends above COSMO-DE model top."
	567	CALL abort('find_vertical_neighbours_and_weights', message)
	568
[2696]	569	ELSE IF (point_is_below_grid) THEN
	570
	571	palm_grid % kk(i,j,k,1:2) = 0
	572	palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp
	573
	574	ELSE
	575	! cycle through intermediate levels until current
	576	! intermediate-grid cell overlaps with current_height
[3182]	577	DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
	578	k_intermediate = k_intermediate + 1
[2696]	579
[3182]	580	h_top = palm_intermediate % h(i,j,k_intermediate+1)
	581	h_bottom = palm_intermediate % h(i,j,k_intermediate)
[2696]	582	point_is_in_current_cell = ( &
	583	current_height >= h_bottom .AND. &
	584	current_height < h_top &
	585	)
	586	END DO
	587
[3182]	588	! k_intermediate = 48 indicates the last section (indices 48 and 49), i.e.
	589	! k_intermediate = 49 is not the beginning of a valid cell.
	590	IF (k_intermediate > nlev-1) THEN
	591	message = "Index " // TRIM(str(k_intermediate)) // &
	592	" is above intermediate grid range."
[2696]	593	CALL abort('find_vertical_neighbours', message)
	594	END IF
	595
[3182]	596	palm_grid % kk(i,j,k,1) = k_intermediate
	597	palm_grid % kk(i,j,k,2) = k_intermediate + 1
[2696]	598
	599	! copmute vertical weights
	600	weight = (h_top - current_height) / (h_top - h_bottom)
	601	palm_grid % w_verti(i,j,k,1) = weight
	602	palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight
	603	END IF
	604
	605	END DO
	606
	607	END DO
	608	END DO
	609
	610	END SUBROUTINE find_vertical_neighbours_and_weights
	611
	612	!------------------------------------------------------------------------------!
	613	! Description:
	614	! ------------
	615	!> Compute the four weights for horizontal bilinear interpolation given the
	616	!> coordinates clon(i,j) clat(i,j) of the PALM-4U grid in the COSMO-DE
	617	!> rotated-pole grid and the neightbour indices ii(i,j,1-4) and jj(i,j,1-4).
	618	!>
	619	!> Input parameters:
	620	!> -----------------
	621	!> palm_grid % clon : longitudes of PALM-4U scalars (cell centres) in COSMO-DE's rotated-pole grid [rad]
	622	!>
	623	!> palm_grid % clat : latitudes of PALM-4U cell centres in COSMO-DE's rotated-pole grid [rad]
	624	!>
	625	!> cosmo_grid % lon : rotated-pole longitudes of scalars (cell centres) of the COSMO-DE grid [rad]
	626	!>
	627	!> cosmo_grid % lat : rotated-pole latitudes of scalars (cell centers) of the COSMO-DE grid [rad]
	628	!>
	629	!> cosmo_grid % dxi : inverse grid spacing in the first dimension [m^-1]
	630	!>
	631	!> cosmo_grid % dyi : inverse grid spacing in the second dimension [m^-1]
	632	!>
	633	!> Output parameters:
	634	!> ------------------
	635	!> palm_grid % w_horiz(:,:,1-4) : weights for bilinear horizontal interpolation
	636	!
	637	! COSMO-DE grid
	638	! -------------
	639	! jj, lat
	640	! ^ j
	641	! \| \ i
	642	! jj(i,j,2/3) + ... 2 ---\--------/------ 3
	643	! \| \| ^ \ / \|
	644	! \| \| \|wp \ / \|
	645	! \| \| v \ / \|
	646	! latpos + ............ o/ (i,j) \|
	647	! \| \| : \|
	648	! \| \| :<----wl---->\|
	649	! jj(i,j,1/4) + ... 1 -------:----------- 4
	650	! \| : : :
	651	! \| : : :
	652	! \| : lonpos :
	653	! L-----+--------+------------+------> ii, lon
	654	! ii(i,j,1/2) ii(i,j,3/4)
	655	!
	656	SUBROUTINE compute_horizontal_interp_weights(cosmo_lat, cosmo_lon, &
[3182]	657	palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)
[2696]	658
	659	REAL(dp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon
[3182]	660	REAL(dp) :: cosmo_dxi, cosmo_dyi
[2696]	661	REAL(dp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon
	662	INTEGER, DIMENSION(0:,0:,1:), INTENT(IN) :: palm_ii, palm_jj
	663
	664	REAL(dp), DIMENSION(0:,0:,1:), INTENT(OUT) :: palm_w_horiz
	665
	666	REAL(dp) :: wl, wp
	667	INTEGER :: i, j
	668
[3182]	669	cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
	670	cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
	671
[2696]	672	DO j = 0, UBOUND(palm_clon, 2)
	673	DO i = 0, UBOUND(palm_clon, 1)
	674
	675	! weight in lambda direction
	676	wl = ( cosmo_lon(palm_ii(i,j,4)) - palm_clon(i,j) ) * cosmo_dxi
	677
	678	! weight in phi direction
	679	wp = ( cosmo_lat(palm_jj(i,j,2)) - palm_clat(i,j) ) * cosmo_dyi
	680
	681	IF (wl > 1.0_dp .OR. wl < 0.0_dp) THEN
	682	message = "Horizontal weight wl = " // TRIM(real_to_str(wl)) // &
	683	" is out bounds."
	684	CALL abort('compute_horizontal_interp_weights', message)
	685	END IF
	686	IF (wp > 1.0_dp .OR. wp < 0.0_dp) THEN
	687	message = "Horizontal weight wp = " // TRIM(real_to_str(wp)) // &
	688	" is out bounds."
	689	CALL abort('compute_horizontal_interp_weights', message)
	690	END IF
	691
	692	palm_w_horiz(i,j,1) = wl * wp
	693	palm_w_horiz(i,j,2) = wl * (1.0_dp - wp)
	694	palm_w_horiz(i,j,3) = (1.0_dp - wl) * (1.0_dp - wp)
	695	palm_w_horiz(i,j,4) = 1.0_dp - SUM( palm_w_horiz(i,j,1:3) )
	696
	697	END DO
	698	END DO
	699
	700	END SUBROUTINE compute_horizontal_interp_weights
	701
	702
	703	!------------------------------------------------------------------------------!
	704	! Description:
	705	! ------------
	706	!> Interpolates u and v components of velocities located at cell faces to the
	707	!> cell centres by averaging neighbouring values.
	708	!>
	709	!> This routine is designed to be used with COSMO-DE arrays where there are the
	710	!> same number of grid points for scalars (centres) and velocities (faces). In
	711	!> COSMO-DE the velocity points are staggared one half grid spaceing up-grid
	712	!> which means the first centre point has to be omitted and is set to zero.
	713	SUBROUTINE centre_velocities(u_face, v_face, u_centre, v_centre)
	714	REAL(dp), DIMENSION(0:,0:,0:), INTENT(IN) :: u_face, v_face
	715	REAL(dp), DIMENSION(0:,0:,0:), INTENT(OUT) :: u_centre, v_centre
	716	INTEGER :: nx, ny
	717
	718	nx = UBOUND(u_face, 1)
	719	ny = UBOUND(u_face, 2)
	720
	721	u_centre(0,:,:) = 0.0_dp
	722	u_centre(1:,:,:) = 0.5_dp * ( u_face(0:nx-1,:,:) + u_face(1:,:,:) )
	723
	724	v_centre(:,0,:) = 0.0_dp
	725	v_centre(:,1:,:) = 0.5_dp * ( v_face(:,0:ny-1,:) + v_face(:,1:,:) )
	726	END SUBROUTINE centre_velocities
	727
	728
	729	FUNCTION phirot2phi (phirot, rlarot, polphi, pollam, polgam)
	730
	731	REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole
	732	REAL(dp), INTENT (IN) :: pollam !< longitude of the rotated north pole
	733	REAL(dp), INTENT (IN) :: phirot !< latitude in the rotated system
	734	REAL(dp), INTENT (IN) :: rlarot !< longitude in the rotated system
	735	REAL(dp), INTENT (IN) :: polgam !< angle between the north poles of the systems
	736
	737	REAL(dp) :: phirot2phi !< latitude in the geographical system
	738
	739	REAL(dp) :: zsinpol, zcospol, zphis, zrlas, zarg, zgam
	740
	741	zsinpol = SIN(polphi * TO_RADIANS)
	742	zcospol = COS(polphi * TO_RADIANS)
	743	zphis = phirot * TO_RADIANS
	744
	745	IF (rlarot > 180.0_dp) THEN
	746	zrlas = rlarot - 360.0_dp
	747	ELSE
	748	zrlas = rlarot
	749	END IF
	750	zrlas = zrlas * TO_RADIANS
	751
	752	IF (polgam /= 0.0_dp) THEN
	753	zgam = polgam * TO_RADIANS
	754	zarg = zsinpol * SIN (zphis) + &
	755	zcospol * COS(zphis) * ( COS(zrlas) * COS(zgam) - &
	756	SIN(zgam) * SIN(zrlas) )
	757	ELSE
	758	zarg = zcospol * COS (zphis) * COS (zrlas) + zsinpol * SIN (zphis)
	759	END IF
	760
	761	phirot2phi = ASIN (zarg) * TO_DEGREES
	762
	763	END FUNCTION phirot2phi
	764
	765
	766	FUNCTION phi2phirot (phi, rla, polphi, pollam)
	767
	768	REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole
	769	REAL(dp), INTENT (IN) :: pollam !< longitude of the rotated north pole
	770	REAL(dp), INTENT (IN) :: phi !< latitude in the geographical system
	771	REAL(dp), INTENT (IN) :: rla !< longitude in the geographical system
	772
	773	REAL(dp) :: phi2phirot !< longitude in the rotated system
	774
	775	REAL(dp) :: zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1
	776
	777	zsinpol = SIN(polphi * TO_RADIANS)
	778	zcospol = COS(polphi * TO_RADIANS)
	779	zlampol = pollam * TO_RADIANS
	780	zphi = phi * TO_RADIANS
	781
	782	IF (rla > 180.0_dp) THEN
	783	zrla1 = rla - 360.0_dp
	784	ELSE
	785	zrla1 = rla
	786	END IF
	787	zrla = zrla1 * TO_RADIANS
	788
	789	zarg1 = SIN(zphi) * zsinpol
	790	zarg2 = COS(zphi) * zcospol * COS(zrla - zlampol)
	791
	792	phi2phirot = ASIN(zarg1 + zarg2) * TO_DEGREES
	793
	794	END FUNCTION phi2phirot
	795
	796
	797	FUNCTION rlarot2rla(phirot, rlarot, polphi, pollam, polgam)
	798
	799	REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole
	800	REAL(dp), INTENT (IN) :: pollam !< longitude of the rotated north pole
	801	REAL(dp), INTENT (IN) :: phirot !< latitude in the rotated system
	802	REAL(dp), INTENT (IN) :: rlarot !< longitude in the rotated system
	803	REAL(dp), INTENT (IN) :: polgam !< angle between the north poles of the systems
	804
	805	REAL(dp) :: rlarot2rla !< latitude in the geographical system
	806
	807	REAL(dp) :: zsinpol, zcospol, zlampol, zphis, zrlas, zarg1, zarg2, zgam
	808
	809	zsinpol = SIN(TO_RADIANS * polphi)
	810	zcospol = COS(TO_RADIANS * polphi)
	811	zlampol = TO_RADIANS * pollam
	812	zphis = TO_RADIANS * phirot
	813
	814	IF (rlarot > 180.0_dp) THEN
	815	zrlas = rlarot - 360.0_dp
	816	ELSE
	817	zrlas = rlarot
	818	END IF
	819	zrlas = TO_RADIANS * zrlas
	820
	821	IF (polgam /= 0.0_dp) THEN
	822	zgam = TO_RADIANS * polgam
	823	zarg1 = SIN(zlampol) * (zcospol * SIN(zphis) - zsinpolCOS(zphis) &
	824	(COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) - &
	825	COS(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) + &
	826	COS(zrlas) * SIN(zgam) )
	827
	828	zarg2 = COS (zlampol) * (zcospol * SIN(zphis) - zsinpolCOS(zphis) &
	829	(COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) + &
	830	SIN(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) + &
	831	COS(zrlas) * SIN(zgam) )
	832	ELSE
	833	zarg1 = SIN (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis) + &
	834	zcospol * SIN(zphis)) - &
	835	COS (zlampol) * SIN(zrlas) * COS(zphis)
	836	zarg2 = COS (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis) + &
	837	zcospol * SIN(zphis)) + &
	838	SIN (zlampol) * SIN(zrlas) * COS(zphis)
	839	END IF
	840
	841	IF (zarg2 == 0.0_dp) zarg2 = 1.0E-20_dp
	842
	843	rlarot2rla = ATAN2(zarg1,zarg2) * TO_DEGREES
	844
	845	END FUNCTION rlarot2rla
	846
	847
	848	FUNCTION rla2rlarot ( phi, rla, polphi, pollam, polgam )
	849
	850	REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole
	851	REAL(dp), INTENT (IN) :: pollam !< longitude of the rotated north pole
	852	REAL(dp), INTENT (IN) :: phi !< latitude in geographical system
	853	REAL(dp), INTENT (IN) :: rla !< longitude in geographical system
	854	REAL(dp), INTENT (IN) :: polgam !< angle between the north poles of the systems
	855
	856	REAL (KIND=dp) :: rla2rlarot !< latitude in the the rotated system
	857
	858	REAL (KIND=dp) :: zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1
	859
	860	zsinpol = SIN(polphi * TO_RADIANS)
	861	zcospol = COS(polphi * TO_RADIANS)
	862	zlampol = pollam * TO_RADIANS
	863	zphi = phi * TO_RADIANS
	864
	865	IF (rla > 180.0_dp) THEN
	866	zrla1 = rla - 360.0_dp
	867	ELSE
	868	zrla1 = rla
	869	END IF
	870	zrla = zrla1 * TO_RADIANS
	871
	872	zarg1 = - SIN (zrla-zlampol) * COS(zphi)
	873	zarg2 = - zsinpol * COS(zphi) * COS(zrla-zlampol) + zcospol * SIN(zphi)
	874
	875	IF (zarg2 == 0.0_dp) zarg2 = 1.0E-20_dp
	876
	877	rla2rlarot = ATAN2 (zarg1,zarg2) * TO_DEGREES
	878
	879	IF (polgam /= 0.0_dp ) THEN
	880	rla2rlarot = polgam + rla2rlarot
	881	IF (rla2rlarot > 180._dp) rla2rlarot = rla2rlarot - 360.0_dp
	882	END IF
	883
	884	END FUNCTION rla2rlarot
	885
	886
	887	SUBROUTINE uv2uvrot(u, v, rlat, rlon, pollat, pollon, urot, vrot)
	888
	889	REAL(dp), INTENT (IN) :: u, v !< wind components in the true geographical system
	890	REAL(dp), INTENT (IN) :: rlat, rlon !< coordinates in the true geographical system
	891	REAL(dp), INTENT (IN) :: pollat, pollon !< latitude and longitude of the north pole of the rotated grid
	892
	893	REAL(dp), INTENT (OUT) :: urot, vrot !< wind components in the rotated grid
	894
	895	REAL (dp) :: zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm
	896
	897	zsinpol = SIN(pollat * TO_RADIANS)
	898	zcospol = COS(pollat * TO_RADIANS)
	899	zlonp = (pollon-rlon) * TO_RADIANS
	900	zlat = rlat * TO_RADIANS
	901
	902	zarg1 = zcospol * SIN(zlonp)
	903	zarg2 = zsinpol * COS(zlat) - zcospol * SIN(zlat) * COS(zlonp)
	904	znorm = 1.0_dp / SQRT(zarg1zarg1 + zarg2zarg2)
	905
	906	urot = u * zarg2 * znorm - v * zarg1 * znorm
	907	vrot = u * zarg1 * znorm + v * zarg2 * znorm
	908
	909	END SUBROUTINE uv2uvrot
	910
	911
	912	SUBROUTINE uvrot2uv (urot, vrot, rlat, rlon, pollat, pollon, u, v)
	913
	914	REAL(dp), INTENT(IN) :: urot, vrot !< wind components in the rotated grid
	915	REAL(dp), INTENT(IN) :: rlat, rlon !< latitude and longitude in the true geographical system
	916	REAL(dp), INTENT(IN) :: pollat, pollon !< latitude and longitude of the north pole of the rotated grid
	917
	918	REAL(dp), INTENT(OUT) :: u, v !< wind components in the true geographical system
	919
	920	REAL(dp) :: zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm
	921
	922	zsinpol = SIN(pollat * TO_RADIANS)
	923	zcospol = COS(pollat * TO_RADIANS)
	924	zlonp = (pollon-rlon) * TO_RADIANS
	925	zlat = rlat * TO_RADIANS
	926
	927	zarg1 = zcospol * SIN(zlonp)
	928	zarg2 = zsinpol * COS(zlat) - zcospol * SIN(zlat) * COS(zlonp)
	929	znorm = 1.0_dp / SQRT(zarg1zarg1 + zarg2zarg2)
	930
	931	u = urot * zarg2 * znorm + vrot * zarg1 * znorm
	932	v = - urot * zarg1 * znorm + vrot * zarg2 * znorm
	933
	934	END SUBROUTINE uvrot2uv
	935
	936	END MODULE
	937

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |