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

Last change on this file since 3716 was 3716, checked in by eckhard, 5 years ago

inifor: bugfix: removed dependency on soilmoisture input files; added netcdf preprocessor flag

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