Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |