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

Last change on this file since 3648 was 3618, checked in by eckhard, 6 years ago

inifor: Prefixed all INIFOR modules with inifor_ and removed unused variables

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