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

Last change on this file since 3614 was 3614, checked in by raasch, 5 years ago

unused variables removed, abort renamed inifor_abort to avoid intrinsic problem in Fortran

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