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

Last change on this file since 3678 was 3678, checked in by eckhard, 3 years ago

inifor: bugfix: avoid empty averaging regions, check if all input files are present

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