Home

Context Navigation

transform.f90 @ 3421

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

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

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

Context Navigation

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

Download in other formats: