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

Last change on this file since 3615 was 3615, checked in by raasch, 3 years ago

bugfix for last commit: abort replaced by inifor_abort

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